text:
# Hubble Space Telescope Expansion Parallaxes of the Planetary Nebulae NGC 6578, NGC 6884, NGC 6891, and IC 2448Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555.

## 1 Introduction

Few planetary nebulae (PNe) in our galaxy have accurately determined distances. This has hampered our ability to use PNe as distance indicators in other galaxies, as well as limiting the accuracy of derived PN properties which are strongly dependent on distance, such as luminosity or nebular mass. Even the spatial distribution of PNe in the Milky Way is poorly determined, leading to uncertainties in the total number of PNe, and the birth rate. Because PN progenitors are sources of carbon and nitrogen, the PN population strongly influences the chemical evolution of the galaxy as a whole (see, for example, Martins & Viegas 2000). In particular, these progenitors contribute to the chemical evolution of carbon in the interstellar medium (e.g., Iben 1985; Palla et al. 2000).

Most PN distances have been derived statistically, and are accurate only when taken as a whole (not for individual members of the population). Prior to 1995, statistical distance scales relied on an assumption of a universal physical property such as the ionized nebular mass, luminosity of the central star, or uniform extinction along the line of sight. These assumptions often introduce large errors. For example, PNe are formed from progenitor stars with masses ranging over an order of magnitude, all of which shed enough mass to drop to between $``$0.5 and $`1M_{}`$ during the PN phase. These statistical scales generally gave results accurate to a factor of two at best for individual targets (see, for example, Shklovsky 1956; O’Dell 1962; Daub 1982; Maciel & Pottasch 1980; Cahn, Kaler, & Stanghellini 1992). See Terzian (1993) for a review of the difficulties inherent in finding distances to PNe.

More recently, Zhang (1995) devised two additional statistical methods to determine distances to PNe, based on empirically determined relationships between the ionized mass and the radius of the object (mass-radius relation), and between the radio continuum surface brightness and the radius (surface brightness-radius relation). Using both of these relations, he claims an accuracy of $`35\%50\%`$. These distances are still not sufficiently accurate to enable a solution to several different PN problems. In part, the poor accuracy of the statistical distances is a result of calibration to a small sample of nearby PNe with distances that have been determined by other means. Most recently, several efforts have been underway to improve this nearby sample.

Trigonometric parallax has been used in recent years to find the distances to PNe, as CCDs have improved the astrometric measurements. Harris et al. (1997) report distances to 11 PNe accurate to 20%, and a further 5 with distances accurate to 50%, all from ground-based trigonometric parallax. The Hipparcos satellite measured parallax distances to a small sample of PNe, most to less than 50% accuracy (Acker et al. 1998). Pottasch & Acker (1998) compared these distances to previously determined spectroscopic distances to three PNe, and found that the parallax distances were all much smaller than the spectroscopic distances. Gutierrez-Moreno et al. (1999) found ground-based trigonometric parallaxes for three objects, and determined that there was no correlation between the previously determined statistical distances and the distances that they find from this fundamentally direct method. They do find, however, that distances determined from other methods, such as expansion parallax or sodium absorption, are highly correlated with their results. They interpret this as a fundamental failing of the statistical methods.

Ciardullo et al. (1999) have used HST to discover close companions of the central stars of PNe. They used main-sequence photometric parallax to derive reliable distances to nine objects with companions. Distances to a further three objects are given in that paper, but the association of the binary pair is less certain. Bond & Ciardullo (1999) also reported a ground-based distance to NGC 246 based on photometry of the resolved binary companion.

Expansion-parallax distances have been determined using the Very Large Array (VLA) for eight PNe (Masson 1989a,b; Gomez, Rodriguez, & Moran 1993; Hajian, Terzian, & Bignell 1993; Hajian, Terzian, & Bignell 1995; Hajian & Terzian 1996). These distances are precise to roughly 10-20% in some cases, provided that the assumptions of elliptical morphology and constant expansion velocity are correct. Lower limits on about six more objects have been determined using this method.

The high spatial resolution offered by the Hubble Space Telescope (HST) makes it potentially very useful for determining expansion parallaxes of PNe. Prior to the work described in this paper, Reed et al. (1999, hereafter Paper I) have used HST to determine the distance to one PN, NGC 6543, using the expansion-parallax method .

The expansion-parallax method relies on the assumption that the expanding nebula is a prolate ellipsoid. Then the minor axis expansion velocity in the plane of the sky ($`V_m`$) is related to the line of sight velocity ($`V`$) at the position of the central star. In the case of complex nebulae, such as point-symmetric, or extremely bipolar nebulae, the method fails. Fortunately, if a nebula is close enough that $`V_m`$ is observable over a few years, it is also close enough that the morphology may be accurately determined.

We also assume, in this method, that the spectroscopy and the imaging are sampling the same physical regimes within the object. It is possible that the angular expansion shows the motion of the ionization front, while the radial velocity measures the bulk motion of the gas. This will lead to an underestimate of the distance to the nebula.

A second assumption is that the nebula has been expanding at the same velocity for the elapsed time over which the parallax is being measured. The VLA and HST have sufficient angular resolution and fidelity to measure an expansion in just four years. It is unlikely that the expansion speed of the nebular gas will be dramatically altered in this short period of time.

The distance is given by the small-angle formula:

$$D=\frac{V_m}{\dot{\theta }},$$
(1)
where $`V_m`$ is the expansion velocity of the emitting material on the minor axis relative to the central star, and $`\dot{\theta }`$ is the angular expansion of the radius of the nebula as measured on the sky. Converting to convenient units gives:

$$D_{pc}=211\frac{V_m(km/s)}{\dot{\theta }(mas/yr)}.$$
(2)
Expansion velocities are derived from Doppler shifts along the line of sight to the central star, corrected for geometric projection effects.

We determined the angular expansion of the nebulae from the $`V`$ band HST images with two techniques: the gradient method and the magnification method, both described below. In this paper, we report the direct measurement of expansion parallax distances to three nebulae. §2 describes the observations and the data reduction, including the alignment of the HST images. §3 describes the three dimensional modeling technique used to determine the ratio of the line of sight velocity to the velocity along the minor axis. §4 describes in more detail the methods for finding the angular expansion of the nebula, and discusses the pros and cons of each method. In §5-8, we discuss the results and derived properties for each of the four PNe. Finally, in §9, we discuss possible sources of systematic errors, and future plans for the application of the methodologies described here.

## 2 Observations

### 2.1 Imaging with HST: Wide Field Planetary Camera 2

This paper describes results specifically for NGC 6578, NGC 6884, NGC 6891 and IC 2448. The first-epoch observations of these objects were obtained in a $`V`$-band (F555W filter) Wide Field Planetary Camera (WFPC2) snapshot survey by Bond and R. Ciardullo (Cycle 5, GO program 6119). The target objects were placed in the center of the planetary camera (PC) chip. Because this survey was optimized for observations of central stars, it was carried out in the broad $`V`$ band, rather than in a nebular emission line. Thus, comparison to the radial velocities obtained from \[O III\] spectroscopic data must be done with some care. However, since \[O III\] is almost always the main contributor to the $`V`$-band flux from PNe, comparing $`V`$-band imaging with \[O III\] spectroscopy will not be a problem. Other possible contributions to the $`V`$ band come from H$`\beta `$, He I (4921, 5015, 5876, and weaker lines), and usually negligible contributions from various Fe and N lines, as well as the nebular continuum (which will contribute more to the $`V`$-band image than to the \[O III\], since the spectral window is wider).

¿From the survey of Bond and Ciardullo, we chose the 24 most spherically symmetric objects, and reobserved many of them with WFPC2 on HST in Cycles 8 and 9. NGC 6578, NGC 6884, NGC 6891 and IC 2448 are the first four objects with second-epoch images. In the second epoch, we obtained \[O III\] (F502N) and \[N II\] (F658N) data in addition to new $`V`$-band images. In some cases, the \[N II\] images allow us to investigate in more detail complementary nebular volumes, and speculate about the nebular evolution. Table 1 summarizes the imaging observations.

#### 2.1.1 Imaging Calibrations

Cosmic rays were removed using the IRAF task ‘crrej’. Correction for optical camera distortions was performed with the IRAF/STSDAS task ‘drizzle’ (Fruchter & Hook 2002) with the Trauger coefficients. After these corrections, the chief problem was aligning images from different epochs properly. Despite our best efforts to repeat the first-epoch observations exactly, the images from the two epochs were not quite perfectly aligned. These small errors could be due to a number of effects, such as proper motions of guidestars, or drifts in the instruments. Translations of fewer than five pixels in each direction and rotations of a few tenths of a degree were required to register the two epochs in each case. Such translations aren’t important unless the optical distortion corrections are in error; even so, they are not likely to produce a false positive expansion for the entire nebula.

Using IDL, we made a first order translation correction by aligning the images relative to a “pivot” star. We found the centroids, and determined the offsets between epochs in the $`x`$-direction ($`\delta x`$) and $`y`$-direction ($`\delta y`$) by subtracting these coordinates.

We then chose three stars in each image, found the centroids in the two epochs (\[$`x_1`$,$`y_1`$\] and \[$`x_2`$,$`y_2`$\]), and solved for the rotation required to align the first image with the second,

$$\theta _{rot}=\mathrm{arctan}[(y_1\delta y)/(x_1\delta x)]\mathrm{arctan}(y_2/x_2).$$
(3)
We compared the three values of $`\theta _{rot}`$ from the three stars, to check for discrepancies due to proper motion, geometric distortions or other complicating factors. If the agreement was good, the three results were averaged to determine the rotation between the two images. If the agreement was bad, three new stars (or a new pivot) were chosen, and the process was repeated.

Finally, the translation of the images relative to the central star of the nebula was determined. For the case of an overexposed central star, we located the star by fitting lines to each of the diffraction spikes (which make an “X” on the image), rather than just finding the centroid.

