text:
# Angular dimensions of planetary nebulae

## 1 Introduction

The angular dimensions belong to the most fundamental observational data for planetary nebulae (PNe). They are required if one wants to study important PN parameters such as linear dimensions, lifetimes, surface brightnesses. Several methods for determining distances or estimating nebular masses rely on observational measurements of the PN diameters. Also studies of the observed evolution of PNe and their central stars often use the angular PN diameters as one of the crucial observational data.

PN dimensions have been measured since the very begining of research interest in these objects and the list of papers with resultant values is very long. Consequently, for almost all the known PNe we have information about their angular dimensions, see e.g. the ESO-Strasbourg Catalogue (Acker et al. gpnc (1992)). However the data from many existing catalogues and databases should be used with caution. They have usually been compiled from numerous different sources based on different observational techniques and methods. Therefore the resulting data for different objects are not always comparable. For many PNe the diameters taken from different observational sources can differ by a factor of a few.

The problem is that different observational techniques and/or different methods for measuring dimensions give different results for many PNe. In the case of well resolved and fairly symmetric nebulae with a well defined outer rim the problem is simple and different methods give consistent results. Most of the PNe have, however, complex structures, soft outer boundaries, and low surface brightness emission extending to large distances. If one measures the maximum extension of the nebular emission in an image, as it has often been done while determining PN dimensions, than the result depends on how deep the image has been made. A reasonable solution in order to get results comparable for all objects is to measure nebular extension at a given level of surface brightness, e.g. at 10% of the maximum value.

Another problem arises in the case of small PNe when the nebular diameter is comparable to the angular resolution of the instrument. The observed image is then a convolution of an intrinsic nebular image and the instrumental profile. It has to be deconvolved in order to get a meaningful result. It is usually done using the so-called Gaussian deconvolution. The disadvantage of this method is that one has to adopt a shape of the intrinsic surface brightness distribution and the result depends on this assumption. The problem of determining dimensions of partially resolved nebulae has recently been discussed by van Hoof (vh (2000)). He has made a thorough model analysis of different methods for measuring partially resolved nebulae.

The question that arises, and which is important for using the measured dimensions in extensive, systematic studies of PNe, is to what extent the results from different methods are comparable. For instance, do the results from direct measurements at the 10% contour for extended nebulae mean the same as the values from a deconvolution method for compact objects? If the PNe were close to model nebulae studied by van Hoof (vh (2000)) the answer would be: yes, they measure more or less the same. However, real nebulae are often much different from the model ones and the answer is not a priori clear. We attempt to answer the above question in this paper.

In recent years, a number of observational surveys for studying PN structures and dimensions have been done. As surveys in radio wavelengths those of Aaquist & Kwok (ak (1990)) and Zijlstra et al. (zpb (1989)) should be mentioned. Over 400 PNe have been observed in these two surveys, mostly at 5 GHz. PN diameters have been determined either from a 10% contour level (for well resolved nebulae) or from a Gaussian deconvolution (for partially resolved objects). Recent optical surveys can be found in Schwarz et al. (scm (1992)), Manchado et al. (mgss (1996)), and Górny et al. (gsc (1999)) where about 600 PNe have been observed. However, the PN dimensions given in these three catalogues usually refer to the maximum extension of the recorded emission. As the ratio of the maximum to faintest recorded surface brightness can vary by a large factor from one image to another these diameters do not measure the same for all the objects in the catalogues.

We have measured angular dimensions for 312 PNe. One of the goals is to provide a homogeneous set of PN diameters, i.e. values which mean the same, as much as possible, for all the objects in the sample. We have applied direct measurements at 10% of the peak surface brightness as well as deconvolution methods. The results have been used to investigate reliability of the methods and for comparative studies of the results from different methods. This helped us to obtain a homogeneous set of the PN diameters for the whole sample without important systematic differences between compact and extended nebulae.

## 2 Observational material

The observational material used in this paper is, in principle, the same as in Schwarz et al. (scm (1992)) and Górny et al. (gsc (1999)). These are images of mostly southern PNe taken with the NTT and other telescopes at the ESO. For our measurements we have used only images in H$`\alpha `$ (or H$`\alpha `$ \+ \[NII\]). This has been done in order to have an uniform set of data and to avoid possible ionization stratification effects in the results which could have appeared if we had used, for instance, \[OIII\] images.

## 3 Methods of measurements

We have applied three methods of measuring PN dimensions. First one is a direct measurement of the nebular diameter at a level of 10% of the peak surface brightness, $`\mathrm{\Theta }_{10\%}`$, from the PN image. As the nebulae in majority are not spherically symmetric we have measured two diameters along perpendicular axes. One of the axes was directed along the longest extension of the nebular emission.

Second method is the Gaussian deconvolution. The full width at half-maximum (FWHM) values, $`\mathrm{\Phi }`$, have been obtained from a two-dimensional Gaussian fitting to the observed PN image. Before fitting the background has been subtracted from the image. The same applied to field stars allowed us to derive the FWHM of the resolution profile, $`\mathrm{\Phi }_\mathrm{b}`$. The deconvolved FWHM diameter, $`\mathrm{\Phi }_\mathrm{d}`$, can then be obtained from

$$\mathrm{\Phi }_\mathrm{d}=\sqrt{\mathrm{\Phi }^2\mathrm{\Phi }_\mathrm{b}^2}.$$
(1)
In order to get an estimate of the real nebular diameter, $`\mathrm{\Phi }_\mathrm{d}`$ has to be multiplied by a conversion factor which is a function of the observational resolution, $`\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}`$, and the adopted intrinsic surface brightness distribution of the nebula. Details of calculating the conversion factor for a constant surface brightness disc and constant emissivity shells can be found in van Hoof (vh (2000)).

Third method is the so-called second-moment deconvolution. This is an analogous method to the Gaussian deconvolution but the FWHM is determined from the second moment of the surface brightness distribution. The advantage of this method is that, unlike Gaussian deconvolution, the conversion factor is independent of the resolution of the observations. Details of applying second-moment deconvolution can be found in van Hoof (vh (2000)).

It is well known that the first method fails for small objects and therefore it is usually applied to well resolved nebulae. The deconvolution methods, on the contrary, are usually applied to small, partially resolved nebulae. In this paper we have, however, attempted to apply all the three methods to all the PNe in our sample. This has been done in order to be able to compare different methods and to study systematic differences between results from them.

## 4 Results of measurements

Table 1 presents the results of direct measurements done on our PN images. Following suggestions of van Hoof (vh (2000)) we give raw data from the deconvolution methods, i.e. the deconvolved FWHM, $`\mathrm{\Phi }_\mathrm{d}`$, and the image resolution, $`\mathrm{\Phi }_\mathrm{b}`$. Our final angular diameters of the PNe are given in Table 2. They are based on the measurements reported in Table 1 as well as on the discussion of the results from different methods presented in Sect. 5.

Column (1) in Table 1 gives the PNG number of the object while its usual name is in column (2). Column (3) gives the number of a figure showing the PN image in the original papers: ”s” and ”g” stand for Schwarz et al. (scm (1992)) and Górny et al. (gsc (1999)), respectively. The deconvolved FWHM, $`\mathrm{\Phi }_\mathrm{d}`$, and the resolution size, $`\mathrm{\Phi }_\mathrm{b}`$, from the Gaussian fit are in columns (4) and (5), respectively. The same but from the second-moment method can be found in columns (6) and (7). The results of direct measuremets at 10% level are given in column (8). Some PNe have relatively deep images and the nebular emission extends well beyond the 10% contour. In this cases we have also measured the dimensions of the $`3\sigma `$ emission level above the background. The results are in column (9). Column (10) gives the percentage of the $`3\sigma `$ level in the peak surface brightness. All the measurement results are in arcsec.

For some objects the deconvolution methods do not work well and Table 1 does not give results in these cases. It has usually happened when the nebula was extended with complex structures and/or had a strong central star or several field stars superimposed on the PN image. In these cases the numerical procedure fitting Gaussian or calculating moments either has failed or has given spurious results, or has just measured a star (central or one of the field stars). Note that the Gaussian fit fails more often that the second-moment method, particularly for larger nebulae.

Not all the objects published in Schwarz et al. (scm (1992)) and Górny et al. (gsc (1999)) appear in Table 1. Some PNe were too faint in the images or larger than the image frame so meaningful measurements have not been possible. On the other hand, images for a few objects (those with no entry in column (3) in Table 1) have not been published. For most of these objects their images (although obtained with different equipement) can be found in Manchado et al. (mgss (1996)).

For a few objects we had two images from independent observations. We give results from both images separately in two consecutive rows in Table 1 in these cases.

## 5 Comparison of the methods

In this section we compare the results from the three methods used to measure the PN dimensions. The aim is to study systematic effects between different methods. This will allow us to work out a procedure for deriving a consistent set of PN dimensions for all the objects, i.e. those well resolved as well as partially resolved, in our sample.

### 5.1 Deconvolution methods: Gaussian versus second-moment

As discussed above, both, Gaussian fitting and second-moment calculations have been applied to measure the FWHM of the resolution profile, $`\mathrm{\Phi }_\mathrm{b}`$, on our images. The results are given in columns (5) and (7) of Table 1. Both methods should give the same result if the instrumental profile is Gaussian. As can be seen from Table 1 this is essentially the case. The mean value and the standard deviation of the ratio $`\mathrm{\Phi }_\mathrm{b}(\mathrm{Gauss})/\mathrm{\Phi }_\mathrm{b}`$(s-m) calculated from the data in Table 1 is 0.965$`\pm `$0.094. Thus the basic assumption of the deconvolution methods, namely Gaussian resolution profile, is satisfied in our observational material.

The final results of the Gaussian and second-moment deconvolution methods are compared in Fig. 1. The figure plots the ratio of the angular diameters, $`\mathrm{\Theta }`$, derived from the two methods versus the ratio of the deconvolved FWHM, $`\mathrm{\Phi }_\mathrm{d}`$, to the resolution size, $`\mathrm{\Phi }_\mathrm{b}`$, (i.e. parameter $`\beta `$ in van Hoof vh (2000)) obtained from the Gaussian fit. The values plotted in Fig. 1 are geometric means of the values obtained from the two-dimensional measurements reported in Table 1. While deriving $`\mathrm{\Theta }`$ a surface brightness distribution from a constant emissivity sphere model (shell 0.0 in van Hoof’s notation) has been adopted for calculating the conversion factors. Note, however, that the diagram shown in Fig. 1 depends only weakly on the adopted surface brightness distribution.

