text:
# Supernova pointing with low- and high-energy neutrino detectors

## I Introduction

Observing a galactic supernova (SN) is the holy grail of low-energy neutrino astronomy. The question “how well can one locate the SN in the sky by the neutrinos alone?” is important for two reasons. Firstly, the MeV-neutrino burst precedes the optical explosion by several hours so that an early warning can be issued to the astronomical community Scholberg:1999tm ; snews2 , specifying the direction to look for the explosion. Secondly, in the absence of any SN observation in the electromagnetic spectrum, a reasonably accurate location in the sky is crucial for determining the neutrino Earth-crossing path to various detectors since the Earth matter effects on SN neutrino oscillations may well hold the key to identifying the neutrino mass hierarchy Dighe:1999bi ; Lunardini:2001pb ; Lunardini:2003eh ; Takahashi:2001dc ; Dighe2003 ; Dighe:2003jg .

Nearly contemporaneously with the optical explosion an outburst of TeV neutrinos is expected due to pion production by protons accelerated in the SN shock Waxman2001 . This neutrino burst could produce of order 100 events in a future km2 high-energy neutrino telescope, allowing for a pointing accuracy of a few tenths of a degree. Apart from the precise SN pointing, the detection of high-energy neutrinos would be important as the first proof that SN remnants accelerate protons.

The optical signal from a SN, if observed, can give the most accurate determination of its position in the sky. Apart from the observations at the optical telescopes, the multi-GeV to TeV photons associated with the accelerated protons in the shock could be detected on the ground by air Cherenkov telescopes after the SN environment becomes transparent to high-energy photons. If a suitable x- or $`\gamma `$-ray satellite is in operation at that time, the SN would be visible in these wavebands starting from the optical explosion. A satellite like INTEGRAL could resolve the SN with an angular resolution of 12 arc-minutes integral .

However, it is possible that the SN is not seen in the entire electromagnetic spectrum. This can be the case if it is optically obscured, no suitable x- or $`\gamma `$-ray satellite operates, and the air Cherenkov telescopes are blinded by daylight or the satellites and telescopes simply do not look in the right direction at the right time. It is also possible that not every stellar collapse produces an explosion so that only neutrinos and perhaps gravity waves can be observed. In such a scenario, the best way to locate a SN by its core-collapse neutrinos is through the directionality of $`\nu e^{}\nu e^{}`$ elastic scattering in a water Cherenkov detector such as Super-Kamiokande Beacom:1998fj ; sato . Much less sensitive methods include the time-of-arrival triangulation with several detectors Beacom:1998fj ; triang  or the systematic dislocation of neutrons in scintillation detectors that measure $`\overline{\nu }_epne^+`$ and the subsequent neutron capture Apollonio:1999jg .

The pointing accuracy of Super-Kamiokande or a future megaton detector such as Hyper-Kamiokande or UNO is strongly degraded by the inverse beta reactions $`\overline{\nu }_epne^+`$ that are nearly isotropic and about 30–40 times more frequent than the directional scattering events. Recently it was proposed to add to the water a small amount of gadolinium, an efficient neutron absorber, that would allow one to detect the neutrons and thus to tag the inverse beta reactions vagins . Evidently this would greatly improve the pointing: a tagging efficiency of 90% would double the pointing accuracy vagins . At high tagging efficiency, however, the nearly isotropic oxygen reaction $`\nu _e+{}_{}{}^{16}\mathrm{O}\mathrm{X}+e^{}`$ remains as the dominant background limiting the pointing accuracy. In this paper we analyze the realistic pointing accuracy of a water Cherenkov detector as a function of the neutron tagging efficiency.

The directionality of the elastic scattering reaction is primarily limited by the angular resolution of the detector and to a lesser degree by the kinematical deviation of the final-state positron direction from the initial neutrino. Extracting information from “directional data” is a field in its own right FLE ; Mardia . An efficient method is the “brute force” maximum likelihood estimate of the electron events, taking into account the angular resolution function of the detector on top of a nearly isotropic background. For a large number of events, the accuracy with this method in fact asymptotically approaches the minimum variance as given by the Rao-Cramér bound Rao ; cramer . However, for a small number of signal events, $`N_s200`$, we find that a fraction of the information content of the data as measured by the Fisher information fisher-paper  cannot be extracted by even the maximum likelihood method. Using the Rao-Cramér bound therefore overestimates the pointing accuracy of an experiment for small $`N_s`$. In this paper, we determine the realistic accuracy by using a concrete and nearly optimal estimation method.

Since the angular resolution depends on the event energy, the likelihood method requires as input the functional form of the neutrino energy spectra that are only poorly known. It is difficult to systematically take into account the errors introduced by a wrong choice of the fit function for the likelihood method. Therefore we look for “parameter-free” methods like harmonic analysis that use only the information contained in the data and exploit the symmetry of the physical situation. We discuss the efficiency of two methods closely related to the harmonic analysis and find a simple iterative procedure making them nearly as efficient as the maximum likelihood approach. We use the most efficient method thus obtained for analyzing the simulated events at a detector. This makes our estimations realistic and even a bit conservative, since the existence of a better parameter-free method is not excluded. We also study the dependence of the pointing accuracy on the neutrino mixing parameters and the initial neutrino spectra, and use the “worst case” scenario in order to estimate the accuracy.

Finally, we briefly study the pointing accuracy in high-energy neutrino telescopes that can pick up the TeV neutrino burst expected around the time of the optical explosion. While an accurate pointing is easy for any of the existing and future neutrino telescopes if the neutrinos are observed through the Earth, it is far more difficult against the atmospheric muon background from above. Even this would be possible for future km2 detectors such as IceCube or Nemo that could detect around $`100`$ SN events with TeV energies within about one hour.

We begin in Sec. II with a discussion of the statistical methodology for extracting information from directional data using a toy model. In Sec. III we study the realistic SN pointing accuracy of water Cherenkov detectors as a function of the neutron tagging efficiency, using realistic SN neutrino spectra. In Sec. IV we turn to the pointing accuracy of high-energy neutrino telescopes. Sec. V is given over to conclusions.

## II Analyzing directional data

### II.1 Pointing with maximum likelihood estimate

As a first case we study SN pointing with the $`\nu `$-$`e`$ elastic scattering events in the absence of any other background, thus obtaining a bound on the realistic pointing accuracy. To this end we use the toy model introduced in Ref. Beacom:1998fj , i.e. we imagine a directional signal that is distributed as a two-dimensional Gaussian on a sphere. This choice is motivated by the observation that the scatter of signal event directions is dominated by the angular resolution of the detector, and by the assumption that the angular resolution function is Gaussian. Later in Sec. III we consider a more realistic approximation to the angular detector response.

The angular width of the assumed Gaussian distribution is denoted by $`\delta _s`$, where $`s`$ stands for “signal.” As a further simplification we assume $`\delta _s\pi /2`$, allowing us to approximate the sphere by a plane. Taking the signal to be centered at $`\vartheta _0=0`$, the probability distribution function (pdf) of the signal events is

$$f_s(\vartheta ,\varphi )d\vartheta d\varphi =\frac{1}{𝒞}\mathrm{exp}\left(\frac{\vartheta ^2}{2\delta _s^2}\right)d\mu $$
(1)
with $`d\mu =\mathrm{sin}\vartheta d\vartheta d\varphi `$. Here $`𝒞𝑑\mu \mathrm{exp}[\vartheta ^2/(2\delta _s^2)]`$ is a normalization constant taking the value $`𝒞=2\pi \delta _s^2`$ for planar geometry.

In the case of Super-Kamiokande, around 300 elastic scattering events constitute the directional signal. Assuming the mean electron energy to be 11 MeV, a cone with opening angle $`\mathrm{}_{68}25^{}`$ around the true direction contains 68% of the reconstructed directions Nakahata:1998pz . Solving

$$_0^{2\pi }𝑑\varphi _0^{(\mathrm{}_{68}\pi /180^{})}𝑑\vartheta f_s(\vartheta ,\varphi )=0.68,$$
(2)
we find $`\delta _s17^{}`$. For $`N_s1`$ signal events without background, the Central Limit Theorem implies that the mean reconstructed direction is within $`\delta _s/\sqrt{N_s}`$ of the true direction for 68% of all SN obervations. This quantity, which is $`1^{}`$ for 300 events, gives the absolute lower bound on the pointing accuracy in the absence of all backgrounds. We note that Ref. Beacom:1998fj  uses $`\delta _s25^{}`$ and hence obtains $`1.5^{}`$ for the pointing accuracy. We further note that our $`\delta _s/\sqrt{N_s}1^{}`$ implies that 95% of all SN reconstructions lie within a circle of angular radius $`2.4^{}`$ of the true direction.

The main degradation of the pointing accuracy is caused by the nearly isotropic inverse beta decay background. The extent of this degradation was addressed for the first time in Ref. Beacom:1998fj . The pdf on a sphere that represents $`N_s`$ signal events distributed like a Gaussian around the direction $`(\vartheta _0,\varphi _0)`$ as well as the $`N_b`$ isotropic background events is

$`f(\vartheta ,\varphi |\vartheta _0,\varphi _0)d\vartheta d\varphi `$
$`={\displaystyle \frac{d\mu }{N_b+N_s}}\left[{\displaystyle \frac{N_b}{4\pi }}+{\displaystyle \frac{N_s}{𝒞}}\mathrm{exp}\left({\displaystyle \frac{\mathrm{}^2}{2\delta _s^2}}\right)\right],`$ (3)
where

$$\mathrm{}\mathrm{cos}^1[\mathrm{cos}\vartheta \mathrm{cos}\vartheta _0+\mathrm{sin}\vartheta \mathrm{sin}\vartheta _0\mathrm{cos}(\varphi \varphi _0)]$$
(4)
is the angular distance between the direction of an incoming neutrino and the experimentally measured direction of the Cherenkov cone. We have introduced here the usual notation $`f(x|x_0)`$ for the pdf to stress the dependence of $`f`$ on the data $`x=(\varphi ,\vartheta )`$ and the parameters $`x_0=(\varphi _0,\vartheta _0)`$.

The maximum likelihood estimate (MLE) method is the most efficient way to extract information from statistical data. For the pdf of Eq. (II.1) the likelihood function for $`N`$ events is

$$L(\vartheta _0,\varphi _0)\underset{\alpha =1}{\overset{N}{}}f(\vartheta ^{(\alpha )},\varphi ^{(\alpha )}|\vartheta _0,\varphi _0),$$
(5)
where $`(\vartheta ^{(\alpha )},\varphi ^{(\alpha )})`$ are the coordinates of the $`\alpha ^{\mathrm{th}}`$ event.

One commonly uses the Fisher information matrix fisher-paper  to estimate a lower bound on the uncertainty of the parameters extracted by the MLE method. In our case the Fisher matrix is defined as

$$F_{ij}\frac{^2\mathrm{ln}L(\vartheta _0,\varphi _0)}{\mathrm{\Theta }_i\mathrm{\Theta }_j},$$
(6)
where $`i,j=1,2`$, $`\mathrm{\Theta }_1\vartheta _0`$, $`\mathrm{\Theta }_2\varphi _0`$ and $`\mathrm{}`$ denotes an average with respect to $`fd\vartheta d\varphi `$. Since the off-diagonal elements of this matrix vanish, the error in the measurement of the two angles is simply given by

$$\mathrm{\Delta }\mathrm{\Theta }_i=\sqrt{1/F_{ii}}.$$
(7)
This lower bound on the pointing error is also known as the Rao-Cramér bound Rao ; cramer .

For the sake of definiteness, we have chosen the SN direction to lie in the equatorial plane so that $`\mathrm{\Delta }\vartheta _0=\mathrm{\Delta }\varphi _0`$, and define the pointing error $`\mathrm{\Delta }\vartheta `$ as