To apply the transformations, the image from the first epoch was magnified by a factor of ten, the translations were applied, and the image was de-magnified to its original size. (IDL allows only shifts of integer pixel size, and uses a bilinear interpolation when magnifying, and neighborhood averaging when de-magnifying.) The second-epoch image was magnified, rotated, and de-magnified. Then the two images were subtracted and the alignment visually inspected. The point-spread function is rotated, and so the stars never subtract completely. We used the nebula itself to check the alignment, looking for patterns of offsets which could result from a shift between one epoch and the next. This could result from, for example, proper motions of the background stars. If necessary, one image was “nudged” relative to the other, with these small translations further improving the alignment.

To check that we are not introducing significant errors in this process, we applied two tests. First, we transformed the final images back and subtracted from the originals. The residuals were small in all cases (less than 4% when integrated across the entire nebula), and appeared mainly in places where the nebula was faint, and close to the level of the noise in the images. Second, we applied the rotation to the first-epoch image, and the translation to the second-epoch image, and subtracted the resulting difference map from the difference map derived above. Again, there was good agreement. This method of aligning the images results in only small artifacts, due mainly to the interpolations carried out when regridding the data near under-sampled features.

### 2.2 Spectroscopy

High spectral resolution echelle data were obtained from the four-meter telescopes at Kitt Peak and CTIO. In both cases, the cross-disperser was replaced with a mirror, so that only one order, centered around \[O III\], was obtained. The instrument resolution was 0.06 Å/pixel (3.7 km/s), and a narrow-band interference filter was used ($`70`$ Å). Errors in the expansion velocity obtained from \[O III\], a strong line in all cases, are of order 10%. The spectroscopic observations are summarized in Table 2.

#### 2.2.1 Spectroscopy Calibrations

The echelle spectra were bias-subtracted and flat-fielded. Thorium-Argon lamp spectra were used to solve for geometric distortions by iteratively using the IRAF tasks identify, re-identify, fitcoords and transform. In some instances, we used sky lines in the images to test the accuracy of these corrections. There were no significant deviations at the 0.1 pixel level. These lamp spectra were also used to calibrate the wavelength scale.

The expansion velocity along the line of sight to the central star was found from these data by fitting Gaussians to the double-peaked velocity profile, subtracting the velocities of the peaks, and dividing by two.

## 3 Modeling

The nebula is assumed to be a pinched oblate ellipsoid inclined to the line of sight, (see Figure 9). It is assumed to be expanding ballistically, so that the ratio of the axes is equal to the ratio of the velocities along those axes:

$$r_p/r_m=v_p/v_m,$$
(4)
where $`r_p`$ and $`v_p`$ are the radius and velocity along the long (polar) axis, and $`r_m`$ and $`v_m`$ are the radius and velocity along the short (minor) axis. Then, the ratio of the line of sight velocity, $`V`$, to the true minor axis velocity, $`v_m`$, is given by

$$\alpha =\frac{V}{v_m}=1+(r_p/r_m1)|\mathrm{cos}^n\varphi |,$$
(5)
where $`\varphi `$ is the azimuthal angle, and $`n`$ is functionally a shape parameter, defining how sharply the ‘waist’ is pinched.

We derive $`n`$, $`r_p/r_m`$ and the inclination angle from code written in IDL that simulates a ballistically expanding pinched-waist elliptical nebula. Once this nebula is created, a “slit” is placed across it, and a velocity spectrum derived, similar to that observed in the echelle. Both the projected nebula and the spectrum are overplotted on the observed data, and the parameters are varied manually until the shapes of the data are recovered. Because we have two sets of echelle data, one along the short axis, and one along the long axis, the parameters are uniquely determined once a set of values is found which matches the image and both spectra. The fits to each nebula are shown in Figures 9-9, and the derived parameters are given in Table 3. These derived parameters are used to determine the true velocity along the minor axis from the line of sight velocity:

$$v_m=\frac{V}{\alpha }.$$
(6)
## 4 Determining the Angular Expansion

### 4.1 The Gradient Method

In the simplest case, when a nebula expands self-similarly, a difference image (epoch 2 $``$ epoch 1) shows a bright outline surrounding a dark outline. When the expansion is resolved, features may expand through several pixels, and the distance between the peak and the trough of the difference map will give the angular expansion. When the expansion is unresolved, features move through less than a pixel. If the flux is high, we can still find the angular expansion by dividing the difference in the flux ($`\delta f`$) by the gradient ($`df/dr`$), thus recovering the angular expansion ($`\delta r`$). Dividing this angular expansion by the time elapsed gives an expansion rate. This method gives two estimates for the angular expansion, one at either end of the minor axis. This method is described in greater detail in Hajian, Terzian & Bignell (1993).

Operationally, we used IDL to make these measurements. After choosing a cut along the minor axis, we overplot the flux measurement along that axis, and the derivative ($`df/dr`$). Using these plots, we can find the difference in the flux and the derivative at the location of the peak, and therefore the angular expansion. This is done several times, for each end of the minor axis, and the average value is reported in Table 5. Errors in this method are found from the standard deviation of these measurements.

### 4.2 The Magnification Method

Following Paper I, we also determined the angular expansion of a nebula by magnifying the first-epoch image by a factor, M, and subtracting it from the second. The process was repeated for a reasonable range of M ($`1<M<1.2`$), to find the value of M which minimized the residuals. This method assumes that the entire nebula expands at the same rate, but since we do not require the entire nebula to disappear when we subtract the magnified epoch, but rather search for an M which gives the best fit to the second-epoch image, the effect of this assumption is minimized. The angular expansion rate of the radius is given by