As can be seen from Fig. 1, the two deconvolution methods do not give the same result for the bulk of objects. The second-moment method tends to give larger diameters than the Gaussian deconvolution. In the case of compact nebulae, i.e. when the Gaussian deconvolution gives $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.0}`$, the difference becomes particularly important and increases with decreasing PN size (as measured from the Gaussian deconvolution). Note that the vertical scatter of the points in Fig. 1, tends to become larger for more extended nebulae. On average $`\mathrm{log}\mathrm{\Theta }_{\mathrm{Gauss}}/\mathrm{\Theta }_{\mathrm{s}\mathrm{m}}`$ for objects with $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.0}`$ is $`0.096\pm 0.071`$ (mean and standard deviation). When adopting a shell-like model this value would go down and for an infinitely thin shell (limit 1.0 in van Hoof’s notation) it would become $`0.125\pm 0.070`$. If a constant surface brightness disc is adopted the figures would be $`0.108\pm 0.071`$.

### 5.2 Gaussian deconvolution versus direct measurements at the 10% level

Fig. 2 compares the results from the Gaussian deconvolution with the 10% contour measurements. Similarly as Fig. 1, it shows the ratio of the PN diameters versus the Gaussian $`\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}`$ ratio. Again, the plotted values are geometric means from the measurements in Table 1. While obtaining the diameters from the Gaussian method the constant emissivity sphere model has been adopted. Note that the diagram in Fig. 2 depends on the adopted surface brightness distribution more significantly than Fig. 1. The general pattern of the diagram remains the same but the points can shift vertically.

Fig. 2 shows an increasing discrepancy between the 10% diameter and the diameter from the Gaussian deconvolution with the decreasing value of $`\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}`$. This is an expected effect and is due to an overestimate of the PN diameter from the direct measurements for partially resolved PNe. Note, however, that a tight correlation seen in Fig. 2 for $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.0}`$ is partly spurious. For $`\mathrm{\Phi }`$ approaching $`\mathrm{\Phi }_\mathrm{b}`$ the y-axis in Fig. 2 tends to plot a reversal of the value plotted in the x-axis.

For small nebulae the diameters from the 10% contour would be more accurate if determined from deconvolved images. We have not attempted to apply this rather cumbersome procedure to our images. Instead we have found that a significant improvement of the 10% diameter can be obtained using a formula analogous to the standard FWHM Gaussian deconvolution (Eq. (1)) but relevant to a 10% level of the maximum, i.e.

$$\mathrm{\Theta }_{10\%,\mathrm{d}}=\sqrt{\mathrm{\Theta }_{10\%}^2(1.823\mathrm{\Phi }_\mathrm{b})^2}.$$
(2)
We will call $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ the deconvolved 10% diameter.

Fig. 3 shows the same as Fig. 2 but using $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ instead of $`\mathrm{\Theta }_{10\%}`$. Unlike in Fig. 2, the distribution of the points in Fig. 3 does not show any significant trend with $`\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}`$. This demonstrates that Eq. (2) allows to significantly improve the reliability of the PN dimension estimate from the 10% contour measurements for small, partially resolved nebulae. $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ seems to be an estimate of compact PN diameters as reliable as that from the Gaussian deconvolution.

However, as can be seen from Fig. 3, the Gaussian deconvolution tends to underestimate the PN dimensions compared to the 10% deconvolved diameters. Note that the value of this systematic difference depends on the surface brightness distribution adopted in the Gaussian method. It is smallest for the constant emissivity sphere, i.e. for the case presented in Fig. 3, and largest for the infinitely thin shell model.

Fig. 3 suggests that the two methods give the most consistent results for rather compact but not too small nebulae, i.e. having $`0.2\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$. The mean value of $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Theta }_{\mathrm{Gauss}}`$ for these objects in Fig. 3 is 0.028$`\pm `$0.041. For an infinitely thin shell (limit 1.0 in van Hoof’s notation) and a constant surface brightness disc the result would be 0.159$`\pm `$0.043 and 0.085$`\pm `$0.042, respectively.

A greater discrepancy of the results for $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}>\mathrm{\hspace{0.17em}0.3}`$ can be seen from Fig. 3. The mean values of $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Theta }_{\mathrm{Gauss}}`$ derived for these objects are: 0.052$`\pm `$0.062 (sphere), 0.195$`\pm `$0.062 (thin shell), and 0.115$`\pm `$0.062 (disc). These nebulae show resolved structures on the images and their surface brightness distributions are usually far from a Gaussian. Therefore a Gaussian fitting to such objects gives uncertain, sometimes even spurious, results.

Because of low statistics, no meaningful analysis can be done for the region $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}<0.2`$. A large scatter of the points is certainly due to increasing sensitivity of the results to the observational uncertainties (quadratic subtraction of two nearly equal values).

### 5.3 Second-moment deconvolution versus direct measurements at 10% level

Figs. 4 and 5 compare the results from the 10% level measurements and the second-moment deconvolution. Fig. 4 is analogous to Fig. 2 and uses ”raw” measurements from the 10% contour. A tendency of overestimating diameters for compact PNe from the 10% measurements can be seen in Fig. 4. It is, however, much less conspicuous than in Fig. 2. The reason is that, as discussed in Sect. 5.1, for partially resolved PNe the second-moment deconvolution results in larger diameters than the Gaussian deconvolution.

The deconvolved 10% diameter, using Eq. (2) and $`\mathrm{\Phi }_\mathrm{b}`$ from the second-moment calculations, have been applied in the diagram shown in Fig. 5. As the deconvolution of the 10% diameter most affects the results for compact nebulae the positions of these objects are significantly shifted down in Fig. 5 comparing to Fig. 4.

Fig. 5 shows that the second-moment deconvolution usually overestimates the PN size compared to the 10% contour measurements. For rather large nebulae, i.e. having $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$ in Fig. 5, $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Theta }_{\mathrm{s}\mathrm{m}}`$ is $`0.033\pm 0.061`$. The discrepancy increases for less resolved nebulae and for $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}<\mathrm{\hspace{0.17em}0.3}`$, $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Theta }_{\mathrm{s}\mathrm{m}}`$ is $`0.132\pm 0.115`$.

As the conversion factors in the second-moment deconvolution are independent of $`\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}`$ adopting the infinitely thin shell model would shift all the points in Fig. 5 upwards by factor 0.111. For the uniform disc model the shift would be of factor 0.048. By the same factors would respectively increase the mean values of $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Theta }_{\mathrm{s}\mathrm{m}}`$ (obviously, the standard deviations would remain the same).

### 5.4 Testing the methods on degraded PN images

One of the main problems in the present study is to find which method(s) gives reliable diameter estimates for partially resolved nebulae. The analysis discussed in the previous subsections shows that the deconvolved 10% diameter and the Gaussian deconvolution give relatively consistent results down to rather poorly resolved objects, although the Gaussian deconvolution systematically underestimates the PN diameter compared to the measurements of the 10% contour. The value of this underestimate depends on the model nebula adopted in the Gaussian deconvolution. On the other hand, the second-moment deconvolution seems to become unstable for compact objects as suggested by an increasing scatter of points in Fig. 5 for less and less resolved nebulae. However, we cannot be sure that this interpretation is correct. Nor we can say which method and when breaks to give reliable results. The reason is that we can compare only the results from different methods but we do not know what real diameters of poorly resolved nebulae are.

We have attempted to clarify the above problems by applying the methods of diameter measurements to some of our images artificially degraded to worse and worse seeing conditions, i.e. convolved with a larger and larger Gaussian profile. From the original images we know the true PN sizes so we can compare them to the results from different methods applied to less and less resolved PN images.

Although the idea is simple it could have been applied to a limited number of our images only. Obviously we have to consider rather well resolved nebulae but not too much. The object has to be several times smaller than the image size. Otherwise the image frame can impose an important cutoff when the image is convolved with a Gaussian having FWHM comparable to the PN size. The nebula has to be well exposed so its image keeps a good contrast while beeing spread over larger and larger surface. There cannot be too many bright field stars, otherwise they merge with the PN image in the convolution process. It is also desirable that among the objects investigated different morphologies are represented.

We present the results of the procedure applied to images of 5 PN, i.e. (in paranthesis usual names and morphological descriptions are given) 021.8–00.4 (M 3-28; bipolar, butterfly-like), 174.2–14.6 (H 3-29; ring-like but strongly asymmetric with one strong blob), 285.7–14.9 (IC 2448; round with a halo), 286.3–04.8 (NGC 3211; ring-like with two blobs), and 345.2–08.8 (Tc 1; round with an irregular halo). The selected objects have diameters between 10 and 20 arcsec. The original images have been obtained with a seeing varying from 0.5 to 2.5 arcsec. We have used our deconvolved 10% diameters from original images as true PN diameters for comparison with the results from degraded images (note that the deconvolution has introduced only minor corrections to the 10% contour measurements in these cases). The original images have been convolved with a Gaussian having FWHM increasing from 2.5 to 20.0 arcsec. Then our three methods have been applied to the resulting images.

We note that the results from measurements on the degraded images reproduce very well the general trends seen in Figs. 15.