$$\mathrm{\Delta }\vartheta \sqrt{\frac{1}{N1}(\overline{\vartheta }\overline{\vartheta })^2},$$
(8)
where $`\overline{\vartheta }`$ is the estimate for $`\vartheta _0`$ from a given method and $`N`$ is the number of simulations used. The Rao-Cramér bound corresponds to the inequality

$$(\mathrm{\Delta }\vartheta )^2(\mathrm{\Delta }\vartheta )_{\mathrm{Fisher}}^21/F,$$
(9)
where $`F`$ may be $`F_{11}`$ or $`F_{22}`$.

The accuracy of the MLE method approaches asymptotically the Rao-Cramér bound for a sufficiently large number of events stat-book . For a finite number of data points, the MLE saturates this bound only if the pdf used in the MLE can be written in the product form

$$f(x|x_0)=g(x_0)h(x)\mathrm{exp}[A(x_0)B(x)].$$
(10)
If the pdf cannot be written in this form, the efficiency of the MLE is significantly smaller than unity at low values of $`N_s`$.

In Fig. 1 we show the efficiency $`\epsilon `$ of the MLE in “Fisher units”,

$$\epsilon (\mathrm{\Delta }\vartheta )_{\mathrm{Fisher}}^2/(\mathrm{\Delta }\vartheta )^2,$$
(11)
as a function of the number $`N_s`$ of signal events while keeping the background-to-signal ratio $`N_b/N_s`$ fixed at 30, which is the expected ratio of the inverse beta decay events to the elastic scattering events in a water Cherenkov detector. Though the MLE efficiency tends asymptotically to the Rao-Cramér bound for large values of $`N_s`$, this bound overestimates the MLE pointing accuracy by $`10`$% for $`N_s200`$.

We stress that the Rao-Cramér bound as an estimate of the MLE pointing efficiency is useful only in the limit where the spherical nature of our problem can be approximated by a planar geometry, i.e. when the angular resolution of the detector is much better than $`90^{}`$. This condition is satisfied in our case. For the general case of a pdf defined on a sphere the MLE still provides an optimal method to determine the pointing accuracy, but the Fisher information matrix no longer provides a direct asymptotic bound on the pointing accuracy. We are not aware of a generalized version of the Rao-Cramér bound that would relate directly to the variance of the pointing estimate in the case of a truly spherical problem.

### II.2 Efficiencies of parameter-free methods

The MLE is an optimal method to extract information from experimental data if the probability distribution function is known. This is not the case in our situation, where the exact forms of the neutrino spectra are needed and these are only poorly known. It is therefore worthwhile to look for other methods which may be less efficient, but which do not depend on the exact form of the pdf. In our case of SN pointing we wish to consider methods that do not depend on prior knowledge of the exact neutrino energy spectra.

Let us consider two pointing methods that exploit the symmetries of our physical situation, but are independent of the exact details of the pdf. If the efficiency of such a method turns out to be comparable to the maximum likelihood method for the toy model of the previous section, then that method may be expected to be efficient also in the realistic case where the exact pdf is not known and the MLE method cannot be employed.

One obvious approach is the “center of mass” (CM) method. The center of mass $`𝐒`$ of the events,

$$S_i\underset{\alpha =1}{\overset{N}{}}m^{(\alpha )}x_i^{(\alpha )},i=1,2,3,$$
(12)
is taken to be the estimator of the true center of the distribution, where $`𝐱(\mathrm{sin}\vartheta \mathrm{cos}\varphi ,\mathrm{sin}\vartheta \mathrm{sin}\varphi ,\mathrm{cos}\vartheta )`$. For an ideal detector the event weights $`m^{(\alpha )}`$ are equal and can be set to one. More realistically, the detection probability $`p`$ is a function of the detection angles $`\vartheta ,\varphi `$ and the weight is $`m^{(\alpha )}=1/p(\vartheta ^{(\alpha )},\varphi ^{(\alpha )})`$.

A second approach is the “orientation matrix” (OM) method FLE . Here, the major principal axis of the orientation matrix

$$T_{ij}\underset{\alpha =1}{\overset{N}{}}m^{(\alpha )}x_i^{(\alpha )}x_j^{(\alpha )}$$
(13)
is taken to be the estimator.

These two methods are equivalent to a harmonic analysis,

$$a_{lm}=\underset{\alpha =1}{\overset{N}{}}Y_{lm}(\vartheta ^{(\alpha )},\varphi ^{(\alpha )}),$$
(14)
restricted to the first moment for CM, and up to the second moment for OM. This equivalence can be seen either by a direct evaluation of Eq. (14) or by identifying $`S_i`$ with a dipole moment and relating $`T_{ij}`$ to a quadrupole moment $`Q_{ij}`$. Since $`3T_{ij}=Q_{ij}+N\delta _{ij}`$, the direction of the major principal axis is identical for both ellipsoids. Note that the orientation matrix is a reducible tensor and therefore contains information from the first as well as the second moment.

Neither of these methods requires any prior knowledge of the neutrino spectra or cross sections. However, they involve some loss of information and hence will give larger pointing errors than the MLE. In order to quantify the efficiency of these methods we generate a data sample according to the pdf of Eq. (3) and show the respective pointing errors in Fig. 2. We keep the number of signal events fixed at $`N_s=300`$, and show the pointing error $`\mathrm{\Delta }\vartheta `$ as a function of $`N_b/N_s`$. Note that in the absence of neutron tagging this ratio is expected to be around 30–40.

Figure 2 shows that the error of MLE is almost the same as the Rao-Cramér (RC) bound. However, the errors of CM and OM are much larger. One may also notice that OM is more accurate than CM. This difference may be attributed to the fact that whereas CM tries to exploit the spherical symmetry of the background, OM exploits the cylindrical symmetry of the background about the arrival direction, which is broken much more weakly by statistical fluctuations in the background. Moreover, in terms of a harmonic analysis, OM involves information from both $`l=1`$ and 2 while CM involves only $`l=1`$.

In order to increase the efficiency of CM and OM, we use the physical input that the signal is concentrated within a small region around the peak. Cutting off the events beyond a certain angular radius would then increase the signal to background ratio and the above methods may be applied iteratively to this new data. This procedure converges quickly and gives a much better estimate of the incoming neutrino direction. The optimal value of the angular cut has a very weak dependence on the number of events and the background-to-signal ratio. It depends mainly on the value of $`\delta _s`$ and is found to lie between $`2\delta _s`$ and $`3\delta _s`$. Within this range, the efficiency depends only weakly on the exact value of the angular cut. We tried both a sharp cutoff and a Gaussian weight function; both choices give practically identical results.

The optimal value of the cut also increases slowly with decreasing background-to-signal ratio, and in the limit of zero background, the method without cut is clearly more efficient than the method with cut since the latter now cuts off signal but no background. However, the variation due to changing the cut is of the order of only a few percent. Therefore, we keep the value of the cut to be constant at $`40^{}`$ for the analysis in this section. Our results for the pointing accuracy will therefore be somewhat conservative.

We denote the CM and OM methods with this cutting procedure by CMc and OMc, respectively. Figure 2 shows how the pointing error decreases drastically with the cutting procedure. It may also be observed that the accuracy of CMc stays close to that of MLE for low values of the background-to-signal ratio, while the accuracy of OMc is close to MLE in the entire $`N_b/N_s`$ range.

The efficiencies of all methods depend on both $`N_s`$ and $`N_b/N_s`$. In Fig. 1 we show the efficiency $`\epsilon `$ for different methods as a function of $`N_s`$, keeping $`N_b/N_s`$ fixed at 30. For a large number of signal events, $`N_s300`$, all methods tend to their asymptotic efficiencies $`\epsilon _{\mathrm{}}`$. The OMc method is close to its asymptotic efficiency $`\epsilon _{\mathrm{}}0.77`$ already for $`N_s200`$ whereas CMc needs $`N_s300`$ events to reach $`\epsilon _{\mathrm{}}0.71`$. Since the OMc method turns out to be more efficient that CMc in all the parameter ranges, henceforth we continue using only the OMc method for further estimations.

With neutron tagging, the value of $`N_b/N_s`$ decreases, and that increases the efficiency of the OMc method. In Fig. 3, we show the efficiency of OMc at different values of $`N_b/N_s`$. At $`N_b/N_s=10`$, which corresponds to the tagging efficiency $`\epsilon _{\mathrm{tag}}67\%`$, the asymptotic efficiency of OMc already increases to 0.85, and at $`N_b/N_s=5`$, corresponding to $`\epsilon _{\mathrm{tag}}83\%`$, it reaches 0.90. Moreover, with decreasing $`N_b/N_s`$ the asymptotic value is reached at lower and lower number of signal events. For higher tagging efficiencies, the optimal value of the angular cut increases. In fact, as noted before, in the limit of no background the OM method without the cut is more optimal. However this limit is physically not reached due to the presence of oxygen events.

The OMc method thus sacrifices less than 25%, and at higher tagging efficiency, even less than 15%, of the pointing accuracy of the MLE method. On the other hand, it has the great advantage of being independent of the detailed neutrino spectra and cross sections. Therefore, this method can be extremely useful for a fast analysis of the SN signal. After all, an early warning would depend on a quick and simple data analysis while later one can certainly optimize by fitting detailed energy spectra to the observed signal.

We stress that our preference for a parameter-free method over MLE in this analysis is strongly influenced by the current status of our knowledge regarding the pdf of the angular distribution. It may indeed be possible to use MLE and include all the systematic uncertainties, perhaps giving a better estimate for the pointing accuracy. However, faced with a tradeoff between model independence and higher efficiency, we give more weight to the former. If in future we understand the primary spectra much better than we do now, this preference may change.

## III Supernova pointing accuracy of water Cherenkov detectors

We now apply our method to a more realistic representation of the SN signal in a water Cherenkov detector. We shall limit our analysis of the pointing accuracy to our best parameter-free method, i.e. the Orientation Matrix method with an angular cut (OMc).

In order to determine the pointing accuracy numerically we simulate a large ensemble of SN signals in a water Cherenkov detector, assuming different efficiencies for neutron tagging. To this end we assume that the SN is at a distance $`D=10`$ kpc and releases the neutron-star binding energy $`E_b=3\times 10^{53}`$ erg in the form of neutrinos. Details of the assumed neutrino spectra and fluxes are given in Appendix A. The spread in the predicted neutrino spectra has been taken care of by using two models, a model from the Garching group (model G) garching-model  and a model from the Livermore group (model L) livermore-model  as described in the same appendix. We take into account the effects of neutrino flavor conversions by considering the three mixing scenarios, (a) normal mass hierarchy and $`\mathrm{sin}^2\mathrm{\Theta }_{13}10^3`$, (b) inverted mass hierarchy and $`\mathrm{sin}^2\mathrm{\Theta }_{13}10^3`$, and (c) any mass hierarchy and $`\mathrm{sin}^2\mathrm{\Theta }_{13}10^3`$. The six combinations of the models and neutrino mixing scenarios are represented by G-a, G-b, G-c, L-a, L-b, L-c. We use $`\mathrm{sin}^2(2\mathrm{\Theta }_{})=0.9`$ for the solar neutrino mixing angle.

As reaction channels we use elastic scattering on electrons $`\nu e^{}\nu e^{}`$, inverse beta decay $`\overline{\nu }_epne^+`$, and the charged-current reaction $`\nu _e+{}_{}{}^{16}\mathrm{O}X+e^{}`$, while neglecting the other, subdominant reactions on oxygen. The cross sections for these reactions are summarized in Appendix B. The oxygen reaction is included because it provides the dominant background for the directional electron scattering reaction in a detector configuration with neutron tagging where the inverse beta reaction can be rejected.

For the detector we assume perfect efficiency above an “analysis threshold” of 7 MeV, and a vanishing efficiency below this energy. The actual detector threshold may be as low as 5 MeV. Though lowering the threshold increases the ratio of elastic scattering events and the inverse beta events, it also introduces a background from the neutral-current excitations of oxygen (see Appendix B). In order to avoid additional uncertainties from the cross section of these oxygen reactions, we use the higher analysis threshold. We have checked that the net improvement by lowering the threshold to 5 MeV is less than 10% in all cases.