$$\dot{\theta }=\frac{(M1)\theta }{\delta T}$$
(7)
where $`\theta `$ is the angle between the center of the nebula and the bright rim being used to determine the best fit, and $`\delta T`$ is the time elapsed between epochs of observation. This method shows graphically how the nebular expansion varies from the assumed uniform expansion. Errors are found by comparing the difference images at various magnifications. The largest step for which there is no observable change between difference images is taken as the accuracy of the measurement. For example, if $`M=1.001`$ and $`M=1.002`$ yield indistinguishable difference maps, but $`M=1.0025`$ is obviously different than either, then the error in the measurement is 0.001. (The error is determined in both directions, i.e. in the preceding example, the difference map at $`M=1.0005`$ is also distinguishable from the difference maps at $`M=1.001`$ and $`M=1.002`$. Table 4 summarizes the magnification results.

## 5 NGC 6578

### 5.1 Morphology

Figure 9 shows a color composite of the \[N II\], \[O III\] and $`V`$ images, while Figure 9a, b and c show the $`V`$ band, \[O III\] and \[N II\] images 0. All of these are images from the second epoch. In the $`V`$ image, NGC 6578 exhibits a bright central core ($`\theta =6\mathrm{}.3`$), surrounded by a fainter halo (about 11 arcseconds in diameter). The inner core appears to have two “blowout” bulbs along a single axis passing through the central star, from slightly west of north to slightly east of south. These blowout lobes are co-located with halo regions of decreased brightness. Perhaps these are places where the inner core is not confined as tightly by the halo as the rest of the core.

In the \[N II\] second-epoch image, there are a few bright knots which seem to be associated with the southeast blowout region. Two of these knots are very close together, and may be physically connected. All of these knots appear to have downstream “tails,” which point away from the central star, and away from the blowout region, similar to those in the Helix nebula, but with many fewer knots.

Along a second axis, approximately from east to west, is a much larger \[N II\] emission axis, with the lower \[N II\] emission region containing many approximately spherical blobs of \[N II\] emission. None of these knots appear to have tails of the type exhibited by the southeast knots. This may be an orientation effect, or may be real—there is no “blowout” appearance along this axis. The outer edge of this set of blobs is sharply bounded. The diffuse \[N II\] emission also appears to drop to zero at approximately the same location. While the diffuse emission extends from the central core to the outer edge of the blobs, the \[N II\] blobs appear to be entirely exterior to the bright inner core.

The inner core has a bright rim, and a blotchy inner appearance, reminiscent of a heap of soapsuds, with faint regions surrounded by bright rims. These bright rims appear to be outlines, not striations. They do not in general cross each other, but where one rim is perpendicular to another, they make a “T”. We may be viewing a tightly packed bundle of smaller hollow bubbles. When viewed in projection, they give this blotchy appearance.

The inner core of this nebula appears to be a waisted ellipsoid with axial ratio 1.3:1 and an inclination angle of $`25^{}`$ (see Table 3). However, the inner core is quite asymmetric, and this is the poorest fit of all the nebulae studied here.

### 5.2 Distance/Size Determination

The time between observations of NGC 6578 was 4.2 years. Figure 9d is the difference image from the epoch subtraction in $`V`$ band, and frames 9a, 9b and 9c are the $`V`$ band, \[O III\] and \[N II\] images respectively. These $`V`$ and \[O III\] images show a strong correspondence, implying that the kinematic \[O III\] expansion can be compared to the angular expansion measured in $`V`$. The echelle spectrum of \[O III\] gives a line of sight expansion velocity of $`19.2\pm 0.5`$ km/s.

The gradient method (as always, along the minor axis) yields an angular expansion rate of $`1.95\pm 1.0`$ mas/year, corresponding to a distance of $`1.63\pm 1.07`$ kpc. The magnification method gives a magnification factor of $`1.002\pm 0.0005`$, which corresponds to an angular expansion rate of $`1.6\pm 0.4`$ mas/yr, and a distance of $`2.00\pm 0.5`$ kpc.

## 6 NGC 6884

### 6.1 Morphology

Figure 9 shows the image data for this object. This nebula has an “s-shaped” inner core, point-symmetric around the central star, similar to NGC 7009. This inner core is surrounded by a filamentary region, perhaps a round ring inclined at an angle of approximately $`45^{}`$. This is in excellent agreement with the $`40^{}`$ inclination angle derived from the model of this nebula. Outside this ellipse is faint \[O III\] emission, which appears to be constrained to regions perpendicular to the ring.

Quite puzzling, but typical of point-symmetric bipolars, is the presence of faint \[N II\] blobs disconnected from the inner nebula (see the red ‘blobs’ in Figure 9). There are two of these, extending along an axis quite different than any of the other symmetry axes in the nebula, at position angle 135. Unfortunately, our kinematic data were taken along position angles $`111^{}`$ and $`200^{}`$, and so we have no kinematic information to further investigate whether these blobs might be FLIERs.

The inner core of this nebula is fit very well by a prolate ellipsoid with axial ratios 1.6:1, and an inclination angle of $`40^{}`$. The nebula is quite symmetrical, and therefore well-described by this elliptical model.

### 6.2 Distance/Size Determination

Four years elapsed between observations of NGC 6884. Figure 9 shows the image data. Again, the \[O III\] images are a good match to the $`V`$ band data (Figure 9), so we are confident that comparing \[O III\] velocities with $`V`$ band expansions is reasonable. The echelle data for this nebula indicate a line of sight expansion velocity of $`16.6\pm 0.4`$ km/s.

The gradient method expansion velocity of this nebula is $`1.6\pm 1.0`$ mas/year, which gives a distance of $`1.56\pm 0.98`$ kpc. Using the minor axis of the central core to determine the magnification gives an angular expansion velocity of $`1.14\pm 0.3`$ mas/yr. This corresponds to a distance of $`2.2\pm 0.8`$ kpc.

## 7 NGC 6891

### 7.1 Morphology

This nebula has an elliptical inner core (see Figure 9), with an irregular outer halo. This is in distinct contrast with the other three nebulae, where the outer halo is symmetrical, and the inner core is irregularly shaped. In detail, the inner core is fairly irregular, with brighter filaments crossing the surface. These filaments may have the same structure as the filaments in NGC 6578, but they are fainter in this nebula, so the pattern is not as clear. There is faint \[N II\] emission beyond the \[O III\] and $`V`$ band emission. This \[N II\] emission is somewhat clumpy, but it is not possible to determine whether the clumps have tails. There is no evidence of blowouts or collimated outflows.

This nebula is well described as a prolate ellipsoid with axial ratio 1.5:1, inclined about $`35^{}`$. The quite straight edges of this nebula require $`n`$ to be larger than for any other nebula in this study, about 6.0 fits the nebula well.

### 7.2 Distance/Size Determination

The echelle data yields an expansion velocity of 7.67 km/s for this nebula, which is only marginally resolved. The $`V`$ band and the \[O III\] images of this nebula do not correspond well (Figure 9). In the \[O III\] line, the halo is of comparable brightness with the inner core (the two regions are distinguished more by their morphologies than by a change in luminosity). It is probable that the \[O III\] spectroscopy is sampling the halo, rather than the inner core.

As might be expected from the low expansion velocity, the expansion was too slow to be reliably detected by either the gradient method or the magnification method. However, we can place an upper limit on the expansion, of approximately 1.6 mas/yr (the smallest expansion detected in the other objects). Combining this value with the marginally resolved expansion velocity gives a lower-bound of 1.0 kpc for the distance to NGC 6891. Correcting for instrumental smearing, the expansion velocity may be as low as 6.5 km/s, which slightly decreases this lower bound, giving 0.9 kpc.

## 8 IC 2448

### 8.1 Morphology

IC 2448 is an excellent target for this project. The nebula is elliptical and plain, with no complicating inner and outer structure (see Figure 6). The lack of field stars made image alignment difficult, but the final difference image shows a high degree of symmetry, and so we are confident that the alignment is good.

The \[N II\] and \[O III\] emissions in this nebula are coincident, and entirely diffuse, with no bright knots or filaments, consistent with a picture of an old, evolved nebula. The \[O III\] emission dominates the $`V`$ band images, as is made clear by the close correspondence in morphology between the two images (Figure 6). We can be confident that the angular expansion observed in $`V`$ band is a good match to the kinematic expansion observed in \[O III\].

The outer ring of IC 2448 appears to compress over time. That is, it gets thinner, with the inner radius increasing faster than the outer radius. This could be interpreted as an indication that a younger, faster wind is sweeping up material, except that the brightness of the ring does not appear to change. A second interpretation is that the nebula is old, the central star has ceased producing ionizing photons, and the nebula is beginning to fade. In the inner, denser regions, the nebula fades faster than in the outer, more diffuse regions. This interpretation is supported by the evolutionary age of the nebula, which McCarthy, et al (1990) report as 8400 years.

This nebula is very well fit by a prolate ellipsoid with an axial ratio of 1.3:1, and an inclination angle of $`25^{}`$.

### 8.2 Distance/Size Determination

Four years elapsed between observations of IC2448. The kinematic line of sight expansion velocity of IC 2448 is $`17.9\pm 0.3`$ km/s. Figure 9 shows the image data.

The gradient method yields an angular expansion velocity of $`2.19\pm 1.0`$ mas/yr, corresponding to a distance of $`1.41\pm 0.64`$ kpc. The magnification method gives an angular expansion velocity of $`2.25\pm 0.6`$ mas/year, or a distance of $`1.38\pm 0.4`$ kpc. At HST resolution, this nebula has no obvious variation in the proper motion with position angle.

## 9 Conclusions and Future Work

Table 6 summarizes the distance determinations of the statistical methods, as well as the expansion parallax methods under consideration here. In general, there is fair agreement between the three methods, considering the large errors on the statistical distances. IC 2448 is the glaring exception, where the difference between the distances determined in this paper and the statistical distances are more than a factor of two. Martin (1994) found a similar discrepancy via the extinction method. His observations put a firm upper limit of 1.5 kpc on the distance to this nebula, since the nebula is in the foreground of a cloud located between 1.5 and 2 kpc away. Our measurements are consistent with this firm upper limit, while the statistical distances are not.

Interestingly, all of the expansion parallax distances determined here are smaller than the statistical distances. However, we feel that it would be a mistake to generalize this to all nebulae. Not only is the current sample small, but it is also biased in favor of close nebulae. More distant nebulae require more time to observe an angular expansion, and also would have small angular diameters, and are therefore unlikely to be chosen for this study.

The magnification method and the gradient method are in good agreement in this study. Certainly they overlap to within the error bars. However, we believe there are systematics inherent in the magnification method that make the errors much more difficult to quantify. For example, when the nebula expands non-uniformly, determining the best magnification is quite subjective. For the cases in this paper, the nebulae expand fairly uniformly, and so these errors are small. When two different individuals perform the magnification method on these nebulae, their answers differ by at most $``$0.001. Operationally, this should be taken as the error in the magnification method (when propagated through the error equations, even this error is negligible compared with the error in measuring the size of the nebula). In cases where the nebula departs drastically from uniform expansion, the gradient method is more easily quantified, and therefore gives more repeatable results.

We have shown that the expansion-parallax method using the Hubble Space Telescope gives distances to PNe. For nebulae which are quite similar in \[O III\] and $`V`$ band images, and which expand uniformly, the errors are well-determined. Non-uniform expansion can cause unrecognized erroneous results, as does a discrepancy between \[O III\] and $`V`$ band images.

These results are, in general, in only fair agreement with previously determined statistical distances. This is not much of a surprise, since these statistical methods often use dubious assumptions.

Future work includes an analysis of the next twenty targets, determining their expansion parallax distances, and further constraining the systematics. Derived distant-dependent properties, such as the luminosity of the nebula and central star, the density, and the mass, will be calculated in future papers. As more targets are analyzed, we will begin to develop a better idea of the distribution and density of nebulae in the galaxy, improving constraints on chemical evolution and galactic dynamics.

We would like to thank Sean Doyle for the use of his 3-D modelling code for IDL. Support for this work was provided by NASA through grant number GO7501 and GO8390 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555.


format: text

text:
# Cryptography based on neural networks - analytical results

## Abstract

Mutual learning process between two parity feed-forward networks with discrete and continuous weights is studied analytically, and we find that the number of steps required to achieve full synchronization between the two networks in the case of discrete weights is finite. The synchronization process is shown to be non-self-averaging and the analytical solution is based on random auxiliary variables. The learning time of an attacker that is trying to imitate one of the networks is examined analytically and is found to be much longer than the synchronization time. Analytical results are found to be in agreement with simulations.

The study of neural networks was originally driven by its potential as a powerful learning and memory machine. Statistical mechanics methods have been used to analyze the network’s ability and explore its limitations Herz ; EnvB . In a recent paper KanterKinzel , the bridge between the theory of neural networks and cryptography was established. It was shown numerically that two randomly initialized neural networks with one layer of hidden units (so called Parity Machines (PMs)) learning from each other, are able to synchronize. The two parties have common inputs and they exchange information about their output. In the case of disagreement, the two PMs are trained by the Hebbian learning rule on their mutual outputs and they develop a full synchronized state of their synaptic weights. This synchronization procedure can be used to construct an ephemeral key exchange protocol for the secure transmission of secret data. An attacker, who knows the architecture of the two parties, the common inputs, and observes the mutual exchange of information, finds it difficult to imitate the moves of the parties and to reveal the common parameters after synchronization. All parties have secret informations which are not known neither to other members nor to possible attackers: Their initial weights and the current state of their hidden units, which we are noted as internal representations (IRs)

During the last decade, the analysis of learning from examples performed by feed-forward multi-layered networks was exhaustively examined using statistical mechanics methods Herz ; EnvB . An interesting network belonging to this class is the tree PM which is characterized by a superior capacity, as was found by replica calculations pm . The study of the generalization ability of such networks was based on a set of training examples generated by a static teacher network. Here we discuss a case where two or several multilayer networks are trained by their mutual outputs. This scenario has been solved only for perceptrons and only for continuous ones MetzlerKinzelKanter . Here we present an analytic solution for PMs with continuous as well as with discrete weights.

In our cryptosystem, each party in the secure channel is represented by a feed-forward network consisting of $`KN`$ random input elements $`x_{ji}=\pm 1,j=1,\mathrm{},N`$, $`K`$ binary hidden units $`\tau _i=\pm 1,i=1,\mathrm{},K`$ and one binary output unit $`\sigma =\mathrm{\Pi }_i\tau _i`$. For the simplicity of the calculations presented below we concentrate only on the case of a tree PM with $`3`$ binary hidden units feeding a binary output $`\sigma =\tau _1\tau _2\tau _3`$. The hidden units are determined via Boolean functions $`\tau _i=\text{sgn}(_jW_{ji}x_{ji})`$ through three disjointed sets of inputs $`𝐗_i=x_{1i},\mathrm{},x_{Ni}`$. The weights are either discrete or continuous, and the analytical results are derived for $`N1`$.

In this Letter we present: (a) An analytical solution of the mutual learning of two PMs whose weight vectors are updated according to the mismatch between their mutual information - their outputs. Synchronization is achieved in the case of discrete weights, $`W_{ji}=0,\pm 1,\mathrm{},\pm L`$, as well as for continuous weights confined to a sphere, $`_{j=1}^NW_{ji}^2=N`$. (b) Analysis of online adaptation of discrete weights, in which each change of a component is not infinitesimally small, demands different methods than the standard ones PercBiehl , and this is at the center of the discussion below. Surprisingly, synchronization is achieved for the discrete weights at a finite number of steps. (c) Dynamical evolution of the discrete networks cannot be characterized by the time evolution of the standard order parameters, since the overlaps between the weight-vectors are not self averaging Reents  even for large networks. The analytical solution is based on calculation of the evolution of the distribution of the order parameters as a function of the initial set of the weights. (d) The analysis is extended to include a possible attacker.

For simplicity of presentation, we first describe the analytical methods developed for the discrete case where detailed results are presented for particular examined cases. At the end of this Letter results for the continuous case are also briefly summarized.

The definition of the updating procedure between the two parties, $`A`$ and $`B`$, that are trying to synchronize their weights, is as follows. In each time step, output of the two parties is calculated for a common random input. Only weights belonging to the one (or three) hidden units which are equal to their output unit are updated, in each one of the two parties. The updating is done according to the following Hebbian learning rule,

$`W_{ji}^{A+}=W_{ji}^A+K(W_{ji}^Ax_{ji}\sigma ^B)x_{ji}\sigma ^B\theta (\sigma ^A\tau _i^A)\theta (\sigma ^A\sigma ^B),`$ (1)
$`W_{ji}^{B+}=W_{ji}^B+K(W_{ji}^Bx_{ji}\sigma ^A)x_{ji}\sigma ^A\theta (\sigma ^B\tau _i^B)\theta (\sigma ^A\sigma ^B),`$
where $`K(y)=1\delta _{L,y}`$ and $`\delta `$ represents the Kronecker function. The purpose of the operator $`K(y)`$ is to prevent the increment (decrement) of the strength of the weights on the boundary value $`L`$($`L`$).

Two important simulation results are crucial for the analytical description of the mutual dynamics. The first observation is that the synchronization time is finite KanterKinzel . The second is that different runs (set of random inputs) of the above dynamics, but with fixed initial conditions for the two parties, result in different sets of IRs. As a result of these two observations, we realized that the variance of the overlaps between the two parties is finite and does not shrink to zero even in the thermodynamic limit. This unusual scenario of on-line mutual learning is taken into consideration in the analytical equations, by the selection of random IRs following the freedom given by the current analytical overlaps. We find an iterative discrete set of equations for the mutual overlaps between the parties, whose evolution depend on some random but correlated ingredients - the current IRs, $`\{\tau _i^A\},\{\tau _i^B\}`$ (see Eq. 1).

In each time step, $`\mu `$, the mutual state of the two parties is defined by a $`(2L+1)\times (2L+1)`$ matrix, $`\overline{\overline{F^i}}(\mu )`$, where $`i`$ represents the hidden unit. The element $`f_{qr}^i`$ of the matrix stands for the fraction of components in the $`i`$th weight-vector which are equal to $`q(r)`$ in the first(second) party, where $`q,r=0,\pm 1,\mathrm{},\pm L`$. The overlap of the weights belonging to the $`i`$th hidden unit in the two parties, $`R_i^{A,B}=𝐖_i^A𝐖_i^B/N`$, as well as their norms, $`Q_i=𝐖_i𝐖_i`$/N, are given by the matrix elements

$$R_i^{A,B}=\underset{q,r}{}qrf_{qr}^i,Q_i^A=\underset{q,r}{}q^2f_{qr}^i,Q_i^B=\underset{q,r}{}r^2f_{qr}^i.$$

These overlaps and norms fixed the probabilities of deriving the same IR via the normalized overlap, $`\rho _i^{A,B}=R_i^{A,B}/\sqrt{Q_i^AQ_i^B}`$. More precisely, the probability of having different results in the $`i`$th hidden unit of the two parties is given by the well known generalization error for the perceptron $`ϵ_p^i=\mathrm{cos}^1\rho _i/\pi `$ Herz ; EnvB .

Each of the PM consists of a tree architecture and for random inputs each of the $`8`$ IRs appears with equal probability. The joint probability distribution of the $`64`$ different pairs of IRs in both parties is correlated, and can be explicitly expressed using $`\{ϵ_p^i\}`$.

The development of the elements of the matrix $`\overline{\overline{F^i}}(\mu )`$ are calculated directly from Eq. 1, where one has to average over the inputs $`x_{ij}`$. We use auxiliary random variables in order to choose one of the possible IRs following their probabilities given by $`\{ϵ_p^i\}`$. In each step we choose two sets of random numbers which are taken from a flat distribution between $`0`$ and $`1`$: Set I: In the event that the number is smaller than $`ϵ_i`$ we deduce that the two hidden units disagree, otherwise we assume an agreement. Set II: All eight IRs are equally probable in the first party, since the architecture consists of a tree PM. We choose one among the eight using the second set of auxiliary variables $`p_r`$, and the corresponding IR for the second network according to the first set.

To exemplify derivation of the iterative equations for $`\{f_{qr}^i\}`$, let us concentrate on the case where the result of the first random set is that all three hidden units are in disagreement. In two possibilities out of the eight IRs all three hidden units are updated, whereas in the other six possibilities only one is updated (we then have to choose at random one among the three). After taking into account all possible internal scenarios, and accordingly the updates, one can show that the iterative equations for $`\{f_{qr}^i\}`$ away from the boundary, $`q,r\pm L`$, are given by

$`f_{q,r}^{i+}=\theta ({\displaystyle \frac{1}{4}}p_r)({\displaystyle \frac{1}{2}}f_{q+1,r1}^i+{\displaystyle \frac{1}{2}}f_{q1,r+1}^i)+`$
$`\theta ({\displaystyle \frac{i+1}{4}}p_r)\theta (p_r{\displaystyle \frac{i}{4}})({\displaystyle \frac{1}{2}}f_{q+1,r1}^i+{\displaystyle \frac{1}{2}}f_{q1,r+1}^i).`$
On the boundary, similar equations can be derived as well as for other internal scenarios. Taking into account all possible scenarios and the inversion symmetry of our PMs, one has to solve iteratively only $`4`$ classes of equations in a manner similar to the abovementioned long\_paper . Note that the time evolution of the $`f_{qr}`$ and the overlaps depends on time dependent random variables.

Different runs for updating of the equations result in different trajectories of the order parameters. In the inset of Fig. 1, we present the average overlap $`\overline{\rho }=_{i=1}^3\rho _i/3`$, and its standard deviation, obtained from $`500`$ different runs of the analytical equations with $`L=1`$. Results of the averaged overlap (with the same standard deviation) obtained in $`500`$ runs of simulations with $`N=10^4`$ are denoted by circles.

An important quantity is the number of steps required to achieve full synchronization, $`t_{synch}`$, since it can be used by the parties to encrypt/decrypt the information using the known output bit. In simulations the synchronization time is well defined - the first step in which all weight vectors of the parties are in an anti-parallel state. In contrast, in the analytical solution the average overlap of the hidden units tends to zero exponentially with the number of steps. In order to compare analytical results to simulations we need to find a criterion which determines synchronization. We chose the criterion $`\overline{\rho }c_l=1+0.1/(NL)`$ to define full synchronization, since $`c_l`$ is much greater than the maximal possible overlap just before synchronization.

The exponential decay of the overlaps with the number of steps and the claim that synchronization is achieved at a finite number of steps even for $`N1`$ has to be clarified. Our synchronization process is mainly characterized by two regimes: The first $`t_a`$ steps which are characterized by different IRs (in some of the steps) for the two parties. Note that $`t_a`$ is fluctuating from sample to sample. The second is the asymptotic regime, last $`t_b`$ steps, where the IRs of the parties are always the same, and the weights are converging to an anti-parallel state similarly to three perceptrons, $`t_b\mathrm{log}(N)`$ long\_paper . Roughly speaking, the two regimes are characterized by $`ϵ_i^P>1/t_a`$ and $`ϵ_i^P<1/t_b`$, respectively. Our analytical results as well as simulations indicate that $`t_a`$ is independent of $`N`$. Hence as long as $`t_a>t_b`$, the $`\mathrm{log}(N)`$ dependent is invisible. For $`L=3`$, for instance, $`t_{synch}400`$, $`t_a300`$, and $`t_b`$ is expected to be equal to $`t_a`$ only for $`Ne^{200}`$.

In Fig. 1, we present the histogram of the number of steps required to achieve $`t_{synch}`$, $`P(t_{synch})`$, in simulations with $`N=10^4`$ and $`L=1`$ and the initial weights were chosen such that $`\rho _i^{A,B}=0`$. This distribution is in a fairly good agreement with the results obtained by the runs of the iterative equations for $`f_{qr}`$.

Let us now examine a possible attack of a third player, an attacker $`o`$, that tries to imitate one of the parties (let us say $`A`$). We assume that the attacker uses the same algorithm as one of the partners. The attacker updates its own weight-vectors only when an updating step is taken by the parties. The natural move of an attacker in such an event is to follow the rule of the parties

$$W_{ji}^{o+}=W_{ji}^o+K(W_{ji}^ox_{ji}\sigma ^B)x_{ji}\sigma ^B\theta (\sigma ^A\tau _i^o)\theta (\sigma ^A\sigma ^B),$$

indicating that only weight-vectors belonging to the hidden units which are in agreement with the output of party $`A`$ are updated, ( more advanced attacks will be discussed elsewhere long\_paper ). The evolution of the overlap of an attacker depends on the evolution of $`6`$ matrices; three matrices describing the overlaps between the parties and similarly, three matrices describing the overlaps between the attacker and the first party. Note that the dynamics of the attacker depends on moves of the parties which depend on their overlaps. Hence, the time evolution of six matrices gives the full description of the overlaps between the attacker and the first party and between the parties themselves. The mutual dynamics of the three networks, two parties and the attacker, depends on the joint probability distribution of $`8\times 8\times 8`$ IRs, and upon the corresponding updates of the six matrices. The full description of the discrete time evolution of the matrices and the overlaps will be given elsewhere long\_paper .

The analytical solution of the dynamics of the attacker indicates that a full learning is achieved in a finite number of steps, $`t_{learn}`$, where a full learning is defined such that $`\overline{\rho }^{A,o}>c_l`$. In Table 1 $`t_{learn}`$ and $`t_{synch}`$ are compared for various $`L`$.

For $`L=1`$ the average learning time is about twice the synchronization time, and one may reach the wrong conclusion that the synchronization process always terminates before the learning process. In Fig. 1 we present the histogram of the synchronization and the learning processes, and a fairly good fit between analytical and simulation results is apparent. The two distributions, $`P(t_{synch}),P(t_{learn})`$ have a finite overlap, indicating that in a finite fraction of the runs the learning process terminates before the achievement of synchronization (which was indeed observed in a finite fraction of the runs of the simulations). Hence the construction with $`L=1`$ is not a good candidate to build a secure channel.

For $`L3`$ the ratio $`r=t_{learn}/t_{synch}`$ averaged over the runs was found to be $`r1`$ (see Table 1). For $`L=3`$, we did not observe, in simulations over $`10^5`$ runs, a case where $`t_{learn}`$ was faster than $`t_{synch}`$. In Fig. 2 we present the histogram of the probabilities of the ratio, $`r`$, as was found by averaging over different runs of the analytical equations. The minimal value of the ratio was $`r6`$ where the largest ratio was $`r680`$. We found that the largest synchronization times are smaller than $`1000`$ whereas the typical learning time is $`4.510^4`$.

Synchronization in the case of PMs with continuous weights is achievable only with the following modifications. (a) Normalization of the weight vectors belonging to each one of the hidden units after every updating step. The natural normalization we use is the spherical normalization, $`_{j=1}^NW_{ji}^2=N`$. (b) The change in the strength of each weight (before normalization) is $`\eta /N`$, where $`\eta `$ is a constant of order one. The synchronization time is proportional to the size of the input, $`N`$, and therefore the analytical description of the system is given by a coupled differential equations. Some limited results and brief description of the method are presented below. More detailed results will be given elsewhere long\_paper .

Updating of weights of the first party for the spherical case is given by

$$𝐖_i^{A+}=\frac{𝐖_i^A+\frac{\eta }{N}𝐗_i\theta (\sigma ^A\sigma ^B)\theta (\sigma ^A\tau _i^A)\sigma ^B}{|𝐖_i^A+\frac{\eta }{N}𝐗_i\theta (\sigma ^A\sigma ^B)\theta (\sigma ^A\tau _i^A)\sigma ^B|}$$

and similarly the updating rules for the second party and the attacker. The analytical calculation can be simplified in the continuous case by the probability that there is a mismatch between the two PMs given that there is a mismatch between two hidden units, $`P_1^iP(\sigma ^A\sigma ^B|\tau _i^A\tau _i^B)=ϵ_p^jϵ_p^k+(1ϵ_p^j)(1ϵ_p^k)`$ and similarly $`P_2^i=P(\sigma ^A\sigma ^B|\tau _i^A=\tau _i^B)=1P_1^i`$. One can map the mutual process onto that of perceptrons, where the updating of the first party, for instance, is given by

$$𝐖_i^{A+}=(𝐖_i^A+\frac{\eta }{N}𝐗_i\tau ^B\mathrm{\Delta }_i^A)/|𝐖_i^A+\frac{\eta }{N}𝐗_i\tau ^B\mathrm{\Delta }_i^A|$$

and similarly for the second party, where $`\mathrm{\Delta }_i^A=\theta (\tau _i^A\tau _i^B)\theta (\frac{P_1^i}{2}p_a)+\theta (\tau _i^A\tau _i^B)\theta (P_2^ip_b)\theta (\frac{1}{2}p_c)`$ and we use auxiliary variables $`p_a`$, $`p_b`$, $`p_c`$ to specify each run.

The next step consists of the averages over the following two quantities. (a) Averaging over the joint probability distributions of the local fields of the two parties. (b) Average over the auxiliary variables, which is unique to the case of mutual learning. The normalized overlap, $`\rho `$, between weight vectors belonging to each pair of hidden units is found to obey the equation, $`d\rho /d\alpha =\eta [C^2+(1C)^2]((1\rho )/\sqrt{2\pi }\eta C/2)(1+\rho )2\eta (1\rho ^2)C(1C)/\sqrt{2\pi }\eta ^2\rho C(1C)^2`$, where $`C=\mathrm{cos}^1\rho /\pi `$. For $`\eta <\eta _c2.68`$ the points $`\rho =\pm 1`$ are repulsive fixed points of the above equation, where for $`\eta >\eta _c`$ a phase transition occurs to a state of full synchronization.

The equation of motion of the overlap of an attacker with the first party after synchronization, i.e., $`\rho ^{A,B}=1`$, $`\rho ^{A,o}=\rho ^{B,o}`$, is given by $`d\rho ^{A,o}/d\alpha =\eta ^2(1\mathrm{cos}^1\rho ^{A,o}/\pi \rho ^{A,o})/2`$. The fixed point of this equation is $`\rho ^{A,o}=\rho ^{B,o}0.79`$ and is independent of $`\eta `$, indicating that perfect learning is not achievable. Analytical results derived from the last two equations in the case of $`\eta =3`$ are presented in the inset of Fig. 2 and are in good agreement with simulation with $`N=5000`$ and $`20`$ runs.

I.K. acknowledges the support from the Israel Academy of Science.


text:
# Keywords

## Keywords

Solar Flare, GPS, ionosphere

## 1 Introduction

The enhancement of X-ray and ultraviolet (UV) emission that is observed during chromospheric flares on the Sun immediately causes an increase in electron density in the ionosphere. These density variations are different for different altitudes and are called Sudden Ionospheric Disturbances, SID (Davies,1990; Donnelly, 1969). SIDs are generally recorded as the short wave fadeout, SWF (Stonehocker, 1970), sudden phase anomaly, SPA (Ohshio, 1971), sudden frequency deviation, SFD (Donnelly, 1971; Liu et al., 1996), sudden cosmic noise absorption, SCNA (Deshpande and Mitra, 1972), sudden enhancement/decrease of atmospherics, SES (Sao et al., 1970). Much research is devoted to SID studies, among them a number of thorough reviews (Mitra, 1974; Davies, 1990).

Highly informative technique is the method of Incoherent Scatter (IS). The Millstone Hill IS facility recorded a powerful flare on August 7, 1972 (Mendillo and Evans, 1974a). The measurements were made in the height range from 125 to 1200 km. The increase of local electron density $`N_e`$ made up 100% at 125 km altitude and 60% at 200 km.

Using the IS method Thome and Wagner (1971) obtained important evidence of the height distribution of the increase in $`N_e`$ at the time of the May 21 and 23, 1967 flares. A significant increase of $`N_e`$ was recorded in the E-region, up to 200%, which gradually decreases in the F-region with the increasing height, down to 10-30%, and remains distinguishable up to 300 km. The earliest increase of $`N_e`$ begins in the E-region, and at higher altitudes it is observed with a delay which is particularly pronounced at F-region heights.

A sudden increase in total electron content (TEC) can be measured using continuously operating radio beacons installed on geostationary satellites. On August 7, 1972, Mendillo et al. (1974b) were the first to make an attempt to carry out global observations of the solar flare using 17 stations in North America, Europe, and Africa. The observations covered an area, the boundaries of which were separated by $`70^{}`$ in latitude and by 10 hours in local time. For different stations, the absolute value of the TEC increase $`\mathrm{\Delta }I`$ varies from $`1.810^{16}`$ to $`8.610^{16}\text{el}\text{m}^2`$, which corresponds to 15-30% of the TEC. Investigations revealed a latitudinal dependence of the TEC increase value. At low latitudes, it was higher compared with high latitudes. Besides, the authors point out the absence of a connection between the TEC increase value and the solar zenith angle.

The advent and evolution of a Global Positioning System (GPS) and also the creation on its basis of widely branched networks of GPS stations (at least 900 sites at the August of 2001, the data from which are placed on the Internet) opened up a new era in remote ionospheric sensing. High-precision measurements of the TEC along the line-of-sight (LOS) between the receiver on the ground and transmitters on the GPS system satellites covering the reception zone are made using two-frequency multichannel receivers of the GPS system at almost any point of the globe and at any time simultaneously at two coherently coupled frequencies $`f_1=1575.42`$ MHz and $`f_2=1227.60`$ MHz.

The sensitivity of phase measurements in the GPS system is sufficient for detecting irregularities with an amplitude of up to $`10^310^4`$ of the diurnal TEC variation. This makes it possible to formulate the problem of detecting ionospheric disturbances from different sources of artificial and natural origins. The TEC unit (TECU) which is equal to $`10^{16}`$ $`\text{el}\text{m}^2`$ and is commonly accepted in the literature, will be used throughout the text.

Afraimovich (2000a); Afraimovich et al. (2000b, 2001a,b) developed a novel technology of a global detection of ionospheric effects from solar flares and presented data from first GPS measurements of global response of the ionosphere to powerful impulsive flares of July 29, 1999, and December 28, 1999. Authors found that fluctuations of TEC are coherent for all stations on the dayside of the Earth. The time profile of TEC responses is similar to the time behavior of hard X-ray emission variations during flares in the energy range 25-35 keV if the relaxation time of electron density disturbances in the ionosphere of order 50-100 s is introduced. No such effect on the nightside of the Earth has been detected yet.

Afraimovich et al. (2001c) and Leonovich et al. (2001) suggested a technique for estimating the ionospheric response to weak solar flares (of X-ray class C). They obtained a dependence of the ionospheric TEC increase amplitude (during the solar flare) on the flare location on the Sun (on the central meridian distance, CMD). For flares lying nearer to the disk center (CMD $`<40^{}`$), an empirical dependence of the ionospheric TEC increase amplitude on the peak power of solar flares in the X-ray range was obtained (using data from the geostationary GOES-10 satellite).

This paper is a logical continuation of the series of our publications (Afraimovich, 2000a; Afraimovich et al., 2000b, 2001a,b,c, Leonovich et al., 2001) devoted to the study of ionospheric effects of solar flares, based on data from the international GPS network.

A limitation of the GPS method is that its results have an integral character, as a consequence of which from measurements at a single site it is impossible to determine which ionospheric region makes the main contribution to the TEC variation. The objective of this study is to develop a method which would help overcome (at least partially) this problem.

## 2 Method of determining the shadow altitude $`h_0`$ over the ground

The method uses the effect of partial ”shadowing” of the atmosphere by the terrestrial globe. Direct beams of solar ionizing radiation from the flare do not penetrate the region of the Earths total shadow. GPS stations located near the shadow boundary on the ground in the nightside hemisphere are used to investigate the solar flare influence on the ionosphere. The LOS for these stations pass partially through the atmosphere lying in the total shadow region, and partially through the illuminated atmosphere. The altitude over the ground at which the LOS intersects the boundary of the total shadow cone, will be referred to as the shadow altitude $`h_0`$.

Fig. 1 schematically represents the formation of the cone of the Earth s total shadow (not to scale) in the geocentric solar-ecliptic coordinate system (GSE): the axis Z is directed to a north perpendicularly planes of an ecliptic, the axis X - on the Sun, the axis Y is directed perpendiculary to these axes. For definition of the shadow altitude $`h_0`$ it is necessary to know coordinates of a cross point C of the LOS and the shadow boundary.

The primary data are the geographical coordinates of station GPS on the Earth (Fig.1; a point P): an elevation angle and azimuth of LOS on a satellite GPS, toward the north clockwise, for the time (UT) corresponding to the phase of solar flare maximum in the X-ray range. These coordinates are converted to the Cartesian coordinate system where the Cartesian coordinates of the GPS station on the ground and the coordinates of the subionospheric point (at 300 km altitude) are calculated. Next, we use the geocentric solar-ecliptic coordinate system following the technique reported by Sergeev and Tsyganenko (1980). To determine the coordinates of the point C we solve a system of equations: the equation of cone (of total shadow), and the equation of a straight line (LOS) specified parametrically. After that, from the resulting point C we drop a perpendicular to the ground and calculate its length (Fig. 1, line $`h_0`$). The value of $`h_0`$, thus obtained, is just the shadow altitude.

## 3  Method of determining the TEC increase in the ionosphere using data from the global GPS net work

This paper exemplifies an analysis of the ionospheric effect of a powerful solar flare of class X5.7/3B recorded on July 14, 2000 (10:24 UT, N22W07) under quiet geomagnetic conditions (Dst = -10 nT). The time profile of soft X-ray emission in the range 1-8 Å (GOES-10 data) at the time of the flare is presented in Fig. 2a.

To determine the TEC increase in the ionosphere we used the data from the international GPS network. The GPS technology provides a means of estimating the TEC variations $`I_0(t)`$ on the basis of TEC phase measurements made with each of the spatially separated two-frequency GPS receivers using the formula (Calais and Minster, 1996):

$$I_0(t)=\frac{1}{40.308}\frac{f_1^2f_2^2}{f_1^2f_2^2}[(L_1\lambda _1L_2\lambda _2)+const+nL],$$
(1)
where $`L_1\lambda _1`$ and $`L_2\lambda _2`$ are the increments of the radio signal phase path caused by the phase delay in the ionosphere (m); $`L_1`$ and $`L_2`$ stand for the number of complete phase rotations, and $`\lambda _1`$ and $`\lambda _2`$ are the wavelengths (m) for the frequencies $`f_1`$ and $`f_2`$, respectively; $`const`$ is some unknown initial phase path (m); and $`nL`$ is the error in determining the phase path (m).

Input data used in the analysis include series of the oblique value of TEC $`I_0(t)`$, as well as corresponding series of elevations $`\theta `$ and azimuths of LOS to the satellite. These parameters are calculated using our developed CONVTEC program to convert standard (for the GPS system) RINEX-files received via the Internet. Input series of TEC $`I_0(t)`$ are converted to the vertical value following a well-known technique (Klobuchar, 1986).

$$I(t)=I_0cos\left[arcsin\left(\frac{R_E}{R_E+h_{max}}cos\theta \right)\right]$$
(2)
where $`R_E`$ is Earth’s radius; and $`h_{max}`$ is the height of the ionospheric $`F2`$-layer maximum.

Variations of the regular ionosphere, and also trends introduced by the motion of the satellite are eliminated using the procedure of removing the trend defined as a polynomial of the third order on a given temporary interval.

Figs. 2b and 2d present the typical time dependencies of the vertical TEC $`I(t)`$ for sites GPS WDLM (PRN02, shadow altitude $`h_0=17`$ km) and LEEP (PRN07, shadow altitude $`h_0=586`$ km). Time dependencies of the TEC $`\mathrm{\Delta }I(t)`$ response, with the trend removed for these series, are presented in Figs. 2c and 2e, respectively.

## 4  Results and discussion

The TEC response to the solar flare was analyzed for 45 GPS stations. Detailed information about the GPS stations and analysis results is summarized in Table 1: names of GPS receiving stations (Site), number of the GPS satellite from which the signal is received (PRN), shadow altitude above the ground ($`h_0`$), absolute increase of TEC $`\mathrm{\Delta }I`$, relative increase of TEC ($`\mathrm{\Delta }I(t)/\mathrm{\Delta }I_{00}(t)`$, and geographical coordinates of GPS stations (latitude, longitude). The increase of TEC $`\mathrm{\Delta }I_{00}`$ corresponds to the amplitude of the TEC increase measured at the station lying at the shadow boundary on the ground ($`h_0=0`$).

Fig. 3b illustrates examples of time dependencies of the TEC $`\mathrm{\Delta }I_0`$ response for LOS to the satellite which intersect the boundary of the shadow cone at different heights $`h_0`$ during the solar flare of July 14, 2000. Fig. 3b (left) presents the values of these altitudes, and (right) the names of corresponding stations. For a better visualization, the dependencies are drawn by lines of a different thickness. It should be noted that the response remains pronounced when the shadow altitude exceeds significantly the electron density peak height in the ionosphere. For station GUAM (PRN26, height of the shadow boundary $`h_0=885`$ km), the response amplitude exceeds the background oscillation amplitude by more than a factor of 2.

It is evident from Fig. 3 that the wave phase (time of the response maximum) is different at different altitudes $`h_0`$. On the one hand, this phenomenon can be caused by the interference of the response with background fluctuations; on the other, this can be due to the fact that at different heights different wavelengths of ionizing radiation are observed, which, in turn, can have independent time characteristics.

The dependence of the absolute TEC increase on the altitude $`h_0`$ for all the cases under consideration is plotted in Fig. 4a. The dependence of the relative TEC increase $`\mathrm{\Delta }I(t)/\mathrm{\Delta }I_{00}`$ on the altitude $`h_0`$ during the solar flare is shown as a percentage in Fig. 4b. The TEC increase $`\mathrm{\Delta }I_{00}(t)`$ corresponds to the amplitude of the TEC increase measured at the station lying at the shadow boundary on the ground ($`h_0=0`$).

Fig. 4b suggests that about 20% of the TEC increase correspond to the ionospheric region lying below 100 km, about 5% refer to the ionospheric E-region (100-140 km), about 30% correspond to the ionospheric F1-region (140-200 km), and about 30% to regions lying above 300 km. We found that a rather significant contribution to the TEC increase is made by ionospheric regions lying above 300 km.

The estimate obtained is consistent with the findings reported by Mendillo and Evans (1974a); Mendillo et al. (1974b). The authors of the cited references, based on investigating the electron density profile in the height range from 125 km to 1200 km using the IS method, concluded that about 40% of the TEC increase during the powerful flare on August 7, 1972, correspond to ionospheric regions lying above 300 km. However, Thome and Wagner (1971), who used the IS method to investigate the ionospheric effects from two others powerful solar flares, pointed out that an increase in electron density associated with the solar flare was observable to 300 km altitude only. This difference can be explained by the fact that each particular solar flare is a unique event which is characterized by its own spectrum and dynamics in the flare process.

Acknowledgements. Authors are grateful to V.G. Mikhalkovsky for his assistance in preparing the English version of the manuscript. This work was done with support from both the Russian foundation for Basic Research (grant 00-05-72026) and RFBR grant of leading scientific schools of the Russian Federation No. 00-15-98509.

References

Afraimovich, E. L., GPS global detection of the ionospheric response to solar flares, Radio Sci., 35, 1417–1424, 2000a.

Afraimovich, E. L., Kosogorov, E. A., and L. A. Leonovich, The use of the international GPS network as the global detector (GLOBDET) simultaneously observing sudden ionospheric disturbances, Earth Planet. Space, 52, 1077–1082, 2000b.

Afraimovich, E.L., Altyntsev, A.T., Kosogorov, E.A., Larina, N.S., and L. A. Leonovich, Detecting of the Ionospheric effects of the solar flares as deduced from global GPS network data, Geomagnetism and Aeronomy, 41, 208–214, 2001a.

Afraimovich, E.L., Altyntsev, A.T., Kosogorov, E.A., Larina, N.S., and L. A. Leonovich, Ionospheric effects of the solar flares of September 23, 1998 and July 29, 1999 as deduced from global GPS network data, J. Atm. Solar-Terr. Phys., 63, 1841-1849, 2001b.

Afraimovich, E.L., Altyntsev, A.T., Grechnev, V.V., and L. A. Leonovich, Ionospheric effects of the solar flares as deduced from global GPS network data, Adv. Space Res., 27, 1333-1338, 2001c.

Calais, E., and J. B. Minster, GPS detection of ionospheric perturbations following a Space Shuttle ascent, Geophys. Res. Lett., 23, 1897–1900, 1996.

Davies K., Ionospheric radio, Peter Peregrinus, London, 1990.

Deshpande, S. D. and A. P. Mitra,  Ionospheric effects of solar flares, IV, electron density profiles deduced from measurements of SCNA’s and VLF phase and amplitude, J. Atmos. Terr. Phys., 34, 255–259, 1972.

Donnelly, R. F., Contribution of X-ray and EUV bursts of solar flares to Sudden frequency deviations, J. Geophys. Res., 74, 1873 – 1877, 1969.

Donnelly, R. F., Extreme ultraviolet flashes of solar flares observed via sudden frequency deviations: experimental results, Solar Phys., 20, 188–203, 1971.

Klobuchar, J. A., Ionospheric time-delay algorithm for single-frequency GPS users, IEEE Transactions on Aerospace and Electronics System, AES , 23(3), 325–331, 1986.

Leonovich, L.A., Altynsev, A.T., Afraimovich, E. L., and V. V. Grechnev, Ionospheric effects of the solar flares as deduced from global GPS network data. LANL e-print archive, 2001, http://arXiv.org/abs/physics/011006.

Liu, J. Y., Chiu, C.S., and C. H. Lin, The solar flare radiation responsible for sudden frequency deviation and geomagnetic fluctuation, J. Geophys. Res., 101, 10855–10862, 1996.

Mendillo, M., and J.V. Evans, Incoherent scatter observations of the ionospheric response to a large solar flare, Radio Sci., 9, 197-203, 1974a.

Mendillo, M., Klobuchar, J. A., Fritz, R. B., da Rosa, A.V., Kersley, L., Yeh, K. C., Flaherty, B. J., Rangaswamy, S., Schmid, P. E., Evans, J. V., Schodel, J. P., Matsoukas, D. A., Koster, J. R., Webster, A. R., and P. Chin, Behavior of the Ionospheric F Region During the Great Solar Flare of August 7, 1972, J. Geophys. Res., 79, 665–672, 1974b.

Mitra, A. P., Ionospheric effects of solar flares, D.Reidel, Norwell, Mass., 1974.

Ohshio, M., Negative sudden phase anomaly, Nature, 229, 239–244, 1971.

Sao, K., Yamashita, M., Tanahashi, S., Jindoh, H., and K. Ohta, Sudden enhancements (SEA) and decreases (DSA) of atmospherics, J. Atmos. Terr. Phys., 32, 1567–1573, 1970.

Sergeev, V. A., and N. A. Tsyganenko,  The earth’s magnetosphere. Results of researches on the international geophysical projects. ’Nauka’, Moscow, (in Russian), 1980.

Stonehocker, G. H., Advanced telecommunication forecasting technique in AGY, 5th., Ionospheric forecasting, AGARD Conf. Proc., 29, 27-31, 1970.

Thome, G.D, and L. S. Wagner,  Electron density enhancements in the E and F regions of the ionosphere during solar flares, J.Geophys. Res., 76, 6883–6895, 1971.


text:
# References

| hep-ph/0202152 |
| --- |
| 16 February 2002 |

$`\nu _e\nu _s`$ oscillations with large neutrino mass in NuTeV?
Carlo Giunti

INFN, Sezione di Torino,
and
Dipartimento di Fisica Teorica, Università di Torino,
Via P. Giuria 1, I–10125 Torino, Italy
Marco Laveder

Dipartimento di Fisica “G. Galilei”, Università di Padova,
and
INFN, Sezione di Padova,
Via F. Marzolo 8, I–35131 Padova, Italy

## Abstract

We propose an explanation of the NuTeV anomaly in terms of oscillations of electron neutrinos into sterile neutrinos. We derive an average transition probability $`P_{\nu _e\nu _s}0.21\pm 0.07`$, which is compatible with other neutrino data if the mass-squared difference that drives the oscillations is $`\mathrm{\Delta }m^210100\mathrm{eV}^2`$.

The NuTeV collaboration measured recently  a value of the electroweak parameter $`\mathrm{sin}^2\theta _W`$ higher than the standard model prediction obtained from a fit of other electroweak data by

$$\mathrm{sin}^2\theta _W^{\mathrm{NuTeV}}\mathrm{sin}^2\theta _W=\left(5.0\pm 1.6\right)\times 10^3.$$
(1)
In this paper we propose an explanation of this discrepancy in terms of $`\nu _e\nu _s`$ oscillations111 Other possible explanations of the NuTeV anomaly have been discussed in Ref. . , where $`\nu _s`$ is a sterile neutrino that does not participate to weak interactions (see ). Our analysis is approximate and based on the limited information contained in the paper of the NuTeV collaboration . A precise analysis of NuTeV data in terms of neutrino oscillations can be done only by the NuTeV collaboration.

The NuTeV collaboration measured the ratio of neutral current to charged current cross sections for neutrino scattering with isoscalar targets of $`u`$ and $`d`$ quarks

$$R^\nu =\frac{\sigma (\nu _\mu +N\nu +X)}{\sigma (\nu _\mu +N\mu ^{}+X)},$$
(2)
and the corresponding ratio $`R^{\overline{\nu }}`$ for antineutrinos. The Paschos-Wolfenstein relation  allows to derive the value of $`\mathrm{sin}^2\theta _W`$ from $`R^\nu `$ and $`R^{\overline{\nu }}`$:

$$\mathrm{sin}^2\theta _W=\frac{1}{2}\frac{R^\nu rR^{\overline{\nu }}}{1r},$$
(3)
where

$$r\frac{\sigma (\overline{\nu }_\mu +N\mu ^++X)}{\sigma (\nu _\mu +N\mu ^{}+X)}\frac{1}{2}.$$
(4)
The NuTeV collaboration distinguishes CC and NC $`\nu _\mu `$ interactions on the basis of the length of the track in the detector: long tracks correspond to CC $`\nu _\mu `$ interactions and short tracks correspond to NC $`\nu _\mu `$ interactions, plus a contribution from CC and NC $`\nu _e`$ interactions. The number $`N_{\nu _e}^{\mathrm{MC}}`$ of $`\nu _e`$ interactions calculated through a Monte Carlo is subtracted to the number $`N_\mathrm{S}^\nu `$ of short tracks in order to obtain the number of NC $`\nu _\mu `$ interactions used to evaluate $`R^\nu `$:

$$R_{\mathrm{NuTeV}}^\nu =\frac{N_\mathrm{S}^\nu N_{\nu _e}^{\mathrm{MC}}}{N_\mathrm{L}^\nu },$$
(5)
where $`N_\mathrm{L}^\nu =N_{\nu _\mu }^{\mathrm{CC}}`$ is the number of long tracks. If $`\nu _e\nu _s`$ oscillations deplete the number of $`\nu _e`$ interactions with respect to the calculated one, the number $`N_{\nu _e}^{\mathrm{MC}}`$ subtracted in the numerator of Eq. (5) is too high, leading to a smaller value of $`R_{\mathrm{NuTeV}}^\nu `$ with respect to the true one and, from the Paschos-Wolfenstein relation (3), a value of $`\mathrm{sin}^2\theta _W^{\mathrm{NuTeV}}`$ larger than the true $`\mathrm{sin}^2\theta _W`$.

In the presence of $`\nu _e\nu _s`$ transitions, the true number of $`\nu _e`$ interactions is

$$N_{\nu _e}=N_{\nu _e}^{\mathrm{MC}}P_{\nu _e\nu _e},$$
(6)
where

$$P_{\nu _e\nu _e}=1P_{\nu _e\nu _s}.$$
(7)
is the electron neutrino survival probability averaged over the energy spectrum of the NuTeV beam. Writing

$$N_\mathrm{S}^\nu =N_{\nu _\mu }^{\mathrm{NC}}+N_{\nu _e},$$
(8)
from Eq. (5) the ratio $`R_{\mathrm{NuTeV}}^\nu `$ reported by the NuTeV collaboration is given by

$$R_{\mathrm{NuTeV}}^\nu =\frac{N_{\nu _\mu }^{\mathrm{NC}}}{N_{\nu _\mu }^{\mathrm{CC}}}\frac{N_{\nu _e}^{\mathrm{MC}}}{N_{\nu _\mu }^{\mathrm{CC}}}\left(1P_{\nu _e\nu _e}\right)=R^\nu \frac{N_{\nu _e}^{\mathrm{MC}}}{N_{\nu _\mu }^{\mathrm{CC}}}\left(1P_{\nu _e\nu _e}\right),$$
(9)
where $`R^\nu `$ is the true value of the ratio (2). A similar relation holds for the antineutrino ratio:

$$R_{\mathrm{NuTeV}}^{\overline{\nu }}=R^{\overline{\nu }}\frac{N_{\overline{\nu }_e}^{\mathrm{MC}}}{N_{\overline{\nu }_\mu }^{\mathrm{CC}}}\left(1P_{\nu _e\nu _e}\right),$$
(10)
with the same survival probability, because of CPT invariance (we assume that $`\nu _e`$ and $`\overline{\nu }_e`$ have approximately the same energy spectrum).

The NuTeV collaboration reported that

$$\frac{N_{\nu _e}^{\mathrm{MC}}}{N_\mathrm{L}^\nu +N_\mathrm{S}^\nu }\frac{N_{\overline{\nu }_e}^{\mathrm{MC}}}{N_\mathrm{L}^{\overline{\nu }}+N_\mathrm{S}^{\overline{\nu }}}1.7\times 10^2.$$
(11)
In a first approximation we have

$$\frac{N_{\nu _e}^{\mathrm{MC}}}{N_{\nu _\mu }^{\mathrm{CC}}}=\frac{N_{\nu _e}^{\mathrm{MC}}}{N_{\nu _\mu }^{\mathrm{NC}}+N_{\nu _\mu }^{\mathrm{CC}}}\frac{N_{\nu _\mu }^{\mathrm{NC}}+N_{\nu _\mu }^{\mathrm{CC}}}{N_{\nu _\mu }^{\mathrm{CC}}}\frac{N_{\nu _e}^{\mathrm{MC}}}{N_\mathrm{L}^\nu +N_\mathrm{S}^\nu }\left(1+R_{\mathrm{NuTeV}}^\nu \right),$$
(12)
and a similar relation for $`N_{\overline{\nu }_e}^{\mathrm{MC}}/N_{\overline{\nu }_\mu }^{\mathrm{CC}}`$. From Eq. (11) and the values $`R_{\mathrm{NuTeV}}^\nu R_{\mathrm{NuTeV}}^{\overline{\nu }}0.4`$ reported by the NuTeV collaboration , we obtain

$$\frac{N_{\nu _e}^{\mathrm{MC}}}{N_{\nu _\mu }^{\mathrm{CC}}}\frac{N_{\overline{\nu }_e}^{\mathrm{MC}}}{N_{\overline{\nu }_\mu }^{\mathrm{CC}}}2.4\times 10^2.$$
(13)
From Eqs. (7), (9), (10), (13) and the Paschos-Wolfenstein relation (3), the value $`\mathrm{sin}^2\theta _W^{\mathrm{NuTeV}}`$ reported by the NuTeV collaboration is given by

$$\mathrm{sin}^2\theta _W^{\mathrm{NuTeV}}=\frac{1}{2}\frac{R_{\mathrm{NuTeV}}^\nu rR_{\mathrm{NuTeV}}^{\overline{\nu }}}{1r}=\mathrm{sin}^2\theta _W+\frac{N_{\nu _e}^{\mathrm{MC}}}{N_{\nu _\mu }^{\mathrm{CC}}}P_{\nu _e\nu _s}.$$
(14)
Using the numerical values in Eqs. (1) and (13), we finally obtain

$$P_{\nu _e\nu _s}0.21\pm 0.07.$$
(15)
The relatively large transition probability (15) is compatible with the fact that reactor experiments (see Ref. ) did not measure a disappearance of electron antineutrinos if the mass-squared difference that drives the NuTeV oscillations is222 Our results may be compatible with those presented a long time ago in Ref. .

$$\mathrm{\Delta }m^210100\mathrm{eV}^2,$$
(16)
corresponding to an oscillation length in reactor experiments

$$L_{\mathrm{osc}}^{\mathrm{reactors}}=\frac{4\pi E_{\mathrm{reactors}}}{\mathrm{\Delta }m^2}110\mathrm{cm},$$
(17)
where $`E_{\mathrm{reactors}}1\mathrm{MeV}`$. Since the source-detector distance in reactor neutrino oscillation experiments is larger than about 10 m, the average survival probability of electron antineutrinos is a constant and cannot be measured without assuming a precise knowledge of the initial $`\overline{\nu }_e`$ flux.

The oscillation length corresponding to a $`\mathrm{\Delta }m^2`$ in the range (16) for an energy of 100 GeV, typical of the NuTeV experiment, is

$$L_{\mathrm{osc}}^{\mathrm{NuTeV}}110\mathrm{km},$$
(18)
which is appropriate for the observation of transitions with the NuTeV source-detector distance of 1.5 km.

Disappearance of electron neutrinos driven by the $`\mathrm{\Delta }m^2`$ given in Eq. (16) is particularly attractive for measurements in the near future using the novel beta beam technique proposed in Ref. , with a source-detector distance of the order of 10 m for a neutrino energy $`E0.11\mathrm{GeV}`$.

These $`\nu _e\nu _s`$ transitions could also be observed in other high precision accelerator short-baseline neutrino oscillation experiments like NOMAD , with a source-detector distance of the order of 1 km for a neutrino energy $`E10100\mathrm{GeV}`$.

A high precision measurement could be done in the long term at a neutrino factory (see Ref. ) with a short baseline detector located at a source-detector distance of 100 m with a neutrino energy $`E110\mathrm{GeV}`$ .

The $`\mathrm{\Delta }m^2`$ in Eq. (16) could be the same that drives the small $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations observed in the LSND experiment . In this case, $`\nu _\mu \nu _e`$ appearance will be observed in the coming MiniBooNE experiment . Simultaneous $`\nu _e\nu _s`$ disappearance and $`\overline{\nu }_\mu \overline{\nu }_e`$ appearance from neutrinos produced in $`\mu ^+\overline{\nu }_\mu +e^++\nu _e`$ decays333 Or $`\overline{\nu }_e\overline{\nu }_s`$ disappearance and $`\nu _\mu \nu _e`$ appearance from neutrinos produced in $`\mu ^{}\nu _\mu +e^{}+\overline{\nu }_e`$ decays.  driven by the same $`\mathrm{\Delta }m^2`$ given in Eq. (16) could be easily observed at a neutrino factory with a short baseline detector with electron charge discrimination capability . The signature of $`\nu _e\nu _s`$ oscillations will be a depletion in the $`e^{}`$ energy spectrum induced by $`\nu _e`$’s and the presence of $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations will be characterized unequivocally by the appearance of $`e^+`$ events.

$`\nu _e\nu _s`$ transitions driven by the high $`\mathrm{\Delta }m^2`$ in Eq. (16) contribute with an energy-independent probability to the disappearance of solar electron neutrinos and must be taken into account in the analysis of solar neutrino data (see ). The energy dependence of the solar electron neutrino disappearance probability observed in solar neutrino experiments could be due to another very small squared-mass difference, in the framework of four-neutrino mixing schemes compatible also with the experimental evidence of atmospheric neutrino oscillations (see ).

The range (16) for $`\mathrm{\Delta }m^2`$ is compatible with the upper limit $`m_{\nu _e}3\mathrm{eV}`$ found in tritium decay experiments (see Ref. ). Indeed, the averaged probability (15) implies approximately

$$\mathrm{sin}^22\vartheta 0.4,$$
(19)
where $`\vartheta `$ is the mixing angle. Assuming a hierarchy of neutrino masses, the effective electron neutrino mass in tritium beta-decay experiments is

$$m_{\nu _e}\mathrm{sin}^2\vartheta \sqrt{\mathrm{\Delta }m^2}0.31\mathrm{eV},$$
(20)
that is compatible with the experimental upper bound.

If there is a hierarchy of neutrino masses, the effective neutrino mass in neutrinoless double-beta decay is equal to $`m_{\nu _e}`$ in Eq. (20). This value is compatible with the existing experimental upper limits  and curiously overlaps with the evidence in favor of neutrinoless double-beta decay claimed in Ref.  (see, however, the discussion in Refs. ).

The large neutrino mass implied by the $`\mathrm{\Delta }m^2`$ in Eq. (16) is compatible with the cosmological upper bound on light neutrino masses (about 10 eV), but would have important effects on the evolution of the universe (see Ref. ).

In conclusion, we have propose an explanation of the NuTeV in terms of $`\nu _e\nu _s`$ oscillations, with the relatively large average transition probability given in Eq. (15) and a large squared-mass difference given in Eq. (16). These transitions are compatible with the results of other neutrino experiments and could be verified by future high-precision experiments. Let us finally remark again that our analysis has been approximate and based on the limited information available in the paper of the NuTeV collaboration . We hope that the NuTeV collaboration will consider seriously our suggestion and perform a precise analysis of their data, which could lead to a determination of the allowed ranges of $`\mathrm{\Delta }m^2`$ and $`\mathrm{sin}^22\vartheta `$.