In Fig. 6 the results of the deconvolved 10% diameter method applied to the degraded images are compared to the deconvolved 10% diameter from the original (non-degraded) images ($`\mathrm{\Theta }_{\mathrm{true}}`$). From this figure we can conlude that the deconvolved 10% diameter measures very well the PN size down to $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$. This limit is practically equivalent to $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}(\mathrm{Gauss})/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.05}`$. Below this value the method becomes unstable.

Results of the Gaussian deconvolution applied to the degraded images are presented in Fig. 7. A constant emissivity sphere model has been adopted. We see that for well resolved nebulae the Gaussian deconvolution underestimates the PN diameter. This is consistent with the trend seen in Fig. 3. For $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$ the underestimate is smaller, again similar as in Fig. 3. For all the points with $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}<\mathrm{\hspace{0.17em}0.3}`$ in Fig. 7 the mean value and the standard deviation of $`\mathrm{log}\mathrm{\Theta }_{\mathrm{Gauss}}/\mathrm{\Theta }_{\mathrm{true}}`$ are $`0.021\pm 0.032`$. This result is very close to the value $`0.028\pm 0.041`$ of $`\mathrm{log}\mathrm{\Theta }_{\mathrm{Gauss}}/\mathrm{\Theta }_{10\%,\mathrm{d}}`$ found in Sect. 5.2 for objects having $`0.2\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$ in Fig. 3. For $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}0.2`$ in Fig. 7 the scatter increases which suggests that the Gaussian deconvolution becomes then unstable. Comparison of Figs. 7 and 6 seems to indicate that the Gaussian deconvolution gives somewhat more stable results than the deconvolved 10% diameter for objects having $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$ or, equivalently, $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}(\mathrm{Gauss})/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.05}`$. A better statistics would be required to clarify this point.

The results of the test of the second-moment deconvolution are presented in Fig. 8. A constant emissivity sphere model has been adopted. The figure shows that the second-moment deconvolution gives uncertain estimates of the PN diameters, especially for partially resolved nebulae ($`\mathrm{log}\mathrm{\Phi }_\mathrm{d}/\mathrm{\Phi }_\mathrm{b}<\mathrm{\hspace{0.17em}0.3}`$). Adopting a shell or disc model in the deconvolution would decrease the systematic difference between $`\mathrm{\Theta }_{\mathrm{s}\mathrm{m}}`$ and $`\mathrm{\Theta }_{\mathrm{true}}`$ in Fig. 8 but the uncertainity (scatter of the points) would remain the same. Our experience suggests that this method is sensitive to faint outer structures (halo, outer butterfly-like filaments) particularly in partially resolved nebulae.

## 6 PN diameters

One of the primary goals of our work is to obtain a homogeneous set of PN dimensions, a set of results which have, as far as possible, the same physical meaning for all the objects in our sample. As a primary method for deriving the PN dimensions we have adopted the direct measurements at the 10% level of the peak surface brightness. This seems to be the only reasonable method for well resolved nebulae. In order to account for the finite resolution in the case of smaller objects we have used the deconvolved 10% diameters, $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$, instead of $`\mathrm{\Theta }_{10\%}`$. As discussed in Sect. 5.4 $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ is reliable down to $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$. For the smallest objects when results from any method become more and more uncertain we have used $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ and the results from the Gaussian deconvolution. Following discussion in Sect. 5.4 we have prefered not to use the results from the second-moment deconvolution. In more detail our prescription is as follows.

For all the objects with $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$ we have adopted $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$. We have used $`\mathrm{\Phi }_\mathrm{b}`$ from Gaussian fit. For more compact objects, i.e. when $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}<\mathrm{\hspace{0.17em}0.3}`$ we have adopted a mean value from $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ and $`\mathrm{\Theta }`$ resulting from the Gaussian deconvolution adopting the constant emissivity sphere model. In order to correct for the systematic underestimate of the Gaussian method, as discussed in Sects. 5.2 and 5.4, the Gaussian $`\mathrm{\Theta }`$ have been multiplied by factor 1.06 before deriving the mean value.