We assume a fiducial detector mass of 32 kiloton of water. Using the neutrino spectra and mixing parameters from the six cases mentioned above, we obtain 250–300 electron scattering events, 7000–11500 inverse beta decays, and 150–800 oxygen events. The ranges correspond to the variation due to the six different combinations of neutrino mixing scenarios and models for the initial spectra.

The procedure of event generation is described in Appendix C. The angular resolution function of the Super-Kamiokande detector does not follow a Gaussian distribution, rather it is close to a Landau distribution that we use for our simulation. In Fig. 4, the angular distribution of $`\overline{\nu }_epne^+`$ events (green) and elastic scattering events $`\nu e^{}\nu e^{}`$ (blue) of one of the simulated SNe are shown in Hammer-Aitoff projection, which is an area preserving map from a sphere to a plane.

The position of the SN is estimated with the OMc method. As explained in Sec. II, the optimal value of the angular cut depends on the neutron tagging efficiency as well as the neutrino spectra. We use a sharp cutoff with $`30^{}`$ opening angle for the OMc, which may not be optimal, but is observed to be close to optimal in almost the whole parameter range. For low values of $`\epsilon _{\mathrm{tag}}`$, the value of the cut should be lowered whereas for large values of $`\epsilon _{\mathrm{tag}}`$ it should be increased by about $`10^{}`$. The optimal cut depends also on the details of the detector properties and neutrino spectra.

A histogram of the angular distances between the true and the estimated SN position found in $`40000`$ simulated SNe for different neutron tagging efficiencies for the case G-a is shown in Fig. 5. The histogram fits well the distribution

$$f(\mathrm{})d\mathrm{}=\frac{1}{\delta ^2}\mathrm{exp}\left(\frac{\mathrm{}^2}{2\delta ^2}\right)\mathrm{}d\mathrm{},$$
(15)
where $`\mathrm{}`$ is the angle between the actual and the estimated SN direction, and $`\delta `$ is a fit parameter.

Defining the opening angle $`\mathrm{}_\alpha `$ for a given confidence level $`\alpha `$ as the value of $`\mathrm{}`$ for which the SN direction estimated by a fraction $`\alpha `$ of all the experiments is contained within a cone of opening angle $`\mathrm{}`$, we show in Fig. 6 the opening angle for 95% C. L. for the six cases of neutrino parameters. Clearly, the pointing accuracy depends weakly on the neutrino mixing scenario as well as the initial neutrino spectra. Some salient features of this dependence may be understood qualitatively as follows.

The signal events are dominated by $`\nu _e`$. Indeed, nearly half of the elastic scattering events are due to $`\nu _e`$, whereas the remaining half are due to the other five neutrino species. The cross section of electron scattering events increases with energy. Therefore, the more energetic the $`\nu _e`$ arriving at the detector, the larger the number of signal events and the better the pointing accuracy. Though the initial average $`\nu _e`$ energies are equal in the models G and L, the model L gives a much larger average energy for the initial $`\nu _x`$ spectrum. The $`\nu _e`$-$`\nu _x`$ mixing then tends to give more energetic $`\nu _e`$ in the model L. As a result, for each mixing scenario, the model L predicts a better pointing accuracy than the model G.

Within a model, the pointing accuracy is governed by the background-to-signal ratio. Since the cross section of the dominating background $`\overline{\nu }_ep`$ reaction increases with energy, more $`\overline{\nu }_e`$-$`\nu _x`$ mixing tends to give more energetic $`\overline{\nu }_e`$ and hence more background and less pointing accuracy. The ratio of $`\overline{\nu }_e`$-$`\nu _x`$ mixing to $`\nu _e`$-$`\nu _x`$ mixing within a model is the smallest for the mixing scenario (a) and the largest for the scenario (b). Therefore, within a model the scenario (a) always gives the best pointing accuracy and the scenario (b) always gives the worst. Note that in the limit $`\epsilon _{\mathrm{tag}}=1.0`$ when all the $`\overline{\nu }_ep`$ background is eliminated, the scenarios (b) and (c) give identical pointing accuracies since the final $`\nu _e`$ spectra in these two schemes are identical.

If neutrino experiments or supernova simulations have already identified the actual scenario, we could use the confidence limits given by that particular scenario. Indeed, the information for identifying the mixing scenario may already be contained in the observed neutrino spectra themselves which may be extracted by further data analysis Dighe:1999bi ; Lunardini:2001pb ; Lunardini:2003eh ; Takahashi:2001dc ; Dighe2003 ; Dighe:2003jg . However, in the absence of the immediate availability of this information, one has to select the least efficient scenario, G-b, in order to obtain the most conservative limits. This is the scenario with the least $`\nu _e`$-$`\nu _x`$ mixing, the largest $`\overline{\nu }_e`$-$`\nu _x`$ mixing and the lowest initial average energy for $`\nu _x`$.

Since the pointing accuracy is worse when the component of initial $`\nu _x`$ flux in the final $`\nu _e`$ flux is smaller, the “worst case” scenario is expected to be when the $`\nu _e`$ survival probability $`p`$ (see Appendix A) takes its highest possible value, which is nearly 0.45 at 3$`\sigma `$ fit-paper  , corresponding to $`\mathrm{sin}^2(2\mathrm{\Theta }_{})0.99`$. In Table 1 we give the values of the opening angle $`\mathrm{}`$ for this “worst case” scenario for various confidence levels $`\alpha `$ and tagging efficiencies. We also give the values of the fit for $`\delta `$ in Eq. (15), from which the numbers for any confidence level can be read off. For $`\epsilon _{\mathrm{tag}}=0`$, at 95% C.L. the pointing accuracy is $`7.8^{}`$, which improves to $`3.6^{}`$ for $`\epsilon _{\mathrm{tag}}=80\%`$ and $`3^{}`$ for $`\epsilon _{\mathrm{tag}}=1`$. This is nearly a factor of 3 improvement in the pointing angle, which corresponds to almost an order of magnitude improvement in the area of the sky in which the SN is located.

The last column (marked by $``$) in the table shows the pointing accuracy in the limit of no background, i.e. the case where all the inverse beta decay as well as the oxygen events are weeded out. This gives the intrinsic limiting accuracy due to the angle of electron scattering, the angular resolution of the detector and the efficiency of our OMc algorithm. For 95% C.L. the table gives a pointing accuracy of $`2.9^{}`$. This may be compared with the $`2.4^{}`$ that was estimated in Sec. II.1 for the toy model when there was no background. The degradation in the pointing accuracy may be attributed to the loss of information in the OMc method, the 10% smaller number of events in the SN simulation, and the difference between the actual angular distribution and the Gaussian as taken in the toy model.

For a SN at 10 kpc, in the worst case scenario we get nearly 10600 events, out of which the electron scattering signal is $`N_s270`$. We are then already at or just below the asymptotic limit for the efficiency of the OMc method. For a larger detector like Hyper-Kamiokande, the desired accuracy can be calculated simply by rescaling according to the number of signal events. For a detector with 25 times the fiducial volume of Super-Kamiokande, the pointing accuracy is then expected to be $`2^{}`$ without gadolinium and $`0.6^{}`$ with $`\epsilon _{\mathrm{tag}}>90\%`$. If the total number of events is much smaller than 10000, for small $`\epsilon _{\mathrm{tag}}`$ the efficiency of OMc will be smaller and this factor needs to be taken into account for calculating the real pointing accuracy. For $`\epsilon _{\mathrm{tag}}>0.8`$, the scaling with number of events should even work for a very small number of events, as can be seen from Fig. 3.

Note that all the model dependences discussed here, though clearly observable, are only around 10% (see Fig. 6). This indicates that the pointing accuracy estimates are quite model independent, and hence robust.

## IV High-energy neutrino telescopes

### IV.1 High-energy supernova neutrinos

Turning to the putative TeV neutrino burst associated with a SN explosion we note that the shock wave may well accelerate protons to energies up to $`10^{16}`$ eV. This idea is supported by the fact that the cosmic rays below the knee, $`E<10^{16}`$ eV, contain a total amount of energy comparable to that injected by all galactic SNe. The possibility of detecting high-energy neutrinos from a galactic SN has been discussed in the literature. In particular, it was shown that during the first year after the explosion, the SN shock wave will produce a large flux of neutrinos with energies above 100 GeV, inducing more than $`10^3`$ muons in a km2 detector BerezinskyPtuskin . More recently it was suggested that the high-energy neutrino signal would arrive just 12 hours after the SN explosion and would last for about one hour, giving about 100 muon events with $`E>1`$ TeV in a km2 detector Waxman2001 .

Of course, the number of expected events strongly depends on unknown parameters, in particular on the total energy emitted in the form of pions $`E_\pi ^{\mathrm{tot}}`$ at a given time and the maximum energy of the emitted neutrinos $`E_{\mathrm{max}}`$. The neutrino spectrum depends on that of the accelerated protons. However, higher-energy neutrinos have a greater chance of being detected so that we are not interested in the exact spectral shape once it is a typical power-law for shock acceleration, $`dN_\nu /dEE^\alpha `$ with $`\alpha 2`$. If the spectrum is softer, $`\alpha >2`$, detecting the neutrinos is more difficult. Note that the spectrum of protons after shock acceleration can be dominated by high-energy particles in some cases Derishev2003 , or even can be monochromatic if protons are accelerated in a potential gap. We will concentrate on the “standard” case of a $`E^2`$ neutrino spectrum, although our results for other cases will be similar. Following the calculation of Ref. Waxman2001  and assuming $`E_{\mathrm{max}}=1`$ TeV we expect around 50 muon events in a km2 detector during 1 hour at a time about 12 hours after the SN explosion. For a larger maximum energy, $`E_{\mathrm{max}}=1000`$ TeV, the number of muons increases to 200.

### IV.2 Signal from below

In order to discuss the expected signal in different neutrino telescopes Spiering2002 ; Halzen2002  we begin with the case where the SN happens in a part of the sky that a given neutrino telescope sees through the Earth. Future km2 detectors like IceCube at the South Pole ICECUBE , or northern projects like the Gigaton Water Detector at Baikal GVD  and Nemo in the Mediterranean NEMO  can detect the high energy neutrinos. For each event the angular resolution is around one degree. In this case the pointing to the SN can be resolved with an accuracy of about $`1^{}/\sqrt{50}8^{}`$ (arc-minutes). However, this purely statistical error does not include possible systematic effects. Most important is the limited knowledge of the alignment of these detectors. Therefore, the pointing accuracy of a km2 detector is probably larger and around a few tenths of a degree.

Even the existing smaller detectors can see a significant signal. The northern sky is under control of AMANDA-II with an effective area of 0.1 km2 and angular resolution of $`2^{}`$ at TeV energies. AMANDA-II will then detect 5 (20) events for $`E_{\mathrm{max}}=1`$ TeV (1000 TeV) and thus will be able to resolve the SN direction to better than $`1^{}`$. After their completion, the northern projects ANTARES antares  and NESTOR nestor  will be comparable to AMANDA-II.

### IV.3 Signal from above

If the high-energy SN neutrinos arrive “from above,” they are masked by the large background of atmospheric muons. For IceCube this background is about $`5.2\times 10^{10}\mathrm{year}^1`$ from the upper hemisphere at the “trigger” level ice-trigger . This corresponds to nearly $`300\mathrm{hour}^1\mathrm{degree}^2`$. If we note that the angular resolution is about $`1^{}`$, the expected signal of 100 events will be much larger than the background fluctuations in one pixel of the sky. Morever, the expected SN neutrinos will have multi-TeV energies so that energy cuts will reduce the background. Also, the significance of the signal will be enhanced by the angular prior defined by the low-energy signal in Super-Kamiokande.

For AMANDA-II size detectors both background and signal are about 10 times smaller. Moreover, the angular resolution is only about $`2^{}`$. Therefore, the expected signal in one pixel of the sky will be comparable to the background fluctuations. Energy cuts may allow one to detect the SN signal from the “bad” side of the sky.

## V Summary and Conclusions

The MeV neutrinos from the cooling phase of a SN will arrive at the Earth several hours before the optical explosion. These neutrinos will not only give an early warning of the advent of a SN explosion, but they can also be used to determine the location of the SN in the sky, so that the optical telescopes may concentrate on a small area for the observation.

In a water Cherenkov detector like Super-Kamiokande, the $`\nu `$-$`e`$ scattering events are forward peaked and thus can be used for the pointing. The main background comes from the inverse beta decay reactions $`\overline{\nu }_epne^+`$, which is nearly isotropic and has a strength more than 30 times that of the electron scattering “signal” reaction. The reactions of the neutrinos with oxygen also contribute to the nearly isotropic background.

The authors of Ref. Beacom:1998fj  have estimated the pointing accuracy using the Rao-Cramér bound. This may overestimate the accuracy by $`10`$% since even the most efficient method, the maximum likelihood estimate, can reach the minimum variance bound only for a large number of signal events. More importantly, the maximum likelihood method needs as an input the exact form of the fit function, which is not available due to our currently poor knowledge of the neutrino spectra. Taking into account the errors due to a wrong choice of the fit function is difficult. Therefore we choose to calculate the pointing accuracy using a concrete and simple estimation method that is independent of the form of the fit function.

We explore some parameter-free methods that only use the data and exploit the symmetries inherent in the physical situation, and therefore give a model independent estimation of the pointing accuracy while sacrificing some information from the data. We perform a statistical analysis of these methods using a toy model, which has a Gaussian signal on top of an isotropic background. We find that a method that uses the “orientation matrix” with an appropriate angular cut (OMc) is an efficient method that uses more than 75% of the information contained in the data if $`N_s200`$. We argue that this loss of information is well worth the gain of model independence, and continue to use this method for the simulation of the actual angular distribution at a water Cherenkov detector. It turns out that this method is much more efficient than the one used in Ref. sato .

One may add gadolinium to Super-Kamiokande in order to tag neutrons and therefore reduce the background due to inverse beta decay events. We quantify the increase in the pointing accuracy as a function of the neutron tagging efficiency $`\epsilon _{\mathrm{tag}}`$. It is found that the accuracy increases by more than a factor of two with $`\epsilon _{\mathrm{tag}}=0.8`$ and by nearly a factor of 3 for $`\epsilon _{\mathrm{tag}}=0.95`$. For $`\epsilon _{\mathrm{tag}}>0.95`$, the oxygen events act as the major background and that saturates the advantage of increasing the tagging efficiency beyond this value.

The efficiency of the OMc method improves with a smaller background-to-signal ratio, which makes this method even more useful at high tagging efficiencies. It is also observed that at higher $`\epsilon _{\mathrm{tag}}`$ this method attains its maximum efficiency level for a much smaller number of events. The optimal value of the cut increases with increasing $`\epsilon _{\mathrm{tag}}`$, though this dependence is weak and we perform our estimations with a fixed value of the cut, which is near optimal in the whole parameter range. Our estimations are therefore slightly conservative.

With a simulation of a water Cherenkov detector like Super-Kamiokande, we determine the pointing accuracy obtained from a SN at 10 kpc. The accuracy has a weak dependence on the neutrino mixing scenarios and the initial neutrino spectra. We find that the worst case scenario is when the detected $`\overline{\nu }_e`$ spectrum has the largest admixture of the initial $`\nu _x`$ spectrum, and the detected $`\nu _e`$ spectrum has the lowest admixture of the initial $`\nu _x`$ spectrum. This worst case turns out to be the one with the inverted neutrino mass hierarchy, $`\mathrm{sin}^2\mathrm{\Theta }_{13}10^3`$, and the largest possible solar mixing angle. The OMc method gives the pointing accuracy of $`7.8^{}`$ at 95% C.L. without neutron tagging. The accuracy improves to $`3.6^{}`$ at 80% tagging efficiency and to $`3.2^{}`$ at 95% tagging efficiency. Beyond this, the pointing accuracy saturates due to the presence of the oxygen events and the limited angular resolution of the detector. For a larger detector, the expected accuracy may be scaled according to the number of events. At a megaton detector like Hyper-Kamiokande, this gives an accuracy of $`0.6^{}`$ for $`\epsilon _{\mathrm{tag}}>0.9`$.

The SN shock wave may produce a TeV neutrino burst that arrives at the Earth within a day of the initial MeV neutrino signal. This can give about 100 events with $`E>1`$ TeV at a km2 detector like IceCube. Since the angular resolution of this detector is as good as 1, the SN may be located to an accuracy of a few tenths of a degree. The limiting factor here is the alignment error of these detectors. The time correlation and the directionality of the events allows IceCube to detect them even “from above”against the background of atmospheric muons.

###### Acknowledgements.

We thank Michael Altmann for suggesting the Landau distribution as model of the angular resolution function of Super-Kamiokande and Saunak Sen for clarifying some statistics issues. We also thank Mathias Keil and Christian Spiering for useful discussions and comments. This work was supported, in part, by the Deutsche Forschungsgemeinschaft under grant No. SFB-375 and by the European Science Foundation (ESF) under the Network Grant No. 86 Neutrino Astrophysics. M.K. and D.S. acknowledge support by an Emmy-Noether grant of the Deutsche Forschungsgemeinschaft, R.T. a Marie-Curie-Fellowship of the European Commission.

## Appendix A Neutrino fluxes

For the time-integrated neutrino fluxes we assume distributions of the form Keil:2002in

$$F^0=\frac{\mathrm{\Phi }_0}{E_0}\frac{(1+\alpha )^{1+\alpha }}{\mathrm{\Gamma }(1+\alpha )}\left(\frac{E}{E_0}\right)^\alpha \mathrm{exp}\left[(\alpha +1)\frac{E}{E_0}\right],$$
(16)
where $`F^0`$ denotes the flux of a neutrino species emitted by the SN scaled appropriately to the distance travelled from the SN to Earth. Here $`E_0`$ is the average energy and $`\alpha `$ a parameter that relates to the width of the spectrum and typically takes on values 2.5–5, depending on the flavor and the phase of neutrino emission. The values of the total flux $`\mathrm{\Phi }_0`$ and the spectral parameters $`\alpha `$ and $`E_0`$ are generally different for $`\nu _e`$, $`\overline{\nu }_e`$ and $`\nu _x`$, where $`\nu _x`$ stands for any of $`\nu _{\mu ,\tau }`$ or $`\overline{\nu }_{\mu ,\tau }`$.

We consider two models for the initial neutrino fluxes. The first one is the recent calculation from the Garching group garching-model , which we refer to as the model G. It includes all relevant neutrino interaction rates, including nucleon bremsstahlung, neutrino pair processes, weak magnetism, nucleon recoils and nuclear correlation effects. The second one is the Livermore simulation livermore-model , referred to as model L, that represents traditional predictions for flavor-dependent SN neutrino spectra that have been used in many previous analyses. The parameters of these models are shown in Table 2. We take $`\alpha (\nu _e)=\alpha (\overline{\nu }_e)=\alpha (\nu _x)=3.0`$ for both models. The value of $`\alpha `$ is not expected to have any significant influence on the pointing accuracy.

When neutrino mixing is taken into account, the fluxes arriving at a detector are

$`F_{\nu _e}`$ $`=`$ $`pF_{\nu _e}^0+(1p)F_{\nu _x}^0,`$ (17)
$`F_{\overline{\nu }_e}`$ $`=`$ $`\overline{p}F_{\overline{\nu }_e}^0+(1\overline{p})F_{\nu _x}^0,`$ (18)
$`4F_{\nu _x}`$ $`=`$ $`(1p)F_{\nu _e}^0+(1\overline{p})F_{\overline{\nu }_e}^0+(2+p+\overline{p})F_{\nu _x}^0.`$ (19)
Since the four neutrino species $`\nu _x`$ cannot be distinguished at the detectors, we only give the sum of their fluxes, $`4F_{\nu _x}`$. Here $`p`$ and $`\overline{p}`$ are the survival probablilities of $`\nu _e`$ and $`\overline{\nu }_e`$ respectively.

Depending on the mass hierarchy and the value of the mixing angle $`\mathrm{\Theta }_{13}`$, the survival probabilities $`p`$ and $`\overline{p}`$ belong to one of the three mixing scenarios shown in Table 3. We neglect the Earth matter effects and the details of the “transition” region around $`\mathrm{sin}^2\mathrm{\Theta }_{13}10^3`$ Dighe:1999bi ; Lunardini:2001pb . We also neglect terms of order $`(\mathrm{\Theta }_{13})^2`$ in $`p`$ and $`\overline{p}`$. Then $`p`$ and $`\overline{p}`$ depend only on the solar mixing angle as given in the table, and are independent of the values of solar and atmospheric mass squared differences as well as the atmospheric mixing angle.

## Appendix B Neutrino reactions in water

### B.1 Elastic Scattering on Electrons

The differential cross-section of the reaction $`\nu +e^{}\nu +e^{}`$ with $`\nu =\{\nu _e,\overline{\nu }_e,\nu _{\mu ,\tau },\overline{\nu }_{\mu ,\tau }\}`$ is given by

$$\frac{d\sigma }{dy}=\frac{G_F^2m_eE_\nu }{2\pi }\left[A+B(1y)^2C\frac{m_e}{E_\nu }y\right],$$
(20)
where $`y=E_e/E_\nu `$ is the energy fraction transfered to the electron, $`G_F`$ the Fermi constant, and $`m_e`$ the electron mass. The coefficients $`A`$, $`B`$ and $`C`$ differ for the four different reaction channels and are given in Table 4. The vector and axial-vector coupling constants have the usual values $`C_V=\frac{1}{2}+2\mathrm{sin}^2\mathrm{\Theta }_W`$ and $`C_A=\frac{1}{2}`$ with $`\mathrm{sin}^2\mathrm{\Theta }_W0.231`$.

The scattering angle $`\vartheta `$ between the incoming neutrino and final-state electron is implied by

$$y=\frac{2(m_e/E_\nu )\mathrm{cos}^2\vartheta }{(1+m_e/E_\nu )^2\mathrm{cos}^2\vartheta }.$$
(21)
### B.2 Inverse Beta Decay

For $`\overline{\nu }_epne^+`$ we use the differential and total cross-sections Eqs. (5–7) of Ref. Strumia:2003zx , including the leading QED radiative corrections.

### B.3 Oxygen as a Target

Another important reaction is the charged-current $`\nu _e`$ absorption on oxygen Haxton:kc ; Kolbe:2002gk . The dominant channels are

$`\nu _e+{}_{}{}^{16}\mathrm{O}`$ $``$ $`{}_{}{}^{15}\mathrm{O}+p+e^{},`$ (22)
$`\nu _e+{}_{}{}^{16}\mathrm{O}`$ $``$ $`{}_{}{}^{15}\mathrm{O}_{}^{}+p+\gamma +e^{},`$ (23)
$`\nu _e+{}_{}{}^{16}\mathrm{O}`$ $``$ $`{}_{}{}^{14}\mathrm{N}_{}^{}+p+p+e^{}.`$ (24)
While these reactions cause far fewer events in a water Cherenkov detector than inverse beta decay, they do not have final-state neutrons and thus cannot be tagged. Therefore, in a detector configuration with efficient neutron tagging, these reactions provide the dominant background to the directional electron scattering reactions.

The neutrino energy threshold in these reactions is approximately 15 MeV. The total cross section, summed over all channels, has been tabulated for the range $`15E_\nu 100`$ MeV Kolbe:2002gk . Directly above threshold the cross section is very small. We find that for $`25E_\nu 100`$ MeV the tabulated cross sections are nicely represented by the analytic fit