The results are presented in Table 2. Columns (1) and (2) give the PNG numbers and the usual names, respectively. The diameters derived according to the above prescription are given in column (3). In the case of partially resolved objects having $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}<\mathrm{\hspace{0.17em}0.3}`$ only one value is given which is the geometric mean from the two dimension measurements. A colon indicates compact objects for which the Gaussian deconvolution resulted in a diameter different by more than 30% from the deconvolved 10% diameter. Objects practically unresolved in our images are marked with ”stellar”. For a few objects having two independent measurements in Table 1 mean values derived from them are given in Table 2.

We have compared the results from Table 2 with the diameters from our data base resulting from a compilation of the data published in the literature. This data base has originally been made for purposes of the Strasbourg-ESO catalogue (Acker et al. gpnc (1992)). Its updated version has been described and used in Siódmiak & Tylenda (st (2001)). The results of the comparison are shown in Fig. 9. In most cases the agreement is quite good although our measurements tend to give somewhat larger diameters than the values adopted from the literature. The mean value of the ratio of our diameters (from Table 2) to the compiled diameters is 1.24 with a standard deviation of 0.71. In some cases, however, the discrepancy is significant. For small nebulae a large scatter of the points in Fig. 9 is not surprising. Any measurement, VLA or optical, of a partially resolved nebulae gives an uncertain result. Significant discrepancies in the case of larger nebulae usually concern bipolar nebulae having a compact core, e.g. IC 4634 (000.3+12.2), M 3-28 (021.8–00.4), M 2-48 (062.4–00.2), He 2-104 (315.4+09.4), Hb 5 (359.3–00.9). The published diameters from VLA or optical measurements apparently referred to the core whereas the 10% contour from our images encompassed a significant portion of the outer bipolar structures.

In Fig. 10 we have compared the diameters from Table 2 with the results of VLA surveys of Aaquist & Kwok (ak (1990)) and Zijlstra et al. (zpb (1989)). Aaquist & Kwok observed mainly compact PNe. They have applied three methods for deriving the PN diameters but it is not clear which method(s) has been used to measure a particular object. As can be seen from Fig. 10, our diameters are systematically larger than those of Aaquist & Kwok, on average by factor 1.7. The agreement of our results with those of Zijlstra et al. is better especially if nebulae greater than $``$10 arcsec are considered and the diameter of Zijlstra et al. resulted from their 10% contour (full circles in the figure). In four cases the discrepancy is larger than factor 2. All the cases concern bipolar nebulae. For three of them, i.e already mentioned M 3-28, M 2-48, as well as NGC 2440 (234.8+02.4), our diameters are larger and the reason has been explained above. In the case of M 2-9 our diameter is significantly smaller and it seems that Zijlstra et al. measured the angular extend of the elongated, bipolar lobes whereas our value is a geometric mean from the major and minor axis measurements.

## 7 Discussion and conclusions

We have measured angular dimensions for 312 PNe. As a primary definition of the PN size we have adopted the extension of nebular emission up to the 10% level of the maximum surface brightness. This seems to be a reasonable choice for extended, well resolved nebulae. For partly resolved objects deconvolution methods have to be used. So far the Gaussian deconvolution has primarily been applied.

We have, however, found that for compact nebulae the 10% contour measurements can also give satisfactory results provided that they are corrected for the finite resolution using Eq. (2). As has been shown in Sect. 5.2 and 5.4 this simple deconvolution of the 10% contour measurments makes the resulting diameters reliable down to nebulae as compact as $`\mathrm{log}\mathrm{\Theta }_{10\%,\mathrm{d}}/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.3}`$ (or equivalently $`\mathrm{log}\mathrm{\Phi }_\mathrm{d}(\mathrm{Gauss})/\mathrm{\Phi }_\mathrm{b}\mathrm{\hspace{0.17em}0.05}`$). The derivation of $`\mathrm{\Theta }_{10\%,\mathrm{d}}`$ has several advantages over the deconvolution methods. It is simpler to apply. No assumption has to be made about the surface brightness distribution. It is a more universal method as it can be applied to compact and partially resolved PNe as well as to extended and well resolved nebulae where the deconvolution methods usually fail in practice.

The Gaussian deconvolution gives fairly reliable and stable results for partly resolved nebulae. However it systematically underestimates the diameters compared to the 10% contour measurements. The difference depends on the model nebula adopted in the Gaussian deconvolution. It is only 5-7% for the constant emissivity sphere but increases up to 45% for a thin shell. This result means that in catalogues and data bases the dimensions of compact PNe, usually derived from the Gaussian deconvolution, are probably underestimated compared to those of extended nebulae. It is difficult to estimate this systematic difference. The deconvolution procedures used in observational surveys, like in Aaquist & Kwok (ak (1990)) or Zijlstra et al. (zpb (1989)), were significantly simplified compared to that of van Hoof (vh (2000)). Usually the same value of the conversion factor has been used for all the objects. This certainly introduced additional uncertainties in the final results.

The second-moment deconvolution appeared to give results often uncertain and in poor agreement with the results from the other two methods, in particular for compact nebulae.

Our study does not solve definitevely the problem of measuring the PN dimensions. It makes however an important step in the field in the sense that it gives a quantity which means physically the same and is comparable for all the objects in our sample. One has, however, to remember that the (faint) nebular emission can extend well beyond the diameters given in Table 2. This can be seen for many objects in our sample by comparing column (9) and column (8) in Table 1.

###### Acknowledgements.

We are very grateful to the referee (J.J. Condon) whose comments resulted in a significant improvement of the work reported in this paper. This work has partly been supported from a grant No. 2.P03D.020.17 financed by the Polish State Committee for Scientific Research.


format: text

text:
# Error Diagrams and Temporal Correlations in a Fracture Model with Characteristic and Power-Law Distributed Avalanches

## 1 Introduction

The fracture of heterogeneous materials has been the subject of intensive research since many years and is one of the oldest concerns of science hans90 . This is not only due to the evident practical and potential profits, but also because the understanding of fracture processes at a basic level has shed light on other phenomena like earthquake occurrence turcotte ; sornette00 . Relatively simple models such as random resistor networks arcan85 , beam networks roux85 , spring-block models burridge ; ofc1 ; ofc2  and the so-called Fibre Bundle Models (FBMs) hans90 ; sornette00 ; daniels ; newman95 ; hemmer ; moreno00  have guided our way to more complex models of fracture. To date, we actually know many of the failure properties and its dynamics. Yet, the range of phenomena associated with fracture has not been cast into a definite physical and theoretical treatment. Some years ago, there was a burst of activity in what today we know as self-organized criticality bak1 ; bak2 , a theoretical framework that was soon applied to the study of fracture processes and avalanche like phenomena in disordered system. Soon afterwards, some interesting claimings were put out as, for example, the suggestion that the earth’s crust is in a self-organized critical state, being the Gutenberg-Richter (GR) law gr  for the frequency of earthquakes with a given magnitude a result of this self-organization process bak3 ; sornette89 ; ito .

From a more practical point of view, perhaps the most important concern related to fracture processes that remains unsolved is their predictability sornette00 ; kb02 . One is interested in knowing not only under what circumstances a material will fail or how the distribution of energy releases looks like, but also when a material will break down or an event of a given magnitude (usually we are interested in largest quakes) will take place. This problem is far from being solved and different approaches have been adopted in order to increase the predictive power of actual methods kb02 . For instance, in seismicity, there is an unavoidable degree of arbitrariness in classifying earthquakes. When a relative big earthquake occurs, it can be either a foreshock of a larger subsequent event, an aftershock of a preceding large quake or the mainshock itself. Thus there is no way to differentiate these events in real time and many algorithms and methods are useless from a practical point of view as they perform only a posteriori. However, regardless of how a big event may be classified, the forecasting of large quakes is the main goal to be achieved. This is the type of challenge that we face here.

In this paper, we analyze a dissipative fibre bundle model of fracture with only one parameter that characterizes the heterogeneity degree of the system moreno99 . Noticeably, the model exhibits two very distinct failure regimes obeying different avalanche size distributions. This allows the study of the occurrence of large avalanches, for the two phases, within the same model. In one of these regimes, the avalanche size distribution shows a power law regime for small events separated by a gap from the larger events, which are of the order of the system size. This regime is reminiscent of the characteristic earthquake behavior found in some seismic faults wes94 ; sieh96 . On the other hand, in the second regime, the system self-organizes into a critical state characterized by a unique power law avalanche distribution similar to the GR law. Following the Molchan method molchan  we compute the error diagrams that quantify the forecasting of large cascade events and analyze their correlations in time. Our results could be potentially applied to the study of earthquake occurrence.

## 2 Failure Model

The basic ingredients and main features of the model used here are as follow moreno99 :

1. A set of $`N`$ elements is located on the sites of a supporting square lattice of linear length $`L`$. Each element represents a fibre, or in general terms, a small volume fraction of the material under study.
2. To each element $`i`$, one assigns a random threshold strength $`\sigma _{i_{th}}`$ taken from a probability distribution. The threshold values represent the maximum load each of the elements can support before failure occurs. We use henceforth a Weibull distribution $`P(\sigma _{th})=1e^{(\sigma _{th}/\sigma _0)^\rho }`$ to assign the failure thresholds assuming a reference load $`\sigma _0=1`$. $`\rho `$ is called the Weibull index; the bigger it is, the narrower the range of threshold values is moreno99 .
3. At each time step, the load on the system is increased by calculating the minimum load required for one element to break and adding this amount to all the elements, i.e., $`\sigma _i\sigma _i+\nu `$, where $`\nu =min\{\sigma _{i_{th}}\sigma _i\}`$.
4. Once the load acting on an element surpasses its failure threshold, it fails and the load it bears is equally redistributed among all the elements of the set irrespective of their state. This amounts to a global fibre bundle version where long range interactions are assumed note1 .
5. The fraction of load that would correspond to already failed elements is simply dissipated so that the model is nonconservative. Moreover, the degree of dissipation depends on both the avalanche size and its internal structure moreno99 .
6. The redistribution of load to non-broken fibres may trigger more failure events. In such a case, step $`4`$ is repeated until the system reaches a state where all the elements support a load below their respective threshold values.
7. When an avalanche has come to an end, the external driving, step $`3`$, is repeated and the elements that have been broken in the preceding avalanche are healed assigning to them new random strengths and a load equal to zero.

The model is thus constructed by implicitly taking into account some general properties of any fracture process. Two ingredients of our model are worth mentioning. First, we assume dissipation through broken elements. This amounts to consider that the energy dissipated during a given fracture event is proportional to the size and the duration of that event. Besides, we incorporate healing only after an avalanche has ended because when the ongoing avalanche is developing, the elements have no time for healing and broken regions should not accumulate stress in that time interval. We additionally note that the results presented here are valid, strictly speaking, in the limiting case of very slow driving.

There are two different time scales in the model moreno99 . One refers to the external driving of the system and defines the natural time interval between failure events. We measure the time elapsed until a given event $`k`$ as the sum of the incremental amounts of load added at every external loading of the system prior to the event $`k`$. This definition of time assumes that the rate at which the system is loaded can be considered constant, and hence that $`\nu `$ is linearly proportional to the time elapsed since the last event. This approximation holds, for example, in seismology. The second time scale of the process quantifies the avalanche lifetimes, and it is measured as lattice updates. That is, we consider the avalanches to be instantaneous with respect to the external driving and separated by an interval of time whose length is proportional to the amount of energy added to the system during the last external input.

The iteration of the rules sketched above leads the system to two different behaviors depending on the level of heterogeneity assumed when assigning the threshold values moreno99 . Figure 1 shows the distributions of avalanche sizes, $`P(s)`$, for the two failure regimes of the model. The upper panel shows the pattern typical of characteristic quakes which are distributed quasi-periodically in time. This regime, where there are large fluctuations of the load on the system, is obtained when the bundle of fibres is homogeneous (large values of the Weibull index $`\rho `$), that is, when the deviation of the thresholds from a mean value is not large. Note that as the system size is increased, large events are better resolved and there are no intermediate avalanches.

A richer and more complex behavior is found when the system departs from homogeneity and the bundle of fibres gets heterogeneous (small values of $`\rho `$). In this case, irrespective of the system size, a power law distribution for the avalanche sizes is obtained. This is depicted in the lower panel of Figure 1 for $`\rho =2`$ and $`N=50000`$ fibres. As we shall see in the next section, the quasi-periodic nature of the first regime makes the forecasting of large avalanches feasible and eventually, for large system sizes (strictly speaking, in the thermodynamic limit), one can always forecast such events with no errors. This is not anymore the case when the system can be regarded as heterogeneous, as non-trivial temporal correlations and clustering of large events come into play making the forecasting of catastrophic avalanches a hard task.

## 3 Error Diagrams and Temporal Correlations

In order to inspect whether or not we would be able to predict large events, we first compute the errors diagrams for both failure regimes. They are a quantitative measure of the success in forecasting the occurrence of a given event and were introduced in seismology several years ago by Molchan molchan  with the aim of rigorously evaluating different algorithms for large earthquakes prediction, and subsequently used by other authors kb02 ; newman02 ; sornette02 ; vp03 .

An error diagram requires the identification of a target to be forecasted using alarms. Then, one plots the trade-off between the fraction of target events that were not predicted, $`f_e`$, and the fraction of time that the alarm was active, $`f_a`$. Ideally, one would like to get the smallest number of failures to predict at the minimum duration of the alarms. Points close to the origin of an error diagram represent the best predictions. Hence, the efficiency of a given forecasting algorithm is measured by how fast $`f_e`$ decreases when increasing $`f_a`$. In practice, one switches on the alarm at a given time interval $`\tau `$ after the preceding event. If $`\tau =0`$, the fraction of failures to predict is zero, but the fraction of alarm time is 1. On the contrary, if $`\tau `$ is equal to (or greater than) the time interval between two successive events, $`f_a=0`$ but $`f_e=1`$. In this way, the result of a completely random (Poissonian) process is a diagonal line.

Let us now construct the error diagrams for each of the regimes of the fracture model under study. First, we define the targets as the large avalanches whose sizes exceed a given threshold $`s_c`$ and fix the time interval $`\tau `$ at which the alarms are switched on after every large failure event. By varying the length of $`\tau `$ and computing the fractions of failures to forecast and of alarm times, $`f_e(\tau )`$ and $`f_a(\tau )`$ respectively, the error diagrams shown in Fig. 2 are obtained, where the dependency on $`\tau `$ has been eliminated. The figure clearly illustrates how different with respect to successful forecast the two failure regimes are. Open circles represent the results of the algorithm when applied to the prediction of characteristic events in the parameter region when the system is homogeneous. Here, target events occur quasi-periodically and thus setting $`\tau `$ close to the mean time interval between large avalanches ensures a successful prediction with a minimum of alarm time. For example, one can choose $`\tau `$ so that $`f_a0.3`$ and $`f_e0.2`$, a rather good forecasting. Besides, increasing the system size improves the prediction of characteristic events since there appears a gap between small and large events making the definition of targets less noisy. In this case, the error diagram approaches a delta function at the origin as large failure events occur almost periodically in time.

On the contrary, when the system gets heterogeneous, in the error diagram $`f_a`$ is always above $`f_e`$ (full circles). Here, we consider an avalanche as a target when its size falls on the edge of the distribution $`P(s)`$, where finite size effects appear. Reducing the target thresholds move the error diagram close to the diagonal line. In fact, the occurrence of an event of any size is a random process as can be observed in the figure when the threshold is set equal to 1.

Another interesting information that can be extracted from the error diagrams for the second regime is the existence of some degree of non-trivial temporal correlations. The departure from a random process when the target threshold is increased indicates that large events show clustering in time. This is clearly appreciated in Fig. 3, where two sequences of quakes for a system made up of $`N=50000`$ elements and $`\rho =2`$ has been represented. The upper series is the occurrence of events with sizes equal to or greater than $`10^2`$ while the lower sequence is obtained when $`s_c`$ is equal to $`10^3`$ (upper limit for events not affected by finite-size effects).

We also present in Fig. 4 the distribution of inter-occurrence time $`\mathrm{\Delta }t`$ for two values of the target thresholds. We see that when $`s_c`$ is increased so that only large avalanches are taken into account, the distribution $`D(\mathrm{\Delta }t)`$ radically changes from an exponential decay (typical of a random process) to a decaying power law. This is in agreement with the error diagrams previously shown and provides further evidence that small avalanches are not correlated while large ones are. Without providing a figure, we point out that the distribution of inter-occurrence time for the first failure regime is a well-behaved function with a clear mean value and a dispersion that decreases as the system gets larger.

## 4 Discussion and Conclusions

The results shown in the preceding section indicate that the clustering of large avalanches makes the prediction of such events a hard task. The shape of the error diagrams and the high rate of prediction failures can be understood by noticing that the average inter-occurrence time between all large events is not a useful reference in this case. Two other quantities seem to be more fundamental: the inter-cluster average time defined as the time interval between clusters of large avalanches and the mean intra-cluster time defined as the time between large events within the same cluster. Since the former is much larger that the latter, the fraction of failures to predict is raised when $`\tau `$ is large. On the other hand, setting $`\tau `$ too small to match the intra-cluster time produces a significant increment of $`f_a`$. Thus both phenomena contribute to the error diagrams in such a way that the forecasting is worse than a random prediction.

The clustering of large events is due to the way in which energy is dissipated and not to the power-law nature of the avalanche size distribution. The average inter-occurrence time between large quakes can be considered as a measure of the buildup time for correlations in the system. When a series of big cascades occurs, the system ends up in a state where most of the elements are unloaded and a significant input of energy (or a large time interval) is needed in order to start another series. In other words, the way in which energy is dissipated creates stress correlations giving rise on its turn to the observed clustering. This does not hold anymore when avalanches of smaller sizes are considered. In this case, the occurrence of a small avalanche does not significantly alter the state of the system. A small event can not influence the occurrence of any event at short times.

The same kind of temporal correlations is seen in other models displaying self-organized criticality. This is the case of the Olami, Feder and Christensen (OFC) spring-block model for earthquakes ofc1 ; ofc2 . They found some years ago that non-trivial temporal correlations of the sort obtained here show up when the model is made nonconservative. By calculating the coefficient of variation, it was shown that large earthquakes show clustering while small ones are uncorrelated in time. Besides, when the model is conservative, neither spatial nor temporal correlations were found ofc2 . The fact that dissipation appears to be a basic ingredient for a model to show correlations in time and/or space is corroborated in Fig. 5, where the error diagrams for the original Bak, Tang and Wiesenfel (BTW) model bak1 ; bak2  are depicted. As can be clearly seen in the figure, the model can not be distinguished from a random process with regard to its predictability. Changing the target thresholds does not improve the forecasting. Thus, in this model, small avalanches and large ones are completely equivalent and all correlations are ruled out.

We would like to remark that the results obtained in this paper may be applied to the study of earthquake occurrence. One of the few well-known facts in seismology is the clustering of earthquakes in time and space. While small earthquakes seem to be uncorrelated joh85 , large earthquakes display strong clustering kagan ; sergio . On the other hand, by tuning a single parameter in our model, one moves from a characteristic earthquake scenario vp02  to a regime where the avalanche sizes are power law distributed (GR law). The possibility of exploring such different regimes within the same model is another motivation for further research in this direction note2 .

Summing up, we have analyzed a fracture model with non-trivial temporal correlations. The results presented seem to support the idea that models aimed at modeling real earthquakes should incorporate some kind of dissipation and self-organization. On the other hand, the clustering of large quakes does not help to improve the predictability of devastating events since two relevant time scales come into play. It might be interesting to compute the error diagrams with a tunable alarm time in such a way that one can better resolve series of large events and periods of stasis.

###### Acknowledgements.

Y. M. acknowledges financial support from the Secretaría de Estado de Educación y Universidades (Spain, SB2000-0357). M. V-P. is supported by the Ph.D grant B037/2001 of the DGA. This work has been partially supported by the Spanish DGICYT project BFM2002-01798.


text:
# Young isolated neutron stars from the Gould Belt

## 1 Introduction

Observations of isolated neutron stars (INSs) are important for gaining deeper insight on their structure and thermal evolution and ultimately might prove decisive in unveiling the physical properties of matter at ultra-high densities. Up to a decade ago, the only known INSs were active radiopulsars, with the addition of the $`\gamma `$-ray pulsar Geminga. Despite the fact that X-ray emission from some radio-pulsars was already detected by Einstein, it was in the ’90s that ROSAT, thanks to its sensitivity in the 0.1-2 keV band, gave a clearer picture of the faint X-ray emission produced by the cooling surface of the closest INSs. ROSAT, supplemented by more recent observations by Chandra and XMM-Newton, revealed a variety of types of behavior in the X-ray emission from INSs and its relation with radio activity.

In particular, a substantial contribution of ROSAT to this field has been the discovery of a group of seven radio-quiet, close-by, thermally-emitting INSs, dubbed ROSAT INSs (RINSs) and sometimes referred to as the “magnificent seven” (see for recent references Treves et al. 2000, Zampieri et al. 2001, Haberl & Zavlin 2002; see also Popov & Prokhorov 2003). The nature of these sources has been controversial and two interpretations were proposed, either as conventional middle-aged cooling NSs or as very old NSs accreting the interstellar gas (see again Treves et al. 2000). Although no compelling evidence has been brought in favor of either picture as yet, it is now generally believed that at least two of these objects (the brightest ones, RX J1856-3754 and RX J0720.4-3125, e.g. Braje & Romani 2002) are indeed relatively young cooling NSs.

Since the seven RINSs have remarkably similar observed properties, it is quite natural to assume that they belong to the same class, i.e. all of them are close-by, cooling NSs. This, however, poses a major problem. A useful and standard way to study a population of sources is to compute the Log N – Log S distribution starting from a model, and then compare it with observations. This has been done for RINSs by Neuhäuser & Trümper (1999) and Popov et al. (2000b). The main conclusion of the latter investigation is that the typical spatial density of radio pulsars in our Galaxy is too low to explain RINSs. In other words, the assumption that RINSs and ordinary radio-pulsars derive from the same parent population underpredicts the number of observed “coolers”, i.e. NSs which are hot and close enough for their thermal emission to be detected in X-rays.

An obvious solution is to invoke a local (both in space and time) overabundance of NSs with respect to those originating in the Galactic disk and seen now as radiopulsars. The main goal of this paper is to investigate the possible origin of these objects. Here we suggest that the likely birthplace for many of the INSs in the Solar proximity is the Gould Belt (see Popov et al. 2002 for a preliminary discussion). The Gould Belt is a compound of young stellar associations extending in a ring tilted at about $`20^{}`$ from the Galactic plane, roughly centered at the position of the Sun and with a radius of $`1`$ kpc. Originally discovered in the middle of the 19th century, the Belt was studied in detail by Benjamin Gould (see Stothers & Frogel 1974 for a historical overview). Most of the known close-by young star associations (in Ophiuchus, Orion, Perseus, Scorpius etc.) belong to the Belt. Since the Gould Belt is a relatively young stellar system it contains a comparatively large fraction of massive stars. As the Solar proximity is embedded within the Belt, it is likely to have harbored quite a number of massive young stars which are the progenitors of NSs.

In Sect. 2 we perform a population synthesis of close-by young INSs, assuming that they are born both in the Gould Belt and in the Galactic disk with an assigned mass spectrum. The dynamical evolution of the young NSs population in the Galactic potential is then followed together with the thermal evolution. This allows us to compute the present surface temperature and spatial distributions of young, close-by NSs and finally their X-ray count rate, once a model for interstellar absorption has been prescribed. This is used to compute the Log N – Log S distribution which is then compared with observations in Sect. 3. We refer to ROSAT PSPC in evaluating the X-ray count rate and in confronting our model with data. Even if Chandra and XMM-Newton provide much better observations of individual sources and small fields, ROSAT all-sky data are still the most complete by far in soft X-rays. Discussion and conclusions follow in Sect. 4.

## 2 The Model

In this section we describe the method we used to compute the local distribution of young, cooling neutron stars. The main ingredients of our model are: the spatial distribution of NS progenitors, the NS formation rate, the NS cooling history, and the interstellar medium (ISM) distribution (the latter determines the interstellar absorption and hence affects the X-ray count rate). Even if we are concerned with a young population ($`\tau 10^6`$ yrs), its dynamical evolution can not be neglected if the number of sources in a limited volume (size $`1`$ kpc) has to be assessed. In fact, a NS with a typical velocity $`300`$ km s-1 travels a distance $`300`$ pc in its lifetime as a cooler. While not much on a Galactic scale, such a displacement is non-negligible in evaluating the cooler population in the Solar proximity. For this reason we account for the dynamical evolution of NSs in the Galactic potential in our model.

Our calculation proceeds in three steps. First, a spatial distribution of the progenitors is selected and the ratio of NSs born in the Galactic plane to those originating from the Gould Belt is fixed, together with the birth rates (Sect. 2.1). Second, the dynamical evolution is followed assuming that NSs at birth receive a kick velocity drawn from a prescribed velocity distribution (see Sect. 2.2). Finally, we derive the X-ray flux as a function of age and position on the basis of an updated set of cooling curves (Sect. 2.3), and translate this in ROSAT count rate for an assumed model of the interstellar absorption (Sect. 2.4).

We varied the parameters of our model (NS formation rate, NSs mass spectrum, time of calculation, time step, spatial distribution). In the following subsections we focus on those parameters corresponding to the results presented in Sect. 3.

### 2.1 Initial spatial distribution

In our picture NSs are born in the Galactic disk and in the Gould Belt with a constant rate over the entire calculation (see Fig. 1 for a schematic view). This is a reasonable approximation since we are dealing with young objects, less than a few Myrs old. The Gould Belt is modeled as a very thin disk of 500 pc in radius and inclination to the Galactic plane of 18 with its center situated at 100 pc from the Sun in the Galactic anticenter direction. The central region of the Gould Belt is devoid of massive stars (see Pöppel 1997 for detailed description of the Gould Belt, and Torra et al. 2000 for a shorter one). To mimic this, we assume that no NS is born in the central region of the Belt up to a distance of 150 pc from its center.

We assume that the birthrate per unit area is independent of position, i.e. it is constant both in the Belt and in the Galactic disk. The two (constant) values are different, the one referring to the Belt being larger. To derive the NS birthrate, we rely upon direct counting of massive stars which are doomed to end in a supernova event. This gives $`2.9\times 10^{11}`$ massive progenitors per yr per pc2 (Tammann et al. 1994) which implies a NS birthrate of 30 NSs Myr-1 in 0.6 kpc around the Sun, the region where data are available. In relating the NS and supernova birthrates we assumed that core collapse produces a neutron star in 90% of the cases. The above value gives a more reliable estimate of the local NS birthrate, at present, as compared with radio-pulsar and average supernova statistics.

To get the separate contributions to the NS birthrate from the Belt and from the disk for the region inside 0.6 kpc we proceed as follows. In one Myr, out of a total of thirty, twenty NSs are born in the Gould Belt according to the estimate of 20–27 supernova events given by Grenier (2000). The remainder is uniformly distributed in the Galactic disk in a ring extending from 100 to 600 pc, to account for the known fact that the closest Solar neighbourhood is underpopulated by massive stars (e.g. Maíz-Apellániz 2001). This is in agreement with the results of Torra et al. (2000), who estimate that about 2/3 of the massive stars inside 0.6 kpc belong to the Gould Belt, and the rest to the Galactic disk. Direct counting of massive progenitors within 0.6 kpc from the Sun accounts for all NSs originating in the Belt. However, the disk clearly extends well beyond 0.6 kpc from the Sun. Although NSs born far away are very dim their contribution may become significant at low fluxes ($`0.1\mathrm{cts}\mathrm{s}^1`$). For this reason we decided to include also NSs born in the Galactic disk from 0.6 up to 3 kpc from the Sun. The NS birthrate per unit area in this region is assumed to be again constant and coincident with that in the disk inside 0.6 kpc. With this choice we find an agreement with existing data on the supernova rate in a region of 1 kpc around the Sun (Tammann et al. 1994).

In summary, we adopt a uniform birthrate of 20 NS Myr-1 in the Gould Belt and 250 NSs Myr-1 in the Galactic plane up to a limiting distance of 3 kpc from the Sun. Our approach underestimates the number of distant ($`1`$ kpc) newborn NSs in the direction of the Galactic center because the rate of supernova events is likely to increase there. However, interstellar absorption towards the Galactic center is very high in the Galactic plane (see Sect. 2.4), and this hinders the detection of sources at low Galactic latitude. We note that the local overabundance of young NSs is probably due not only to the enhanced number of massive stars harbored in the Belt but reflects also the fact that the star formation history in the Belt favors an enhanced supernova rate in the present epoch.

### 2.2 Dynamical evolution

A detailed description of the evolutionary code may be found in Popov et al. (2000a, b). Typically we calculate about a thousand evolutionary tracks in each run (up to $`10^4`$ per production run) and then normalize our results to the actual number of NSs born in the volume ($`1000`$) during the 4.25 Myrs time interval, as discussed above. An evolutionary track is calculated for specified initial position and velocity of a NS. All tracks are used over their whole duration and are applied to NSs of different masses, i.e. for each track we have several different cooling histories, the number of which depends on the number of different masses used in the run. As in Popov et al. (2000a,b), each track actually represents a population of NSs of different masses continuously born at the specified location with a prescribed spatial velocity distribution, during the entire time interval of the calculation.

NSs are assumed to be born with a Maxwellian kick velocity distribution with an average of $`225`$ km s-1. Varying this parameter does not change our results significantly because we are probing young objects. The Galactic potential is assumed to be axisymmetric, resulting from the sum of three contributions: the disk, bulge and halo (see Paczynski 1990). The Sun is in the Galactic plane at a distance 8.5 kpc from the Galactic Center.

### 2.3 Cooling curves and flux calculation

To calculate the cooling of NSs we use the results of the St. Petersburg group (see Kaminker et al. 2002, and the review by Yakovlev et al. 1999). Cooling curves are for NSs of masses from $`1.1M_{}`$ to $`1.8M_{}`$ with a step of $`0.1M_{}`$ (see Fig. 2). Curves take into account all important processes of neutrino emission. The equation of state (EOS) used in Kaminker et al. (2002) was introduced by Prakash et al. (1988). More precisely, it corresponds to Model I of Prakash et al. for symmetry energy and compression modulus of saturated nuclear matter $`K=240`$ MeV. The maximum NS mass in this model is 1.977 $`M_{}`$. Since we base our calculation on a definite EOS, the star radius $`R`$ is known once the mass $`M`$ is fixed. For $`1.1M_{}<M<1.8M_{}`$ it is $`12.2\mathrm{km}<R<13.2\mathrm{km}`$, where $`R`$ is circumferential radius. The EOS used here corresponds to matter composed of neutrons, protons, and electrons (no muons, hyperons, and no exotic particles). Neutron superfluidity was ignored.

We start to compute the emitted flux when a NS has an age of 10000 yrs. This choice is motivated by the fact that the Vela pulsar (the youngest close-by thermally emitting NS) is slightly older than that and also because the NS birth rate used in our calculation corresponds to an event every $``$ 10 000 yrs in the $`1.7`$ kpc region around the Sun. For each NS the calculation is stopped when the temperature drops below $`10^5`$ K. This happens at an age of 4.25 Myrs for the lightest NSs ($`M=1.1M_{}`$) or less for more massive stars. NSs that cold could have been detected by ROSAT PSPC within a distance of only 10 pc with a count rate $`5\times 10^3`$ cts s-1. We assume that emission comes from the entire star surface. This appears reasonable in the light of the relatively low pulsed fractions ($`15`$-20%) detected so far in the majority of isolated NSs (see Haberl & Zavlin 2002). We do not take into account any reprocessing of the blackbody surface emission in the NS atmosphere.

NSs are expected to have a mass spectrum (see discussion in Woosley et al. 2002). In the calculations we present here NSs are taken to have a flat mass spectrum between 1.1 and $`1.8M_{}`$. The mass spectrum (and clearly the mass of each star) are assumed to be constant during the calculations. With the above values (maximum age 4.25 Myrs and time step 10000 yrs) each evolutionary track represents $`425\times 8`$ NSs of eight different masses born with a time step $`10^4`$ yrs at the same place with the same initial velocity in a period 4.25-0.01 Myrs ago.

As one can see from Fig. 2, the main contribution to X-ray bright sources comes from NSs with masses $`<1.5M_{}`$. Young isolated NSs which are observed in several SN remnants as compact X-ray sources also should be relatively low-mass objects (Kaminker et al. 2002). Observations of binary radio pulsars suggest that NSs masses are strongly peaked around 1.35 $`M_{}`$ (Thorsett & Chakrabaty 1999). We repeated our calculations using a similar distribution and found that it produces similar results, although the number of observable isolated NSs increases by $``$ 30 %. Calculations which take into account the realistic mass function of massive progenitors in the Solar vicinity will be the subject of a future paper.

### 2.4 ISM and absorption

Since young cooling NSs are expected to emit most of their luminosity at UV/soft X-ray energies ($`20200`$ eV, corresponding to temperatures $`10^5`$$`10^6`$ K), interstellar absorption plays a crucial role with respect to their observability. Any attempt to derive the number of observable cooling isolated NSs using the unabsorbed flux would result in a substantial overestimate.

For the ISM distribution we use the same prescription as in Popov & Prokhorov (1998). The Local Bubble is modeled as a sphere with a radius of 140 pc and ISM density of 0.1 cm-3. Typical column densities for sources inside the calculated volume are in the range $`N_\mathrm{H}10^{19}`$$`10^{21}`$ cm2. After the column density is evaluated for a current NS position, we calculate the unabsorbed flux corresponding to the temperature, radius of the NS and its distance from the Sun, and apply the standard procedure to derive the ROSAT PSPC count rate. Outside the Local Bubble in directions close to the Galactic plane, absorption starts to play a crucial role. That is why regions closer to the Galactic center ($`<7`$ kpc), where the NS formation rate should be higher than in the Galactic disk in the Solar vicinity, cannot add many sources to our sample. At $`N_\mathrm{H}=3\times 10^{21}`$ cm-2 even very young and hot low-mass NSs with $`kT0.1`$ keV cannot have been detected at a distance 1 kpc in the ROSAT Bright Survey (RBS) at the threshold of $`0.2`$ cts s-1.

Except for the Local Bubble (and even that in a quite simplified way, see Sfeir et al. 1999 for a more complete description), our model does not take into account small-scale irregularities in the ISM distribution. They can be important if one makes an attempt to produce a realistic map of RINSs distribution on the celestial sphere, but in the case of the all-sky averaged Log N – Log S distribution our approximation is adequate.

## 3 Results

As discussed in the previous section, the dynamical evolution of young neutron stars, together with the calculation of their X-ray flux, allows us to derive the Log N – Log S for these sources. In particular, having fixed the number of NSs originating in the Gould Belt and in the Galactic disk, we can compute the separate contributions of these two sub-populations to the total Log N – Log S. Our main results are presented in Fig. 3 where the total (disk+Belt) and the disk-only distributions are compared with observations. All curves refer to the whole sky, i.e. the angular coverage used here is the entire solid angle. The observed Log N – Log S has been derived from ROSAT data of the seven RINSs and six other close-by young isolated NSs (see Table 1). Contrary to RINSs, the latter sources exhibit a composite spectrum with a non-thermal high-energy tail superimposed on the thermal component. The total count rate has been used in this case, but we stress that the non-thermal contribution is sizeable only for the Vela pulsar and PSR 1929+10. Poissonian errors are assumed in both observations and simulations. The computed curves were generated with a sufficiently large number of stars to reduce the statistical noise which is particularly severe at large flux (see §2.2). Statistical errors are shown for the observational points only, being too small to be appreciated in the plot of the simulated distributions. Uncertainties arising from the assumed model parameters are more difficult to assess and could introduce non-negligible deviations. Test runs with different choices of key parameters show that differences are of the order of 50%. Symbols which show observational points (filled diamonds or open circles) correspond to the type of the faintest object (RINS or not) which contributes to the total number at the specified count rate.

Previous investigations have convincingly shown that cooling NSs with the same spatial density of ordinary radio-pulsars cannot explain the observed Log N – Log S distribution for RINSs (Neuhäuser & Trümper 1999, Popov et al. 2000b). In particular, it was stressed that an additional contribution to the local population of young cooling NSs should be invoked to explain the relatively large number of bright sources at fluxes $`0.1`$ cts s-1. The Gould Belt is a natural candidate to provide the missing NSs in the Solar proximity. The population synthesis calculations presented here strongly support this claim and account not only for the observed distribution of RINSs but also for that of other young, close-by isolated NSs observed by ROSAT.

As can be seen from Fig. 3, the Gould Belt provides the major contribution to the local population of young cooling NSs and the theoretical prediction gives a very good fit to the data, within statistical errors. The contribution of NSs born in the Galactic disk is not very important at relatively large fluxes (it leaves rooms for only a few sources with count rate $``$ 0.1 cts s-1), but, as expected, becomes dominant at lower fluxes ($``$ 0.01 cts s-1) where far away stars contribute most. Fig. 3 clearly shows that the curve referring to the Galactic disk alone is always below the observed points with high statistical significance and below the RBS limit. The worsening of the agreement at low fluxes ($``$ 0.05 cts s-1) is to be attributed to the incompleteness of the X-ray sample. The deficit of very bright objects ($`1`$ cts s-1) may be attributed to chance statistics.

Our calculations show that there can be about $`1050`$ unidentified isolated NSs in the ROSAT All-Sky Survey (RASS) at a limiting flux of $`0.015`$ cts s-1 depending on parameters of the model. The number of INSs in the RASS/BSC has been recently investigated by Rutledge et al. (2003), who found that at most 67 sources could have been detected at the 0.05 cts s-1 level and have escaped identification. Our results are well within this limit. Also there may be a few unidentified RINSs at fluxes $`0.1`$ cts s-1 at low Galactic latitudes (see also Schwope et al. 1999). Most sources should be observed at $`\pm `$ 20 from the Galactic plane towards the directions of lower absorption. Some of them may be Geminga-like objects with counterparts among unidentified gamma-ray sources (also connected with the Gould Belt, see Grenier 2000). Identification of these objects may prove important for constraining cooling models and the NS mass spectrum.

Absorption, the flat geometry of NS initial distribution and the finite extension of the Gould Belt naturally explain the very flat (slope flatter than -1) Log N – Log S curves in Fig. 3. Our model contains some degrees of freedom (e.g. the NS mass spectrum, details of the formation rate and surface emission, a simple blackbody was used here) that can be varied. We nevertheless believe that even this simple picture (albeit based on realistic assumptions) gives a quite satisfactory explanation of the observed properties of isolated NSs in the Solar proximity. Preliminary calculations which take into account realistic mass spectrum of NS progenitors, atmospheric effects and different spatial and velocity distributions are in progress and show no qualitative changes.

## 4 Discussion and conclusions

At present about 20 nearby (distance $`<1`$ kpc), young (age $`<4.25`$ Myrs) isolated NSs are known (see Table 1). The local NS population includes objects with different properties: radio-quiet, thermally-emitting NSs (the seven RINSs), Geminga and the Geminga-like object 3EG J1835+5918 (these are probably active pulsars with radiobeams missing the Earth), radio-pulsars with observed thermal X-ray emission (PSR 1929+10 and the “three musketeers”: the Vela pulsar, PSR 0656+14, PSR 1055-52) and seven other pulsars which are not detected in X-rays. The latter objects are relatively old in comparison with the others and lie further away, as can be seen from Table 1 where the wide gap in the estimated age of PSR 1055-52 and PSR J0056+4756 is apparent. For ages of a few Myrs, even very low-mass (and hence slow-cooling) NSs are too cold by now to have been detected with ROSAT at distances $`0.5`$ kpc. The case of PSR 1929+10 is intermediate since it is $`3`$ Myrs old but relatively close ($`300`$ pc) and its X-ray emission is mainly non-thermal.

As has been discussed by Neuhäuser & Trümper (1999) and Popov et al. (2000b), the number of X-ray bright INSs is too large to be explained in terms of the average density of radio-pulsars in the Solar neighborhood. This motivated the suggestion that a sizeable fraction of NSs never become active radio emitters, as suggested at the same time by considerations on young NSs in supernova remnants (Gotthelf & Vasisht, 2000) albeit the origin of the local NS population was left open. In this paper we have shown that thermally-emitting INSs are naturally explained as cooling NSs born in the Gould Belt. The relatively high local spatial density of young NSs is then due to the large number of massive progenitors in the young stellar associations which constitute the Gould Belt. Our analysis lends support to the idea that RINSs (the “magnificent seven”) represent the slowly cooling NSs and hence a clean sample of NSs with $`M1.3M_{}`$ (Kaminker et al. 2002).

The computed Log N –Log S distribution for X-ray thermally-emitting INSs born in the Gould Belt and (to a lesser extent) in the Galactic disk accounts for all bright INSs in the Solar vicinity (see Fig. 3) and leaves room for at most 1-2 undetected sources above $`1`$ cts s-1. The absence of sources brighter than RX J1856 (Log S$`>0.56`$) is consistent with our calculations at the 2 $`\sigma `$ level. Having neglected very young NSs (age $`<10000`$ yrs) does not change our conclusions significantly. Only young stars with $`M<1.35M_{}`$ (about 1/2 of their total number) have temperatures about $`2\times 10^6`$ K (see Fig. 2), which corresponds to a luminosity $`10^{34}`$ erg s-1. With our NS formation rate we expect to have one such NS within a distance of $`2.5`$ kpc. About half of the NSs born in that region should originate inside 1.7 kpc. Placing the source at this distance its unabsorbed flux is $`3\times 10^{11}`$ erg cm-2 s-1. With a typical column density of $`10^{21}`$ cm-2, Web-PIMMS gives $`1.3`$ cts s-1 for ROSAT PSPC. We can then conclude that stars younger than 10000 yrs can contribute only to fluxes $`<1.3`$ cts s-1. This translates into a $`20\%`$ difference in the Log N - Log S curve shown in Fig. 3 in the worst case. Typical uncertainties of the present model are at about the same level.

The local ($`1`$ kpc) INS population was estimated to comprise total about one hundred objects younger than $``$4 Myr, taking into account that some NSs born inside 1 kpc can leave the local space in their lifetime. These NSs are not detected as radio pulsars, but tens of them could be identified in ROSAT surveys as dim sources. The beaming effect can be responsible only for part of these young NSs to be radiosilent. According to Tauris & Manchester (1998) in fact, about 50-70% of young pulsars are not visible from the Earth (see also Brazier & Johnston 1999). In our model this would imply the presence of at least $`30`$ active pulsars with an age $`4`$ Myrs in the Solar vicinity. However, only about 1/3 of such pulsars are observed (see again Table 1). A bimodal velocity distribution, with the high charcteristic velocity $`500`$ km s-1 (Arzoumanian, Chernoff & Cordes 2002), would reduce the discrepancy since a number of young pulsars leave the local volume considered here. However, even barring observational bias against low-luminosity radio-pulsars, our model seems to support the argument by Gotthelf & Vasisht (2000), that at least half of the observed young neutron stars follow an evolutionary path quite distinct from that of the Crab pulsar.

###### Acknowledgements.

We want to thank Dmitry Yakovlev for putting his cooling model to our disposal and for his invaluable help with it. We are also indebted to Vasily Beskin, Matteo Chieregato and Andrea Possenti for many useful discussions. The work of SP was supported by the Russian Foundation for Basic Research (RFBR) grant 02-02-06663 and by RSCI. SP thanks the Universities of Insubria at Como and Milano-Bicocca, where part of this investigation was carried out, for hospitality. The work of MC, AT and RT was partially supported by the Italian Ministry for Education, University and Research (MIUR) under grant COFIN-2000-MM02C71842.


text:
# 𝐾_𝐿→𝛾𝜈𝜈̄ decay beyond the standard model

## Abstract

The decay $`K_L\gamma \nu \overline{\nu }`$ is investigated beyond the standard model. Interestingly, the upper limit of the $`CP`$-conserving and $`CP`$-violating branching ratios of the decay, induced from the possible extensions of the standard model, would be larger than the corresponding branching ratios given in the standard model respectively, and it is expected that the $`CP`$-violating part could be enhanced.

USTC-ICTS-02-05

Rare kaon decays play an important role in studying flavor-changing neutral currents (FCNC) and $`CP`$-violating phenomena in modern particle physics . Within the context of the standard model, they are suggested to test the Cabibbo-Kobayashi-Maskawa (CKM) paradigm. They also provide an ideal place to search for new physics beyond the standard model . Even the new physics may appear in $`B`$-meson decays, it is interesting to study the footprint it leaves in rare kaon decays and especially the $`sd`$ weak transition.

The decays $`K^+\pi ^+\nu \overline{\nu }`$ and $`K_L\pi ^0\nu \overline{\nu }`$ are both governed by the same $`sd\nu \overline{\nu }`$ transition. They are dominated from the short-distance loop contributions containing virtual heavy quarks. It is believed that the long-distance contributions in these decays are much smaller than the short-distance ones, and could be negligible . Due to the absence of possible large theoretical uncertainty from the long-distance, these two modes are very useful both to test the standard model and to explore the new physics. Experimentally, the decay $`K^+\pi ^+\nu \overline{\nu }`$ has been measured by E787 group at BNL

$$B(K^+\pi ^+\nu \overline{\nu })=(1.57_{0.82}^{+1.75})\times 10^{10},$$
(1)
which is consistent with the prediction of the standard model

$$B(K^+\pi ^+\nu \overline{\nu })_{SM}=(0.72\pm 0.21)\times 10^{10}.$$
(2)
The amplitude of $`K_L\pi ^0\nu \overline{\nu }`$ in the standard model is proportional to the CKM factor $`Im(V_{td}^{}V_{ts})`$ and the decay branching ratio is found to be at the level of $`10^{12}`$ . Whereas the current experiment limit is less than $`5.9\times 10^7`$, and this work is ongoing. But from the experimental point of view, it is difficult to find the event, because all the final-state particles are neutral, only the $`2\gamma `$’s from $`\pi ^0`$ can be detectable.

An alternative way is to search the decay $`K\pi \pi \nu \overline{\nu }`$ . However, the decay branching ratio is small and the background for $`\pi \pi `$ is large . So another good choice is to study the decay of $`K_L\gamma \nu \overline{\nu }`$, where there is only one photon at the final state.

The total branching ratio of the decay $`K_L\gamma \nu \overline{\nu }`$ can be divided into CP-conserving part and CP-violating part . The contribution within the standard model to this decay has been investigated by some authors, and it is believed that the mode is short-distance dominated . The purpose of this paper is to examine the possible contribution to this decay from the new physics scenarios beyond the standard model. It will be shown below that there exists the significant model-independent limit, from which the branching ratio of $`CP`$-conserving part could be up to the same order of magnitude as the one predicted within the standard model while the $`CP`$-violating part could be of more than one order of magnitude lager than the one given in the standard model. As an interesting example, we also consider the possible dominant supersymmetric contributions to this decay mode.

The effective Hamiltonian for $`sd\nu \overline{\nu }`$ transition generally takes the form

$`H_{eff}={\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{2\pi \mathrm{𝚜𝚒𝚗}^\mathrm{𝟸}\theta _W}}W_{ds}\overline{d}\gamma _\mu (1\gamma _5)s\overline{\nu }\gamma ^\mu (1\gamma _5)\nu +H.c.,`$ (3)
where the short-distance physics is lumped in $`W_{ds}`$. In the standard model $`W_{ds}`$ is dominated by penguin and box diagrams with intermediate charm and top quark

$`W_{ds}^{SM}={\displaystyle \underset{l=e,\mu ,\tau }{}}[\lambda _cX_{NL}^l+\lambda _tX(x_t)]`$ (4)
where $`x_t=m_t^2/M_W^2`$, and $`\lambda _i=V_{is}^{}V_{id}`$. The functions $`X_{NL}^l`$ and $`X(x_t)`$ correspond to charm and top contributions in the loops with the next-to-leading logarithmic approximation respectively. It has been shown explicitly in Ref.  that the corresponding top contribution can be written as $`X(x_t)=\eta _tX_0(x_t)`$, with the QCD-uncorrected top quark contribution

$`X_0(x_t)={\displaystyle \frac{x_t}{8}}\left[{\displaystyle \frac{2+x_t}{1x_t}}+{\displaystyle \frac{3x_t6}{(1x_t)^2}}\mathrm{𝚕𝚗}(x_t)\right],`$ (5)
and the QCD correction factor $`\eta _t=0.994`$.

Also from table I in Ref., one has

$$X_{NL}^{e,\mu }=11.00\times 10^4,$$
(6)
and

$$X_{NL}^\tau =7.47\times 10^4.$$
(7)
with the central value of the QCD scale $`\mathrm{\Lambda }=\mathrm{\Lambda }_{\overline{MS}}^{(4)}=325\pm 80MeV`$, and the charm quark mass $`m_c=\overline{m_c}(m_c)=1.30\pm 0.05GeV`$.

From the effective Hamiltonian in eq.(3), one can evaluate the hadronic matrix element $`\gamma |J_\mu |K^0`$, where $`J_\mu =\overline{d}\gamma _\mu (1\gamma _5)s`$. Then the corresponding matrix elements are parameterized as

$`\gamma (q)|d\gamma ^\mu \gamma _5s|K^0(p+q)`$ $`=`$ $`e{\displaystyle \frac{F_A}{M_k}}\left[ϵ^\mu (pq)(ϵ^{}p)q^\mu \right],`$ (8)
$`\gamma (q)|d\gamma ^\mu s|K^0(p+q)`$ $`=`$ $`ie{\displaystyle \frac{F_V}{M_k}}\epsilon ^{\mu \alpha \beta \gamma }ϵ_\alpha ^{}p_\beta q_\gamma ,`$ (9)
where $`q`$ and $`p+q`$ are photon and $`K`$-meson four momenta, $`F_A`$ and $`F_V`$ are form factors of axial-vector and vector respectively, and $`ϵ^{}`$ is the photon polarization vector. The form factors $`F_V`$ and $`F_A`$ have been evaluated in  using light front quark model. We find, it is enough to use momentum independent form factors, $`F_A(0)=0.0429`$ and $`F_V(0)=0.0915`$ given by the authors of  to illustrate the numerical estimates in the present paper.

Following the similar way to study $`K_{L,S}\pi \pi \nu \overline{\nu }`$ in Ref. , here we also consider two representative cases to probe possible new physics effects of $`K_L\gamma \nu \overline{\nu }`$ : one is a model-independent way; the other is a special case–low energy supersymmetry, which however represents one of the most interesting and consistent extensions of the standard model.

i) Effective flavor-changing neutral currents (FCNC) interaction. As formulated by Nir and Silverman , it reads

$`^{(z)}={\displaystyle \frac{g}{4\mathrm{𝚌𝚘𝚜}\theta _𝚆}}U_{ds}\overline{d}\gamma _\mu (1\gamma _5)sZ^\mu .`$ (10)
Combining with the coupling of the $`Z`$ boson to neutrino-antineutrino pair, one finds that

$`W_{ds}^{NP}={\displaystyle \frac{\pi ^2}{\sqrt{2}G_FM_W^2}}U_{ds}=0.93\times 10^2U_{ds}`$ (11)
as the new piece of $`W_{ds}`$ in the effective Hamiltonian contributing to the basic transition of $`sd\nu \overline{\nu }`$.

The upper bounds on $`U_{ds}`$ have been determined by other processes involving $`K`$-meson and were summarized in Ref. . After translating them into the bounds on $`W_{ds}`$, we obtain

$`|Re(W_{ds})|`$ $``$ $`0.93\times 10^3,`$ (12)
$`|W_{ds}|`$ $``$ $`2.8\times 10^3,`$ (13)
$`|Re(W_{ds})Im(W_{ds})|`$ $``$ $`1.1\times 10^5,`$ (14)
$`|Im(W_{ds})|`$ $``$ $`0.93\times 10^3.`$ (15)
The bound on $`|W_{ds}|`$ can be improved using the most recent measurement  of $`K^+\pi ^+\nu \overline{\nu }`$, which is $`|W_{ds}|=1.0_{0.31}^{+0.46}\times 10^3`$.

ii) Supersymmetry: Possible interesting supersymmetric effects arise from penguin diagrams involving charged-Higgs plus top-quark intermediate states or squark and chargino intermediate states. These give additional pieces to the effective Hamiltonian of the form

$`W_{ds}^{NP}=\lambda _t{\displaystyle \frac{m_H^2}{M_W^2\mathrm{𝚝𝚊𝚗}^\mathrm{𝟸}\beta }}H(x_{tH})+{\displaystyle \frac{1}{96}}\stackrel{~}{\lambda _t},`$ (16)
where tan$`\beta `$ is the ratio of the two Higgs vacuum expectation values and $`x_{tH}=m_t^2/M_{H^\pm }^2`$. The quantity $`H(x)`$ is given by

$`H(x)={\displaystyle \frac{x^2}{8}}\left[{\displaystyle \frac{\mathrm{𝚕𝚗}x}{(x1)^2}}+{\displaystyle \frac{1}{x1}}\right].`$ (17)
The definition of $`\stackrel{~}{\lambda _t}`$ is shown in Ref. . The $`\stackrel{~}{\lambda _t}`$ part in eq. (16), which comes from chargino-squark diagrams, represents the dominant supersymmetric effect  \[Only small value of tan$`\beta `$ might enhance the first term in eq. (16), which however is not in favor with the present phenomenologies\]. The bound on $`\stackrel{~}{\lambda _t}`$ can be imposed by similar considerations that were used for $`W_{ds}`$. In Ref. , the upper limits $`|Re\stackrel{~}{\lambda _t}|0.21`$ and $`|\stackrel{~}{\lambda _t}|0.35`$ have been obtained from the observed branching ratios of decays $`K_L\mu ^+\mu ^{}`$ and $`K^+\pi ^+\nu \overline{\nu }`$. Using the most recent experiment data , the last limit could be updated as $`|\stackrel{~}{\lambda _t}|0.14`$.

Using the effective Hamiltonian for $`K^0\gamma \nu \overline{\nu }`$ in eq.(1) together with $`K_LK_2=(K^0+\overline{K^0})/\sqrt{2}`$, the amplitude of $`K_L\gamma \nu \overline{\nu }`$ can be divided into the $`CP`$-conserving and $`CP`$-violating parts as follows

$`(K_L\gamma \nu \overline{\nu })=_{CPC}+_{CPV},`$ (18)
where

$`_{CPC}`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{2\pi \mathrm{𝚜𝚒𝚗}^\mathrm{𝟸}\theta _𝚆}}{\displaystyle \frac{2}{\sqrt{2}}}Re(W_{ds})ϵ^{\mu \nu \alpha \beta }{\displaystyle \frac{F_V}{M_K}}ϵ_\mu ^{}q_\alpha p_\beta ^{}\overline{u(p_\nu )}\gamma _\nu (1\gamma _5)v(p_\nu ),`$ (19)
and