$`\sigma \left(\nu _e+{}_{}{}^{16}\mathrm{O}\mathrm{X}+e^{}\right)=`$
$`\mathrm{\hspace{0.17em}4.7}\times 10^{40}\mathrm{cm}^2\left[\left({\displaystyle \frac{E_\nu }{\mathrm{MeV}}}\right)^{1/4}15^{1/4}\right]^6.`$ (25)
For an accurate determination of the detector response one needs the differential distribution of final-state $`e^{}`$ energies and angular directions for a given incident $`E_\nu `$. Ref. Haxton:kc  provides extensive plots of such distributions after folding them with thermal $`E_\nu `$ distributions. This information is too indirect for our purposes. Therefore, we limit our investigation to a schematic implementation of this process where we assume that in every reaction the final-state energy is $`E_e=E_\nu 15`$ MeV. For the angular distribution we assume

$$\frac{d\sigma }{d\mathrm{cos}\vartheta }=1\frac{1+(E_e/25\mathrm{MeV})^4}{3+(E_e/25\mathrm{MeV})^4}\mathrm{cos}\vartheta ,$$
(26)
where $`\vartheta `$ is the angle between incident $`\nu _e`$ and final-state $`e^{}`$. This means that for small energies the angular distribution is proportional to $`1\frac{1}{3}\mathrm{cos}\vartheta `$ while for large energies it is $`1\mathrm{cos}\vartheta `$, i.e. it becomes more backward peaked for high energies. Our schematic approach roughly mimics the behavior shown in Ref. Haxton:kc . Since the oxygen cross section is very energy dependent, the contribution of this reaction to the pointing accuracy depends sensitively on the neutrino energy spectrum and is thus very uncertain anyway.

Another potentially important class of charged-current reactions is Haxton:kc ; Kolbe:2002gk

$$\overline{\nu }_e+{}_{}{}^{16}\mathrm{O}X+e^+.$$
(27)
However, the contribution to the detector signal is somewhat smaller than caused by the above $`\nu _e`$ reactions. Moreover, the $`\overline{\nu }_e`$ reactions typically involve final-state neutrons and thus are rejected by neutron tagging. One exception is

$$\overline{\nu }_e+{}_{}{}^{16}\mathrm{O}{}_{}{}^{16}\mathrm{N}+e^+,$$
(28)
but its contribution is small. Therefore, we neglect this entire class of reactions in our study.

Another class of reactions is the neutral-current excitation of oxygen Langanke:1995he

$$\nu +{}_{}{}^{16}\mathrm{O}\nu +X+\gamma .$$
(29)
Most of these reactions cannot be rejected by neutron tagging. However, the total cross section for neutral-current scattering, including the channels without final-state $`\gamma `$, is smaller than for the charged-current $`\nu _e`$ reaction Kolbe:2002gk . Moreover, the $`\gamma `$ energies are below 10 MeV, and most of them even below our analysis threshold of 7 MeV. Therefore, we also neglect this class of reactions.

## Appendix C Event generation

For the generation of each of the events the following steps are performed:

1. The energy $`E_\nu `$ of the reacting neutrino is chosen according to $`F_\nu (E_\nu )\sigma (E_\nu )`$. The energy $`E_e`$ and the scattering angle $`\vartheta `$ of the outgoing electron/positron is chosen according to the differential cross-section of the particular reaction.
2. The measured energy $`E_{\mathrm{det}}`$ of the scattered electron/positron is determined by adding Gaussian noise with variance $`\sigma (E_e)=\sqrt{E_eE_0}`$, where $`E_0=0.22`$ MeV. If $`E_{\mathrm{det}}<E_{\mathrm{th}}=7`$ MeV, the event is not used in the data analysis.
3. The measured position $`(\varphi ,\vartheta )`$ of the event is simulated according the angular resolution function of Super-Kamiokande.

The angular resolution $`(\mathrm{})d\mathrm{}`$ is defined as the probability that inside two cones with opening angles $`\mathrm{}`$ and $`\mathrm{}+d\mathrm{}`$ around the true direction the reconstructed direction is contained. Reference Nakahata:1998pz  gives numerical values for the opening angle $`\mathrm{}`$ of a cone around the true direction which contains 68% of the reconstructed directions as well as the values of $`(\mathrm{})`$ as a function of $`\mathrm{}`$ for various energies. An inspection by eye shows that $`(\mathrm{})`$ is characterized by a large tail and cannot be well fitted by a Gaussian distribution. Inspired by the Landau distribution for energy losses, we have found that $`(\mathrm{})`$ is well described by

$$(\mathrm{})=𝒞\mathrm{exp}\left(\frac{x+e^x}{2}\right)\mathrm{sin}(\mathrm{}),$$
(30)
where $`𝒞`$ is a normalization constant, and