$`_{CPV}`$ $`=`$ $`{\displaystyle \frac{G_F}{\sqrt{2}}}{\displaystyle \frac{\alpha }{2\pi \mathrm{𝚜𝚒𝚗}^\mathrm{𝟸}\theta _𝚆}}{\displaystyle \frac{2}{\sqrt{2}}}Im(W_{ds}){\displaystyle \frac{F_A}{M_K}}ϵ_\mu ^{}(p^{}qg^{\mu \nu }+p^\mu q^\nu )\overline{u(p_\nu )}\gamma _\nu (1\gamma _5)v(p_\nu )`$ (20)
In $`K_L`$ rest frame the partial decay rate of $`K_L\gamma \nu \overline{\nu }`$ is given by

$`d^2\mathrm{\Gamma }={\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle \frac{1}{8M_K}}||dE_\gamma dE\nu .`$ (21)
By integrating the $`E_\gamma `$ and $`E_\nu `$, we obtain

$`\mathrm{\Gamma }_{CPC}={\displaystyle \frac{\alpha }{15}}\left({\displaystyle \frac{G_F\alpha }{16\pi ^2\mathrm{𝚜𝚒𝚗}^\mathrm{𝟸}\theta _W}}\right)^2|Re(W_{ds})|^2F_V^2M_K^5,`$ (22)
$`\mathrm{\Gamma }_{CPV}={\displaystyle \frac{\alpha }{15}}\left({\displaystyle \frac{G_F\alpha }{16\pi ^2\mathrm{𝚜𝚒𝚗}^\mathrm{𝟸}\theta _W}}\right)^2|Im(W_{ds})|^2F_A^2M_K^5.`$ (23)
Using $`\mathrm{𝚜𝚒𝚗}^2\theta _W=0.23`$ and $`M_K=0.5GeV`$ the decay branching ratios are found to be

$`B(K_L\gamma \nu \overline{\nu })_{CPC}`$ $`=`$ $`6.7\times 10^8|Re(W_{ds})|^2,`$ (24)
$`B(K_L\gamma \nu \overline{\nu })_{CPV}`$ $`=`$ $`1.4\times 10^8|Im(W_{ds})|^2.`$ (25)
The decay branching ratios of $`K_L\gamma \nu \overline{\nu }`$, within the standard model, have been evaluated by Geng, Lih, and Liu

$`B(K_L\gamma \nu \overline{\nu })_{CPC}^{SM}=1.0\times 10^{13},`$ (26)
and

$`B(K_L\gamma \nu \overline{\nu })_{CPV}^{SM}=1.5\times 10^{15}.`$ (27)
For the representative scenarios of physics beyond the standard model, the model-independent upper limits of the $`|Re(W_{ds})|`$ and $`|Im(W_{ds})|`$ have been given in eq.(15). Thus the maximal values of the branching ratios of the decay could be

$`B(K_L\gamma \nu \overline{\nu })_{CPC}=1.8\times 10^{13},`$ (28)
and

$`B(K_L\gamma \nu \overline{\nu })_{CPV}=0.42\times 10^{13}.`$ (29)
Note that a factor 3 has been multiplied in deriving eqs. (28) and (29) since all three flavor neutrinos should be taken into account. Also from the previous discussions on the limits of $`\stackrel{~}{\lambda _t}`$, it is easy to see that supersymmetric extensions of the standard model could allow $`|Re(W_{ds})|`$ and $`|Im(W_{ds})|`$ to reach the upper limit $`0.93\times 10^3`$ without conflict with other constraints.

In conclusion, it is shown that the upper limit value of $`CP`$-conserving branching ratio in eq.(28) can be of the same order of magnitude as the result given in the standard model, while the upper limit value of $`CP`$-violating branching ratio in eq.(29) is nearly a factor of 30 larger than the value within the standard model. Thus it seems that there exists room to look for new physics beyond the standard model by observing the decay $`K_L\gamma \nu \overline{\nu }`$ although this is not a very easy task. In addition to investigations on the decays of $`K\pi \nu \overline{\nu }`$ and $`K\pi \pi \nu \overline{\nu }`$, it is therefore expected that both experimental and theoretical studies of $`K_L\gamma \nu \overline{\nu }`$ could provide complementary information on the basic weak transition $`sd\nu \overline{\nu }`$.

ACKNOWLEDGMENTS

This work is partially supported by NSF of China 90103002 and 19905008.