$`x`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{\sigma }}a,`$ (31)
$`a`$ $`=`$ $`0.7E/\mathrm{MeV}3.7,`$ (32)
$`\vartheta _{\mathrm{max}}`$ $`=`$ $`38^{}\sqrt{\mathrm{MeV}/E},`$ (33)
$`\sigma `$ $`=`$ $`[\vartheta _{\mathrm{max}}+\mathrm{ln}(2\sqrt{3})]/a.`$ (34)
For illustration, we show in Fig. 7 the angular resolution $`(\mathrm{})`$ as a function of $`\mathrm{}`$ for positron energies $`E_e=6.8`$ and 16.1 MeV as implemented in our simulation, together with the measurements extracted from Nakahata:1998pz .


format: text

text:
# A Dynamical Systems Approach to an Inviscid and Thin Accretion Disc

## 1 Introduction

In astrophysical fluid dynamics, studies in accretion processes occupy a position of prominence. Of such processes, critical (transonic) flows — flows which are regular through a critical point — are of great importance . A classic example of such a flow, in the stationary spherically symmetric situation, is the Bondi flow . What is striking about the Bondi solution is that while it is virtually impossible to numerically generate it in the stationary limit of the hydrodynamic flow, it is easy to lock on to the Bondi solution when the temporal evolution of the flow is followed. This dynamic and non-perturbative selection mechanism of the transonic solution is in very favourable agreement with Bondi’s conjecture that it is the criterion of minimum energy that will make the transonic solution the favoured one .

The appeal of the spherically symmetric flow, however, is somewhat pedagogic, limited as it is by its not accounting for the fact that in a reasonably realistic situation, the infalling matter would be in possession of angular momentum — and hence the process of infall should lead to the formation of what is known as an accretion disc. A well understood and quite regularly invoked model of such a system is the inviscid and thin accretion disc . In this work, we show how such a physical system — studied in its stationary limit — gives rise to saddle point like behaviour for at least one of the critical points of the flow. The adverse implications regarding the realizability of a solution passing through such a point is then addressed by considering the system in a real non-linear time evolution. A simplified mathematical analog of such a situation shows that the evolution indeed selects the critical (transonic) solution. Finally, we discuss that in the vicinity of the critical point, the resulting perturbation equation arising out of a linear stability analysis of the stationary solutions, quite closely resembles the metric of an acoustic black hole, and that is probably the signal that trajectories pass through the critical point in a unidirectional fashion.

## 2 The stationary equations of the flow and its critical points

The flow system to be considered is an axisymmetric thin disc of a compressible astrophysical fluid, with a thickness $`H`$ . In the vertical direction the flow is considered to be in hydrostatic equilibrium. The governing equations for the flow should therefore be the equation of continuity and the equations for momentum balance in both the radial and the azimuthal directions. It is customary to write the stationary flow equations in the following notation .

* Equation of continuity:

$$\frac{d}{dR}(\rho RHv)=0$$
(1)
* Radial momentum balance equation:

$$v\frac{dv}{dR}+\frac{1}{\rho }\frac{dP}{dR}+\left(\mathrm{\Omega }_{K}^{}{}_{}{}^{2}\mathrm{\Omega }^2\right)R=0$$
(2)
* Angular momentum balance equation:

$$v\frac{d}{dR}(\mathrm{\Omega }R^2)=\frac{1}{\rho RH}\frac{d}{dR}\left(\frac{\alpha \rho c_{s}^{}{}_{}{}^{2}R^3H}{\mathrm{\Omega }_K}\frac{d\mathrm{\Omega }}{dR}\right)$$
(3)

In the above, $`v`$ is the radial velocity, $`R`$ is the radial distance, $`\mathrm{\Omega }`$ is the angular velocity, $`\mathrm{\Omega }_K`$ is the Keplarian angular velocity defined as $`\mathrm{\Omega }_{K}^{}{}_{}{}^{2}=GM/R^3`$, $`c_s`$ is the velocity of sound defined as $`c_{s}^{}{}_{}{}^{2}=P/\rho `$ and $`\alpha `$ is the effective viscosity of Shakura and Sunyaev .

It is a standard practice to make use of a general polytropic equation of state $`P=k\rho ^\gamma `$ where $`k`$ and $`\gamma `$ are constants, with $`\gamma `$ being the polytropic exponent , whose admissible range $`(1<\gamma <5/3)`$ is restricted by the isothermal limit and the adiabatic limit respectively. The condition of hydrostatic equilibrium along a direction perpendicular to the thin disc, allows for the use of the approximation $`(H/R)(c_s/v_K)`$, where $`v_K=R\mathrm{\Omega }_K`$ . In that case Eq. (1) leads to

$$c_{s}^{}{}_{}{}^{2n+1}vR^{5/2}=\mathrm{constant}$$
(4)
where $`n=(\gamma 1)^1`$. Combining Eqs. (2) and (4) yields,

$$\frac{d}{dR}(v^2)=\frac{2v^2}{R}\left[\frac{5c_{s}^{}{}_{}{}^{2}/(\gamma +1)(\mathrm{\Omega }_{K}^{}{}_{}{}^{2}\mathrm{\Omega }^2)R^2}{v^22c_{s}^{}{}_{}{}^{2}/(\gamma +1)}\right]$$
(5)
and a first integral of Eq. (3) can be written as

$$\alpha \frac{d}{dR}(\mathrm{\Omega }^2)=2\mathrm{\Omega }\left(\mathrm{\Omega }\frac{L}{R^2}\right)\frac{v\mathrm{\Omega }_K}{c_{s}^{}{}_{}{}^{2}}$$
(6)
where $`L`$ is a constant of integration.

We now consider the inviscid limit of Eq. (6) by setting $`\alpha =0`$. This, somewhat simplifying, prescription has found regular favour in accretion literature , and it gives the condition $`\mathrm{\Omega }R^2=L`$, where $`L`$, which can be physically identified as the specific angular momentum of the flow, now becomes a constant of the motion. Such a constraint allows us to fix the critical points of the flow. At such points both the numerator and the denominator in the right hand side of Eq. (5) vanish simultaneously , and this will deliver the critical point conditions as

$`v_{c}^{}{}_{}{}^{2}`$ $`=`$ $`{\displaystyle \frac{2c_{sc}^{}{}_{}{}^{2}}{\gamma +1}}`$
$`{\displaystyle \frac{5c_{sc}^{}{}_{}{}^{2}}{\gamma +1}}`$ $`=`$ $`(\mathrm{\Omega }_{Kc}^{}{}_{}{}^{2}\mathrm{\Omega }_{c}^{}{}_{}{}^{2})R_{c}^{}{}_{}{}^{2}`$ (7)
with the subscripted label $`c`$ indicating critical (fixed) point values. We should be able to find $`v_c`$ and $`R_c`$ from Eqs. (7). The second relation above, which is actually a quadratic in $`R_c`$, leads us to

$$R_c=\frac{\gamma +1}{10}\frac{GM}{c_{sc}^{}{}_{}{}^{2}}\left[1\pm \sqrt{1\frac{20}{\gamma +1}\left(\frac{c_{sc}L}{GM}\right)^2}\right]$$
(8)
and this gives the two critical points at radii $`R=R_{c1}`$ and $`R=R_{c2}`$, with $`R_{c2}>R_{c1}`$. To fix the critical point coordinates $`(R_c,v_c)`$, we need to look at the integral of Eq. (2), which reads as

$$\frac{v^2}{2}+\frac{L^2}{2R^2}\frac{GM}{R}+\frac{c_{s}^{}{}_{}{}^{2}}{\gamma 1}=\mathrm{constant}=\frac{c_{s}^{}{}_{}{}^{2}(\mathrm{})}{\gamma 1}$$
(9)
in which the constant of integration has been determined with the help of the boundary condition that for very great radial distances, the speed of sound approaches a constant “ambient” value $`c_s(\mathrm{})`$, while $`v0`$. In this expression we use the condition $`v_{c}^{}{}_{}{}^{2}=2c_{sc}^{}{}_{}{}^{2}/(\gamma +1)`$ and the value of $`R_c`$ from Eq. (8), to fix $`c_{sc}`$ in terms of the constants of the system. This will in turn fix the critical point coordinates $`(R_c,v_c)`$.

## 3 The disc as a dynamical system

To have any understanding of the nature of the critical points of the inviscid and thin disc being studied here, we need to cast the equation giving the stationary solutions of the disc system, along the lines of a dynamical system. To do that we appeal to Eq. (5) and parametrize it in the form

$`{\displaystyle \frac{d\left(v^2\right)}{d\tau }}`$ $`=`$ $`2v^2\left[{\displaystyle \frac{5c_{s}^{}{}_{}{}^{2}}{\gamma +1}}\left(\mathrm{\Omega }_{K}^{}{}_{}{}^{2}\mathrm{\Omega }^2\right)R^2\right]`$
$`{\displaystyle \frac{dR}{d\tau }}`$ $`=`$ $`R\left[v^2{\displaystyle \frac{2c_{s}^{}{}_{}{}^{2}}{\gamma +1}}\right]`$ (10)
It is to be noted that the parametrization has been carried out in an arbitrary mathematical parameter space $`\tau `$, and in such a space the two equations above represent a dynamical system. To analyse the nature of the fixed points, we expand and linearize about the fixed point coordinates $`(v_{c}^{}{}_{}{}^{2},R_c)`$ in Eqs. (10), and in terms of the perturbed quantities we obtain a set of linear equations given by

$`{\displaystyle \frac{d}{d\tau }}\left(\delta v^2\right)`$ $`=`$ $`2v_{c}^{}{}_{}{}^{2}\left[{\displaystyle \frac{5}{2}}{\displaystyle \frac{\gamma 1}{\gamma +1}}\delta v^2+\left({\displaystyle \frac{64\gamma }{\gamma +1}}{\displaystyle \frac{GM}{R_{c}^{}{}_{}{}^{2}}}{\displaystyle \frac{L^2}{R_{c}^{}{}_{}{}^{3}}}{\displaystyle \frac{73\gamma }{\gamma +1}}\right)\delta R\right]`$
$`{\displaystyle \frac{d}{d\tau }}\left(\delta R\right)`$ $`=`$ $`R_c\left[{\displaystyle \frac{2\gamma }{\gamma +1}}\delta v^2+5{\displaystyle \frac{\gamma 1}{\gamma +1}}v_{c}^{}{}_{}{}^{2}\delta R\right]`$ (11)
Use of solutions of the form $`e^{\lambda \tau }`$ would deliver the eigenvalues $`\lambda `$ — growth rates of $`\delta v^2`$ and $`\delta R`$ — as

$$\lambda ^2=4\frac{\left(5\gamma \right)}{\left(\gamma +1\right)^2}\frac{GM}{R_c}c_{sc}^{}{}_{}{}^{2}\left[\beta \frac{L^2}{L_{Kc}^{}{}_{}{}^{2}}\right]$$
(12)
where $`\beta =(53\gamma )/(5\gamma )`$ and the local Keplerian angular momentum is $`L_{Kc}=\sqrt{GMR_c}`$. It is now quite evident that if $`R_c`$ is the outer fixed point $`R_{c2}`$, then $`\lambda ^2`$ is positive and the fixed point is saddle type, while if $`R_c`$ is the inner point $`R_{c1}`$, then $`\lambda ^2`$ is negative and the fixed point is centre type. We show the trajectories in Fig.1, in which the arrows indicate the direction of flows.

The fixed point $`\mathrm{C}_2`$ at $`R=R_{c2}`$ is the saddle. In a manner of speaking, it may be dubbed the “sonic point” because one of the two trajectories passing through it is transonic, rising from subsonic values far from it to supersonic values when $`R<R_{c2}`$. This trajectory is shown as the heavy solid curve, and is the accretion flow. The heavy dashed curve is the wind. The above demonstration that the fixed point $`\mathrm{C}_2`$ is a saddle has adverse implications regarding the realizability of the transonic trajectory passing through it. Conventional wisdom about the nature of saddle type points  gives us to believe that if we go directly to the static limit of the dynamic evolution, then the transonic trajectory would not be realizable. However, it is generally believed that transonic trajectories do occur in astrophysical systems such as the one being studied here . Our contention is that it is the temporal evolution of such a system that selects the transonic trajectory. This selection mechanism is purely non-perturbative as opposed to the perturbative technique of a linear stability analysis of all physically feasible stationary inflow solutions. The latter method offers us no clue as to the choice of the system for any particular solution, since all stationary solutions have been found to be stable under the influence of a linearized perturbation — considered both as a travelling wave and a standing wave .

As an aside it may also be mentioned here that the nature of the outer fixed point $`\mathrm{C}_2`$ in Fig.1, has a bearing on the possible range of values that the constant specific angular momentum $`L`$ may be allowed. For solutions passing through the outer critical point, it can be well recognized that the flow would be sub-Keplerian . In that case, it may be concluded that the condition $`0<(L/L_{Kc2})^2<1`$ would hold good. In addition to that, from Eq. (12), the saddle type behaviour of the point $`\mathrm{C}_2`$ would also imply that there would be another upper bound on $`L`$, given by $`(L/L_{Kc2})^2<\beta `$. For the admissible range of the polytropic index $`\gamma `$, the possible range for $`\beta `$ would be $`0<\beta <1/2`$. This would then imply that the latter bound on $`L`$ would be more restrictive, as compared to the former. It is interesting that this essentially physical conclusion could be drawn from modelling the disc as a dynamical system in a mathematical parameter space.

## 4 A model for a dynamic selection mechanism

In the previous section we contended that the temporal evolution of the disc system selects a critical solution that passes through the outer critical point. To understand how this comes about, we consider a simplified mathematical example, since it is well known that the equations of a compressible flow cannot be exactly integrated.

The model system that we introduce, describes the dynamics of $`y(x,t)`$ as

$$\frac{y}{t}+\left(yx\right)\frac{y}{x}=y+2x+x^2$$
(13)
whose static limit yields

$$\frac{dy}{dx}=\frac{y+2x+x^2}{yx}$$
(14)
and which, viewed as dynamical system, is seen as

$`{\displaystyle \frac{dy}{d\tau }}`$ $`=`$ $`y+2x+x^2`$
$`{\displaystyle \frac{dx}{d\tau }}`$ $`=`$ $`yx`$ (15)
In the $`yx`$ space, the fixed points $`(x_c,y_c)`$ are to be found at $`(0,0)`$ and $`(3,3)`$. As in the preceding section, a linear stability analysis of the fixed points in $`\tau `$ space gives the eigenvalues $`\lambda `$ by $`\lambda ^2=1+2(x_c+1)`$. It is then easy to see that $`(0,0)`$ is a saddle point while $`(3,3)`$ is a centre type point. The integral curves are

$$y^22xy2x^2\frac{2}{3}x^3=c$$
(16)
and the trajectories passing through the saddle point are the ones with $`c=0`$. The different possible trajectories are shown in Fig.2. The separatrices are shown as the heavy solid curve and the heavy dashed curve. The similarity between Figs.1 & 2 is obvious. We now explore the dynamics.

To obtain a solution to Eq. (13), we need to apply the method of characteristics . This involves writing

$$\frac{dt}{1}=\frac{dx}{yx}=\frac{dy}{y+2x+x^2}$$
(17)
The task is to find two constants $`c_1`$ and $`c_2`$ from the above set and the general solution of Eq. (13) would be given by $`c_1=F(c_2)`$, where the function $`F`$ is to be determined from the initial conditions. It is easy to see that one of the constants of integration is clearly the $`c`$ of Eq. (16). Hence, we write $`c_1=c`$ and from the first part of Eq. (17), we can write

$$𝑑t=\pm \frac{dx}{\sqrt{3x^2+\left(2/3\right)x^3+c}}$$
(18)
which solves the problem in principle. To put this in a usable form, we need to carry out the integration of Eq. (18). This cannot be done exactly. For small $`x`$ (the most important region, since it is near the saddle), we can drop the $`x^3`$ term to a good approximation. Further, we choose only the positive sign in the right hand side of Eq. (18) by the physical argument that we would wish to evolve the system through a positive range of $`t`$ (time) values. Integration of Eq. (18) will then lead us to the result

$$\left(x+\sqrt{x^2+c/3}\right)e^{\sqrt{3}t}=c_2$$
(19)
which will then make the solution of Eq. (13) look like

$$y^22xy2x^2\frac{2}{3}x^3=F\left(\left[x+\sqrt{\frac{\left(yx\right)^2}{3}\frac{2}{9}x^3}\right]e^{\sqrt{3}t}\right)$$
(20)
We wish to study the evolution of the system from the initial condition $`y=0`$ for all $`x`$ at $`t=0`$. For small $`x`$, dropping the $`x^3`$ term again, allows us to determine the form of the function $`F`$ as $`F(z)=3(2\sqrt{3})z^2`$. The solution, consequently, becomes

$$y^22xy2x^2\frac{2}{3}x^3=3\left(2\sqrt{3}\right)\left[x+\sqrt{\frac{\left(yx\right)^2}{3}\frac{2}{9}x^3}\right]^2e^{2\sqrt{3}t}$$
(21)
and as $`t\mathrm{}`$, we select the separatrices given by

$$y^22xy2x^2\frac{2}{3}x^3=0$$
(22)
a result whose remarkable feature is worth stressing. The evolution started far away from the separatrix. In fact, it started as $`y=0`$ for all $`x`$ at $`t=0`$. The evolution proceeded through a myriad of stationary states (all arguably stable under a linear stability analysis) and then locked on to the separatrix. Our claim is that this dynamic mechanism selects the velocity profile in the practical situation as well, and this is entirely in consonance with a similar and well established selection mechanism for the spherically symmetric flow .

## 5 Concluding remarks

In a previous section we discussed that linear stability analysis of the stationary inflow solutions offers no direct clue as to the preference of the disc system for a particular stationary solution. However, in a somewhat indirect approach to this issue we could make use of a linear stability analysis to have a feel for the special feature of the critical solution. A very convenient way of carrying out the linear stability analysis is in terms of the variable $`\varphi =\rho ^{(\gamma +1)/2}vR^{5/2}`$ . In terms of this variable, the perturbation equation for $`\varphi `$ reads

$$\frac{^2\varphi ^{}}{t^2}+2\frac{}{R}\left(v_0\frac{\varphi ^{}}{t}\right)+\frac{1}{v_0}\frac{}{R}\left[v_0\left(v_{0}^{}{}_{}{}^{2}\frac{2}{\gamma +1}c_{s0}^{}{}_{}{}^{2}\right)\frac{\varphi ^{}}{R}\right]=0$$
(23)
in which $`\varphi ^{}`$ is the perturbation on the steady solution of $`\varphi `$, and the subscripted label $`0`$ indicates steady state values wherever in use .

In a spacetime with metric $`g^{\mu \nu }`$, the wave equation for a scalar variable is

$$\frac{}{x^\mu }\left(g^{\mu \nu }\frac{\varphi ^{}}{x^\nu }\right)=0$$
(24)
Ignoring the angle variables, and using the notation that $`\mu =1`$ corresponds to $`t`$ and $`\mu =2`$ corresponds to $`R`$, we have the identification

$$g^{11}=1,g^{12}=g^{21}=v_0,g^{22}=v_{0}^{}{}_{}{}^{2}\frac{2}{\gamma +1}c_{s0}^{}{}_{}{}^{2}$$
(25)
The effective covariant components are

$$g_{11}=1\frac{v_{0}^{}{}_{}{}^{2}}{\stackrel{~}{c_{s0}}^2},g_{12}=g_{21}=\frac{v_0}{\stackrel{~}{c_{s0}}^2},g_{22}=\frac{1}{\stackrel{~}{c_{s0}}^2}$$
(26)
in which $`\stackrel{~}{c_{s0}}=\sqrt{2}(\gamma +1)^{1/2}c_{s0}`$.

The differential distance in the metric space has the form

$$ds^2=\frac{1}{\stackrel{~}{c_{s0}}^2}\left[\stackrel{~}{c_{s0}}^2dt^2\left(dx^iv_{0}^{}{}_{}{}^{i}dt\right)\delta _{ij}\left(dx^jv_{0}^{}{}_{}{}^{j}dt\right)\right]$$
(27)
and this is to be compared with the Painlevé-Gullstrand form of $`ds^2`$ for a rotating black hole of mass $`\stackrel{~}{M}`$

$$ds^2=dt^2\left(dr\pm \sqrt{\frac{2G\stackrel{~}{M}}{r}\frac{L^2}{r^2}}dt\right)^2+r^2d\theta ^2$$
(28)
We now note from Eq. (7) that at the sonic point, $`v_{c}^{}{}_{}{}^{2}=(2/5)[(GM/R_c)(L/R_c)^2]`$, and hence in the vicinity of the sonic point the metric implied by the linear stability analysis of the stationary solution appears quite close in form to the metric of a rotating black hole. The unidirectional nature of trajectories near a black hole suggests that trajectories go through the sonic point, which is equivalent to saying that transonic trajectories will be chosen.

In concluding we would also like to comment on the nature of the possible critical point(s) when we consider the more likely case of a disc with a mechanism for the outward transport of angular momentum. It is possible with the help of Eq. (6) to expand and linearize about the critical point value of the angular velocity $`\mathrm{\Omega }`$, but without any knowledge of an exact dependence of $`\mathrm{\Omega }`$ on $`R`$, there would be no way of fixing the critical point(s). However it can be conceived that as an extension of the inviscid case, there may be multiple parameter-dependent critical points  and in that event at least one of them may be saddle type, with all its associated difficulties within the stationary framework itself. In that event, once again the selection mechanism has to be time evolutionary in nature.

## Acknowledgements

This research has made use of NASA’s Astrophysics Data System. One of the authors (AKR) would like to acknowledge the financial assistance given to him by the Council of Scientific and Industrial Research, Government of India.


text:
# Coronal element abundances of the post-common envelope binary V471 Tauri with ASCA

## 1 Introduction

Solar observations have indicated that coronal abundances are photospheric for high-FIP (First Ionization Potential) elements, but that the coronal plasma is over abundant in low-FIP elements (Feldman 1992). Since low-FIP elements are ionized in the solar chromosphere, whereas high-FIP elements are not (Geiss 1992), it is believed that some mechanism separates ions from neutrals within the chromosphere and transports them to the corona. Hénoux 1995 provides a review of suggested mechanisms.

CCD X-ray detectors and proportional counting spectrographs have revealed multi-temperature stellar coronae with a variety of abundance distributions over ionization potential. Some show a FIP effect similar to the sun, others show the opposite trend - abundance decreasing with FIP relative to the solar photosphere, i.e. an inverse FIP effect. However, because of the moderate spectral resolutions available, combined with the confusion of line blending, much has been made of the possible ambiguity in abundance measurements using spectral models of collisionally-excited plasmas (Drake 1996).

The grating spectrographs onboard EUVE (Drake, Laming and Widing 1996), and recently launched on the XMM-Newton and Chandra satellites (Brinkman et al. 2001; Audard et al. 2001) have confirmed the general results of lower-resolution spectroscopy. Consequently, despite the poorer resolution compared to grating cameras and the possibility of line confusion, confidence in CCD and proportional counter measurements has increased. In this paper we present ASCA spectroscopy of the K star in the close binary V471 Tau, in order to measure its coronal content and contrast this with photospheric abundances.

V471 Tau is a member of the Hyades open cluster. Photospheric abundances of the Hyades members have been measured by Cayrel et al. (1985) and Varenne & Monier (1999) and some specific abundances of V471 Tau by Martín et al. (1997). Accepting model assumptions, we show in this paper that coronal Fe is under-abundant relative to the mean photospheric Hyades content and Ne is over-abundant, at least relative to the solar photosphere. This may be indicating that there is an inverse-FIP effect in the upper atmosphere of the star.

V471 Tau is a 12.5 hour eclipsing binary with a white dwarf and tidally-locked K2 companion (Nelson and Young 1970). At a distance of 47 pc (Werner & Rauch 1997), it is the closest object which has been through a recent common envelope phase of evolution. During the red giant phase of the white dwarf progenitor’s life, the envelope of the giant was large enough to contain the orbit of the K star companion. While the two stars shared a common envelope, tides, friction and mass loss created a significantly smaller binary, with a short orbital period (Paczýnksi 1976). The envelope of the giant has since been ejected, but the binary still loses angular momentum through wind braking and gravitational radiation (Taam 1983). V471 Tau is therefore one of the best candidates for a pre-cataclysmic variable (pre-CV). CVs occur when the orbital period of a white dwarf–red dwarf binary becomes short enough for the main-sequence companion to fill it’s Roche lobe. Quasi-persistent accretion will then occur through Roche lobe overflow. These objects are the sources of dwarf nova outbursts and classical nova eruptions that enrich the galactic ISM (Warner 1995).

Soft X-rays from V471 Tau were discovered by HEAO 1 (van Buren, Charles & Mason 1980) and it was detected as a flaring source by the Einstein Observatory (Young et al. 1983). This suggested that the K star was the source of X-rays, but by determining the white dwarf temperature to be 35,000 K using the IUE satellite, Guinan & Sion (1984) showed that the compact object may also be a significant X-ray emitter. Jensen et al.  (1986) used EXOSAT to find white dwarf eclipses, pulsations, orbital dips and emission from both stellar components. Using the simultaneous EUV/X-ray data from the ROSAT all-sky survey, Barstow et al. (1992) determined that the 555 s pulsations were the result of accretion caps on the white dwarf rather than g-mode oscillations. The interaction of rotating magnetic fields from the white dwarf and K star result in non-thermal emission at radio energies between the two objects (Patterson, Caillault & Skillman 1993; Lim, White & Cully 1996; Nicholls & Storey 1999). Possibly it is this magnetically ejected material that is accreted upon the white dwarf. Wheatley (1998) employed ROSAT to separate the K star from the white dwarf spectrally and pointed out that the K star flux was brighter than the average Hyades member but similar to the younger, faster rotators of the Pleiades. Therefore coronal luminosity in this source is probably related to rotation rate rather than age. The white dwarf is softer spectrally than the ASCA energy band, therefore it is not a factor during the analysis of the current data.

## 2 Observations

The Advanced Satellite for Cosmology and Astrophysics (ASCA; Tanaka, Inoue & Holt 1994) comprised four identical telescopes of nested thin foil reflectors (Serlemitsos et al. 1995). Four cameras were mounted at the focal planes. Two had Solid state Imaging Spectrometers (SIS0 and SIS1; Burke et al. 1991), and two carried Gas Imaging Spectrometers (GIS2 and GIS3) which were scintillation proportional counters (Ohashi et al. 1991). The SIS cameras have an energy range 0.4–10 keV, a spectral resolution of $`E/\mathrm{\Delta }E`$ = 20 at 2 keV and a 22 arcmin field of view. The GIS instruments have an energy range 0.7–10 keV, a spectral resolution of $`E/\mathrm{\Delta }E`$ = 7.5 at 2 keV and a 50 arcmin field of view.

Between MJD 50321.75 and MJD 50322.90 (1996 Aug 26-7), ASCA was pointed at V471 Tau for a total exposure time of 50 ksec, over 17 near-contiguous orbits. The ASCA sequence number of the observations is 24032000. The SIS cameras operated in FAINT mode and the data converted to BRIGHT mode, while GIS data was taken in PH mode with normal bit assignments.

## 3 Analysis

Observations were screened using standard criteria and the FTOOLS software package v5.1. Events were filtered based on the satellite’s proximity to South Atlantic Anomaly passages, Earth elevation angle, pointing stability and rigidity of the geomagnetic field. Hot and flickering pixels were removed from the SIS events and we rejected particle events from the GIS observations. SIS data were additionally screened for proximity to the bright Earth-limb, time since the last passage over the day-night terminator, telemetry saturation and event rate. Only well-calibrated events, with grades 0, 2, 3 or 4 were retained. GIS data were additionally screened spatially to remove background ring and calibration source and further intervals of high particle background were clipped. A total of 41 ksec (SIS0), 44 ksec (SIS1), 48 ksec (GIS2) and 48 ksec (GIS3) of good data remained. The hot white dwarf component is softer spectrally than the ASCA energy band (Wheatley 1998), we find no evidence for it in the current data, and therefore we do not consider it in the following analysis.

### 3.1 Photometry

SIS event rates were determined within two spatially-distinct, circular region masks of diameter 8 arcmin and 5 arcmin for the source and background, respectively. Source events were extracted from the GIS tables within a circular region 12 arcmin in diameter and a background annulus 12–24 arcmin from the target. Background rates were normalized by the appropriate ratio of extraction areas and subtracted from the source events. The total SIS+GIS source event rate during the observation is provided in Fig. 1, averaged without weights into 1 ksec bins. We also provide the hardness ratio in the bands 0.3–1.0 keV and 1.0–10.0 keV.

Phases for each time bin were calculated using the orbital ephemeris of Guinan & Ribas (2001), converted to Modified Julian Date.

$$T=\text{MJD}\mathrm{40\hspace{0.17em}609.564056}(11)+0.521183398(26)E$$
(1)
where $`T`$ is the average MJD of each bin and $`E`$ is the cycle number. The orbital phase, $`\varphi `$, in Fig. 1 corresponds to $`E\mathrm{18\hspace{0.17em}624}`$. Quantities in parentheses correspond to $`1\sigma `$ uncertainties.

Count rate is approximately constant, except for three flares, each of duration a few ksec at orbital phases $`\varphi `$ 1.0, 2.4 and 2.9. Statistically, the shape of the flares are not well-defined and there is only marginal evidence for spectral hardening. Flares of similar duration and amplitude are common to the rapidly rotating stars in tidally-locked binaries, e.g., the RS CVn systems (Franciosini, Pallavicini & Tagliaferri 2003), and to single M dwarf stars (Tsikoudi & Kellett 2000). The observed flares do not occur at preferred orbital phases.

V471 Tau is relatively faint compared to most well-observed coronal sources (e.g. Brinkman et al. 2001, Audard et al. 2001), so in order to maximize counting statistics during spectral fitting, we do not separate flare events from the quiescent data.

### 3.2 Spectroscopy

To improve the accuracy of $`\chi ^2`$ fitting, the SIS and GIS energy channels containing less that 30 events were grouped with their neighbors. Channels flagged as bad were rejected and those with corrected pulse heights of $`0.3>E>10`$ keV were ignored. Best fits to the data were determined using the apec thermal plasma model (Brickhouse et al.  2000), absorbed through a neutral interstellar column with photo-electric cross-sections provided by Balucinska-Church & McCammon (1992). Spectra from all four cameras were employed together to converge on a single model.

The apec code (v1.10) contained within the current version of the spectral fitting package, xspec v11.2, does not include high $`n`$ transitions of Fe xvii–xix. Also there are systematic problems in line wavelengths for the set of Ni lines in this version. Therefore additional tables for Fe xvii; $`n`$ = 6, 7 and 8 $``$ 2, Fe xviii; $`n`$ = 6, 7 $`2`$ and Fe xix; $`n`$ = 6 $``$ 2 are included together with a table of corrected Ni wavelengths, as described in Brickhouse et al. (2000). Since these models are additive, abundances cannot be decoupled from model normalizations in these two cases. Although these models avoid errors in fitting high-$`n`$ Fe and Ni, they cannot be used to measure abundances.

If all abundances are allowed to roam as independent fit parameters, then a statistically significant constraint on any individual element is not possible using this data set. By trial and error we find that significant detections of Mg, Fe, Si, O and Ne (with H-like transitions at the peak of the ASCA response) are made if we adopt some reasonable assumptions concerning the fractions of the other elements contributing to the line spectrum, i.e., C, N, Al, S, Ar and Ca. Consequently we describe below four abundance models. Solar abundance values were taken from Grevesse & Sauval (1998).

Model 1 assumes the abundances of all metals are equal; in this case the abundance is dominated by the L shell lines of Fe. Model 2 allows the abundances of Mg, Si, O and Ne to float as free parameters while the others are floating but assumed equal. The combined metal abundance will again be dominated by Fe. Model 3 is similar to model 2 except that the Fe abundance is decoupled from the combined elements and allowed to float on its own. The remaining abundances of C, N, Al, S, Ar and Ca float, but are equal. Model 4 is identical to model 3 except C, N, Al, S, Ar and Ca are fixed abundances with the mean photospheric value of the Hyades members (\[m/H\] = $`+0.1`$; Martín, Pavlenko and Rebolo 1997, where \[m/H\] is the logarithmic ratio of the element abundance, m, and H abundance, relative to solar).

Realistically, we expect a whole range of plasma temperatures and, unsurprisingly, single-temperature models provide poor fits to the data with $`\chi ^2=571`$ for 337 d.o.f. (model 3). Two-temperature models yield statistically significant fits although they are still undoubtedly a simplification of the temperature structure. Table 1 lists the best fit parameters and $`\chi ^2`$ statistic for each model. All models are acceptable according to the $`\chi ^2`$ statistic so we cannot adopt a best-model and can generalize the results only.

Between Models 2, 3 and 4 we find general consistency between temperatures and emission measures. If the neutral hydrogen column density was left as a free parameter, in all cases it was found to be both consistent with zero and smaller than the total galactic column in that direction (1.58 $`\times \text{10}^{21}`$ cm-2; Dickey & Lockman 1990). However, since the ASCA energy range is not well suited for the measurement of the expected column in front of V471 Tau, we adopt the value derived from EXOSAT measurements by Jensen et al. (1986) of $`N_\mathrm{H}`$ = 5 $`\times \text{10}^{18}`$ cm-2. The flux from the source according to these models is 2.4 $`\times \text{10}^{12}`$ erg s-1 cm-2.

For each model, Fe abundance is most-likely sub-solar, the opposite to the solar corona, whereas Ne is greater than solar. In Fig. 2 we show the best-fit model 3, folded on the energy-resolved count rates from SIS0 and GIS2 and their residuals. Although far from conclusive, the high Ne/Fe ratio is suggestive of an inverse-FIP effect. For models 2, 3 and 4, the coronal abundances of the five significantly-constrained elements are plotted against FIP in Fig. 3.

## 4 Discussion

For various model assumptions, we have determined the element abundances of Mg, Fe, Si, O and Ne in the corona of V471 Tau and demonstrated a high Ne/Fe ratio which could be indicative of the inverse-FIP effect. In the following section, we ask whether we expect to find fossil abundances from the post common envelope epoch. We also discuss the proposed mechanisms behind the FIP effect and, finally, compare the coronal abundances with the measured photospheric abundances of dwarf stars in the Hyades cluster.

### 4.1 Common envelope abundances

During the white dwarf progenitor’s red giant phase of evolution, the two stars shared a common envelope and, therefore, there was element mixing between them. Using LTE atmosphere fits to optical spectroscopy of V471 Tau, and providing NLTE corrections, Martín et al. (1997) found photospheric abundances consistent with the Hyades, \[m/H\] = +0.1, in Al, Ca, Fe and Si. The one exception was Li, where \[Li/H\] = $`2.4\pm 0.3`$. This is greater than a factor $`10^2`$ more abundant than normal Hyades members. Although Li is usually depleted in convection layers, Martín et al.  argue that the rapid rotation of the K star maintains the fossil level of Li after sharing a common envelope phase with the Li-rich giant (de la Reza, Drake & da Silva 1996). Unfortunately there are no Li transitions in the ASCA energy band to test this independently.

Marks & Sarna (1998) calculate the evolution of photospheric abundances on the companion star of cataclysmic variables, including the effects of a common envelope epoch. They indicate that, indeed, photospheric abundances in the K star are expected to remain unchanged through the common envelope phase, and remain constant thereafter until the white dwarf accretion rate is sufficient to drive a cycle of nova eruptions. Therefore we do not expect to find evidence for the common envelope phase in the current measurements, however, peculiar abundances may yet be present if the proposed IUE detection of expanding nova material around V471 Tau by Bruhweiler & Sion (1986) and Sion et al.  (1989) is confirmed. However subsequent HST spectroscopy shows no indication of shell features but does detect transient coronal ejection events (Bond et al. 2001). These may well have been mistaken for an expanding shell before coronal ejections were time-resolved by HST. In light of both the observational and theoretical evidence, we will assume that the photospheric abundances of the companion star in V471 Tau are identical to the average Hyades dwarf. This is useful later for comparing the coronal and photospheric abundances.

### 4.2 FIP bias

Observations have indicated that the coronal, energetic particle and cosmic ray abundances in the sun are all different from the solar photospheric content (Meyer 1985a,b). This is probably equally true for active stars (Brinkman et al. 2001; Audard et al. 2003). While in most cases the photospheric abundances remain uncertain, there is at least one strong piece of evidence for fractionation regions in between stellar photospheres and coronae. It derives from the fact that coronal abundances are correlated with elemental FIP. In the solar atmosphere the ratio of coronal to photospheric abundances decreases with FIP. When comparing stellar coronal abundances to the solar photosphere some sources follow this trend while others show the “inverse FIP effect”, i.e. an increasing ratio with increasing FIP (Brinkman et al. 2001). The Ne/Fe ratio in V471 Tau is suggestive of an inverse-FIP trend.

FIP effects are thought to be directly associated with the process of element fractionation in stellar atmospheres (Hénoux 1995). The ionized fraction of each element will depend on the plasma temperature in the fractionation region. Some process is required to decouple the ionized and neutral plasmas. Possibly the charged particles are accelerated along field lines in loops above the stellar photosphere (Wang 1996). Flare ejection provide another mechanism to propel material into the corona (Schmelz 1993). The abundances of coronal elements as a function of FIP has been observed to change during X-ray flares. Audard, Güdel & Mewe (2001) find an inverse FIP affect in HR 1099 during quiescence but a quenching of this effect during flares.

### 4.3 Abundances in the Hyades

Abundances in most stellar photospheres are not well known. Instead, it is common practice to compare stellar coronal abundances with the solar photosphere. Rather than adopting solar photospheric abundances, ideally we would prefer to have direct measurements of photospheric abundances in the K companion of V471 Tau in order to make an unbiased comparison between photosphere and corona. These are poorly determined generally because of uncertainties in galactic element mixing and the details of stellar spectral synthesis models. Fortunately, though, V471 Tau is a member of the nearest open cluster, the Hyades (Werner & Rauch 1997), where photospheric abundances are considerably better determined than most coronal X-ray sources.

Metallicities in open clusters are indicators of the quantity and rate of chemical mixing in the galaxy. However stellar abundance determinations are generally plagued by a wide range of biases (see e.g. Griffin & Holweger 1989). Various techniques have consequently supplied abundances for the Hyades cluster ranging from \[Fe/H\] = $`0.09`$ (Tomkin & Lambert 1978) to \[Fe/H\] = $`+0.42`$ (Gustafsson & Nissen 1972).

Although not unequivocally free of bias, Cayrel, Cayrel de Strobel & Campbell (1985) determine a mean \[Fe/H\] enhancement over solar of $`+0.12\pm 0.03`$ using high dispersion spectroscopy of a sample of 12 Hyades dwarfs distributed tightly about the solar spectral type. This result was determined by comparison of the curve of growth over the linear and saturated branches of Fei lines with photospheric spectral models. A large surface fraction of active regions over these stars would probably result in an under-estimate of the abundance, but an identical estimate of \[Fe/H\] = $`+0.14\pm 0.01`$ from a smaller sample of Feii lines suggests that any bias due to chromospheric activity is small among the measured cases. While systematic errors may still be present, this result has been confirmed with a sample of F stars with $`v\mathrm{sin}i<30`$ km s-1 by Boesgaard (1989), finding \[Fe/H\] = $`+0.13\pm 0.03`$. Consequently, by comparison with the ASCA results, Fe abundance is higher in the photosphere compared to the corona of V471 Tau, within the uncertainties.

Cayrel et al. (1985) provide coronal abundances for Sii; \[Si/H\] = $`+0.16\pm 0.04`$. Varenne & Monier (1999) present photospheric O abundances for 26 F stars in the Hyades with an average value of \[O/H\] = $`+0.07`$. In these cases there are no significant differences between the coronal and photospheric populations, but the coronal abundances are poorly defined using the current data. However, as in other sources (e.g. Brinkman et al. 2001), the suggestion of a coronal FIP bias in V471 Tau hinges crucially on the measurement of Ne.

Critically, photospheric Ne has not bean measured in the Hyades because the low photospheric temperatures do not provide an environment for noble gas emission. We combine the Hyades abundances from Cayrel et al. and Varenne & Monier in Fig. 4. We do not include the two stars in the Cayrel et al. sample that are thought to be active sources. Data points from a sample of cluster members have been averaged and uncertainties obtained by calculating the standard deviations for each species. A linear least-squares fit yields \[m/H\] = $`0.069_{0.140}^{+0.140}+0.006_{0.017}^{+0.018}E`$, where $`E`$ is the FIP in eV. Uncertainties are the 90 percent confidence limits, i.e. the abundances are consistent with a flat distribution across the range of FIP.

We fit the slope of the X-ray FIP distributions from Fig. 3, including the Ne data. Clearly the FIP-trend is consistent with a constant if Ne is ignored. Least-square fitting yields slopes consistently of $`0.04\pm 0.03`$ for models 2, 3 and 4. Uncertainties are 90 percent confidence limits. Consequently, with Ne included, the coronal FIP distribution is steeper than the Hyades mean photospheric distribution.

## 5 Conclusion

We have measured model-dependent coronal abundances of Mg, Fe, Si, O and Ne in the corona of the post-common envelope binary V471 Tau. A single-temperature plasma model does not fit the data adequately. Binary evolution calculations have predicted that we would see no symptoms of element mixing during the common-envelope epoch and we indeed find no evidence for it; i.e., abundances are not unusual compared to the coronae of single stars or wide binaries. In all likelihood, the coronal Fe abundance is significantly less than the Hyades photospheric mean. This is in direct contrast to the ratio of Fe in the solar corona which is overabundant relative to the solar photosphere. There is evidence for an inverse FIP bias, although this relies entirely on our measurement of coronal Ne abundance. Due to the lack of low-ionization Ne lines in the optical band, the photospheric abundance of this element cannot be directly compared to the coronal value, adding further uncertainty. While CCD spectroscopy is not currently an ideal method for abundance determination, these data indicate that V471 Tau is a viable and interesting target for the grating instruments on-board the Chandra and XMM-Newton observatories for such a purpose.

This paper was based on data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center. We thank Nancy Brickhouse for supplying additional Fe and Ni atomic data.