text:
# High Speed Phase-Resolved 2-d 𝑈𝐵𝑉 Photometry of the Crab pulsar Based on observations taken at SAO, Karachai-Cherkessia, Russia

## 1 Introduction

As the youngest known isolated rotation powered neutron star, and having the highest known spin-down flux density, the Crab pulsar is not surprisingly the most efficient source of high energy magnetospheric emission. All theoretical models attempting to elucidate the processes of nonthermal emission must satisfy the empirical constraints provided by the extensive observational datasets available for this object from $`\gamma `$-rays to the infrared.

Contemporary theoretical frameworks in existence to explain the nonthermal emission from this and other isolated rotation-powered neutron stars may be broadly divided into two schools. The first places the source of emission close to the polar cap region, e.g. Sturrock (1971), and the second that place the source at a considerable distance above the neutron star surface, in the outer magnetosphere (e.g. Cheng, Ho & Ruderman 1986). In recent times, the development of these models has been to some extent restricted to the problem of $`\gamma `$-ray emission. These efforts have been predominantly analytical in nature, as in this energy regime primary emission is being sampled, which avoids complications associated with the lower energy sources of emission. In these latter cases, the emission is secondary in nature and a result of interactions that require numerical representation. Thus the polar cap school (e.g. Daugherty & Harding 1982) concluded with correct flux estimates but poor light curves, and the outer gap school (Cheng, Ho & Ruderman 1986) with reasonable light curves and approximate flux distributions but requiring extraordinary lines-of-sights and magnetospheric physics to justify the emission.

The group at Stanford University lead by R. Romani is generally credited by being the first to rigourously put the existing models to the test numerically. Their conclusions in many ways matched the intuitive premonitions regarding the two model frameworks, with perhaps the ‘evolved’ outer gap model of Cheng, Ho & Ruderman (1986) being the more likely candidate (Chiang & Romani 1992; Romani & Yadagiroglu 1995). Recently Romani (1998) has indicated that it is only by both the development of more advanced numerical simulations that incorporate the full emission physics (from $`\gamma `$-rays to the IR) in a truly self-consistent manner and the accurate characterisation of pulsar light curves in the lower energy regimes that one can hope to disentangle the problem of nonthermal pulsar emission.

Optically, spectroscopic, integrated and high speed single-pixel photometry have produced some critically important results in this regard. However, true phase-resolved data acquisition with accurate isolation of the pulsar’s emission from that of the nebula requires high time resolved ($``$ 1 $`\mu `$s) 2-d photometry to confirm unambiguously earlier conclusions and to determine the true nature of the magnetosphere’s emission in the range $``$ 1 - 10 eV.

In a previous paper (Golden et al. 2000), we presented an analysis of such observations of the Crab pulsar in which we isolated and theoretically speculated upon the unpulsed component of the pulsar’s optical emission. These observations were made in January 1996 with the TRIFFID 2-d high speed photometer in three colourbands ($`UBV`$) using the 6m BTA of the Special Astrophysical Observatory in the Russian Caucasus. Here we present a thorough photometric and temporal analysis of the same data. Following a review of the pulsar’s relevant empirical and theoretical characteristics, we briefly outline the observational and data reduction methodology - the reader is encouraged to consult Golden et al. (2000) for a more rigorous treatment of this process. We then detail the results of the photometric and temporal analysis, and conclude with a discussion on the implications of these results for current thinking on this pulsar, with wider ramifications for pulsar emission theory in general.

## 2 Review

The Crab pulsar is identified as the south-west member of the double stars located in the central region of the Crab Nebula (Cocke et al. 1969). Radio observations, initially by Staelin $`\&`$ Reifenstein (1968), and regularly ever since, have provided perhaps one of the most detailed pulsar rotational datasets to date. The pulsar has a $`P`$ $``$ 33 ms and $`\dot{P}`$ $``$ 4.211$`\times 10^{13}`$ s$`\mathrm{s}^1`$ suggesting a canonical age of $`\tau `$ $``$ 1260 years and surface magnetic field $`B_{surface}`$ $``$ 3.7 $`\times 10^{12}`$ Gauss.

Historically, high speed optical photometry has been generally dominated by the use of single pixel detectors, with variable time resolution. Following the detection of Cocke et al. (1969), numerous targetted observations followed (Kristian et al. 1970; Wampler et al. 1969; Warner et al. 1969; Becklin et al., 1973; Muncaster and Cocke, 1972; Cocke and Ferguson 1974; Groth 1975a,b) with the pulsar generally observed in the $`BVRI`$ wavebands at typical accurate time resolutions of order $``$ ms, and the datasets ‘added’ together in a least-squares fashion - there being no common stabilised ‘clock’ from which to time the observations, nor a definitive radio ephemeris with which to accurately obtain in modulo lightcurves at that time.

The development of more ”temporally consistent” instrumentation in conjunction with other advances in technology allowed for more detailed observations to be taken. Peterson et al. (1978) applied for the time rather novel techniques in the image processing of data obtained via the use of a 6.2 ms time resolved 2-dimensional Image Photon Counting System camera. Significant emission from the unpulsed presumed ‘off’ component was noted in the $`U`$ and $`B`$ bands, with Peterson et al. estimating that it was some 2% of the total pulsed emission.

At an effective time resolution of 20 $`\mu `$s, Smith et al. (1978) were able to make unprecedented observations via a S20 photocathode of the Crab’s main pulse region, indicating that the peak had a plateau of emission $``$ 100 $`\mu `$s in duration. Photon binning at 1 kHz allowed for phase resolved tests to be made of the counting statistics incident on the detector. They yielded statistical distributions that did not deviate significantly from random fluctuations expected, making an allowance for variations in sky transparancy.

Subsequent high speed observations incorporating an appropriate polarisation set-up yielded phase resolved polarimetry (Jones et al. 1981, Smith et al. 1988). Under analysis and corrected for interstellar polarisation, these datasets show evidence for sharp swings in polarisation angle around both principle peaks, in addition to some form of evolution of the polarisation in the bridge area. Remarkably, the results indicated substantial polarisation associated with what was always assumed to be ‘off’ or at the very least, the pulsar at a minimum, as noted earlier by Peterson et al. (1978).

The High Speed Photometer (HSP) on board the Hubble Space Telescope ($`HST`$) was used to observe the Crab pulsar in October 1991 and again in January 1992, as reported by Percival et al. (1993). The passbands used were F400LP and F160LP, which crudely approximate to standard $`V`$ and $`UV`$. Although the on-board clock is accurate to 10.74 $`\mu `$s, accuracy to UTC was typically $``$ $`\pm `$ 1.05 ms per observing run. Each resulting light curve had an effective temporal resolution of 21.5 $`\mu `$s.

The pertinent results of Percival et al. (1993) are as follows. They estimated that the main pulse in the visible leads the radio by -1.47 $`\pm `$ 1.05 ms, and for the UV, -0.91 $`\pm `$ 1.05 ms, with both $`UV`$ and $`V`$ main peaks offset by 0.56 $`\pm `$ 0.92 ms. Via autocorrelation and least-squares analysis of the light curves, the Full Width Half Maximum (FWHM), peak separation and peak ratios were determined. Despite the summing and normalising of the individual light curves which might be expected to introduce an element of bias in any analysis, Percival et al. (1993) obtain consistent and reasonable estimates for these parameters. The implication is that there is a trend whereby the pulse peaks ‘tighten’ and come together as a function of energy. Resolution of the main peak cusp was attempted via a least-squares polynomial fit, and was found to be $``$ 150 $`\mu `$s in extent. Despite a number of caveats outlined previously Percival et al. (1993) can claim to have shown to a reasonable degree that

1. The peaks’ FWHM and separation contract as a function of energy.
2. The $`UV`$ and $`V`$ main peaks appear to be coincident to within $``$ 1 ms.
3. The $`UV`$ and $`V`$ main peaks show signs of leading the radio main peak.
4. A plateau of emission is not definitively resolvable to within 150 $`\mu `$s.
5. The photon statistics would appear to be Poissonian throughout the light curve.

Later observations, using the $`HSP`$ in Polarisation mode, were made of the pulsar (Smith et al. 1996). Qualitatively the results were similar to earlier optical observations. The conlusions reached were no different to the original 1988 paper, other than that any model of the emission regions must be essentially wavelength independent.

Following on from $`IR`$ band observations of Ransom et al. (1994), Eikenberry et al. (1996) made observations in the $`J`$, $`H`$ and $`K`$ bands using a single-pixel Solid State Photomultiplier, set-up to record using 20 $`\mu `$s time bins. Light curves were then obtained in the usual manner.

The analysis proceeded along three avenues. The first was determination of the Half Width Half Maximum (HWHM), essentially a truncated FWHM, with each side of the FWHM bisected by the maximum of the peak in question, yielding a ‘leading’ and ‘trailing’ HWHM component. The second was the determination of phase resolved spectra, and the third an in-depth examination of the colour variation around the main peak. Monte Carlo simulations were used to build up a statistical picture of the optimum fits to the light curve data.

This analysis showed clear evidence that the peaks (morphologically and in terms of flux) change as a function of energy throughout the light curve, and also that time limit over which such rapid spectral changes occured was $``$ 180 $`\mu `$s. Assuming emission occurs in the vicinity of the light cylinder, then one can estimate the typical coherence dimensions for the region, being less than or equal to 1.8 $`\times 10^4`$ s $`\times `$ $`c`$ = 54km - indicating that the emission region is consistent with a localised region of the magnetosphere.

In a subsequent paper (Eikenberry & Fazio 1997) these techniques were applied to a range of light curves, from $`\gamma `$-rays (Ulmer et al. 1995), X-rays (reconstructed from the $`ROSAT`$ archives), $`UV/V`$ (Percival et al. 1993) and $`KHJ`$ (from previous work). The overall energy range spanned 0.47 eV to 10 MeV. The principle intention was to further test the hypothesis that the actual peaks themselves altered morphologically with energy. The actual mechanics of how the datasets were analyzed was similar to the earlier analysis of Eikenberry et al. (1996) \- it will suffice to note their more pertinent conclusions;

1. The ratios of the flux of the bridge and peak 2 to peak 1 show an overall trend to increase with energy, but reverses at low energies ( $`<`$ 1eV ).
2. The peak-to-peak separation shows a ‘nonlinear’ trend to decrease with energy over the range, which is consistent with energy stratification of the emission region.
3. The peaks FWHM show neither a definite trend nor pattern with energy, yet display significant variability over the range.
4. There is a strong energy dependency as regards the evolution of the HWHM for leading and trailing edges of both peaks, with some evidence of a continous functionality over the range - except at the $`IR`$ region. The HWHM behaviour explains the rather puzzling FWHM behaviour.
5. Peak 2 appears to change profoundly as regards morphology between the $`IR`$-optical and X-ray-$`\gamma `$-ray range, showing a fast rise and slow fall for the former, and the converse for the latter. This is not explicable by theory.
6. Many of the pulse shape parameters show maxima or minima at the energies of 0.5-1eV in the $`IR`$ regime, suggesting some interesting phenomena occuring in this waveband.

Eikenberry & Fazio (1997) concluded that the many unusual energy dependent morphological phenomena that are observed under this analysis were incompatible with current model frameworks - in particular the single-pole emitting outer-gap (Romani & Yadigaroglu 1995). In this model, the emission originates from a topologically extended region that effectively maps out a hemisphere about the open field region. Consequently, with the emission profile a result of viewing geometry and crossing caustics to the line of sight, one would expect either no evident effect, or a common effect that should manifest itself in a smooth and continuous manner at all such energies, as the effect would be essentially the averaging of many seperate spatially emission sources. These expectations contrast with the empirically determined results.

## 3 Precession Issues

Free precession of an isolated pulsar, most likely due to some moment of inertia inequality, possibly as a result of irregular distributions of neutron star matter on its surface, or perhaps as a result of internal interactions and effects, would manifest itself as an appropriate modulation in the observed temporal signal. Previous reports by groups examining the variations in the ratios of peak 1 to peak 2 in $`\gamma `$-rays have generally indicated periods of the order of years (Kanbach et al. 1990; Nolan et al. 1993; Ulmer et al. (1994). Previous reports of something akin to a precessional like signal from the Crab pulsar in the radio have tended to suggest a precessional period of order months/years (e.g. Jones 1988).

At higher frequencies, evidence of precession could provide both direct constraints to the condensed matter equation of state, and imply the emission of gravitational radiation. It is well known that a precessing gravitationally compact body rotating at high speeds would be expected to radiate gravity waves, at a frequency twice its rotational frequency (e.g. Misner, Thorne & Wheeler, 1971). Various other emission modes are theoretically possible. Planned ground based gravitational detectors are being designed that could have the capability to detect displacement amplitudes of order $`h`$ $``$ $`10^{23}`$. However do this, these detectors must be ‘tuned’ to a particular tight frequency band for up to 1 year to build up sufficient $`S/N`$. Clearly, one needs a guaranteed signal to lock on to - and pulsars may very well provide such a source.

Hence the great interest in the reported suggestions of a precession signal with a period of order 60 seconds from the Crab by Cadez and Galicic (1996). This was based on a time-series analysis of the High Speed Photometer data of Percival et al. (1993), and their own ground based observational CCD dataset. Despite nearly three decades of radio observations, no such modulation has been reported in the literature. Furthermore, an analysis of high speed optical data by Jones et al. (1980) with early photometers indicated a lack of optical source variations at the 1 % level over a ten year period. It is clearly important to try and test this with independent data, as a confirmation would have important possible implications for other fields, particularly that of gravitational wave detection.

## 4 January 1996 BTA Run - Observations in $`UBV`$ $`\&`$ Data Preparation

The details of these TRIFFID high speed observations of the Crab pulsar, and the subsequent preliminary analysis have been published elsewhere (Golden et al. 2000). We shall thus limit our description accordingly. The observations took place over 5 nights in January 1996 using the Special Astrophysical Observatory’s 6m telescope in the Caucasus. TRIFFID high speed camera used incorporated a 2-dimensional photon counting MAMA detector (Cullum (1990); Timothy & Bybee (1985)) with standard $`U`$, $`B`$ and $`V`$ filters. Timing was accurate to within 400 ns per 10 second period using a GPS receiver. There were 6 $`U`$ band exposures totalling 5546 seconds, 5 $`B`$ band exposures totalling 5750 seconds, and 9 $`V`$ band exposures totalling 6759 seconds.

The photon \[x,y,t\] data stream per dataset was corrected via a Weiner filter modified shift-and-add algorithm for telescope wobble and gear drift (Redfern et al. 1993). Processing the datasets in this way and incorporating the derived flat fields yielded full field images of the inner Crab nebula from which the pulsar and its stellar companions were registered. Figure 1 shows the integrated 96zs2.0.0 dataset following this technique. The average seeing for this dataset was $``$ 1.3$`\mathrm{}`$. The Crab pulsar is evident as the slightly dimmer upper component of the central double star.

As detailed in Golden et al. (2000), phase resolved images were obtained by firstly extracting all photons with a fixed pixel radius of the Crab pulsar, barycentering them, and then geometrically aligning each dataset in \[x,y\] per colour band. This yielded a sequence of time-series per colour band that, when read by a modified epoch-folding algorithm, generated both standard light curves and phase-resolved 2-dimensional image within a certain specified phase range. The latter used the appropriate data taken from the Jodrell Bank Crab Ephemeris (Lyne & Pritchard, 1996).

## 5 Analysis - Temporal Issues

It is absolutely crucial in this analysis to ascertain the temporal quality of our datasets, and in this way, test for possible glitches in the calibrated timing hardware and potential nonlinearities caused as a result of microchannel plate saturation, thus assessing any unusual temporal activity in the pulsar time series.

Regarding the problem of $`MAMA`$ saturation it was determined that, within the bounds of variable sky transparancy, no non-Poisson behaviour at the 99% level was noted, for a range of datasets over the $`UBV`$ wavebands (Golden 1999). This was repeated throughout the light curve, indicating that the emission process is consistent with dominant incoherent synchrotron processes. Giant radio pulses, which are of order 103 times more intense than the average radio pulse, occur in approximately 1 % of radio pulse periods (Lundgren et al. 1994). There were no deviations in the photon statistics accumulated in these observations consistent with such giant pulse activity.

Regarding the temporal stability of the timing system, we previously determined the timing accuracy to within $``$ 1 $`\mu `$s of UTC, and for long (4-5 days) timeseries datasets, it is possible to test the timing against the Crab pulsar itself via epoch folding or Fourier methods. For the former case, this will be as accurate as the timing ephemeris provided - during our observing period the Jodrell Ephemeris was accurate to 50 $`\mu `$s. In the latter, we are restricted to the accuracy with which one can determine to have detected a harmonic at high $`\nu `$. Via epoch folding, the best reduced $`\chi ^2`$ was determined for the given rotational parameters from Jodrell within 100 $`\mu `$s of the given reference time stamp given by the GPS system. Using an Origin2000 supercomputer, detection of the high harmonics (typically $``$ 50) for several daily $`B`$ datasets yielded frequency estimates that were accurate to the corrected ephemeris frequency to within 1 in 107. Consequently, we were satisfied with the temporal integrity of these MAMA datasets taken in Russia in January 1996.

In this vein, the unusual low frequency signal of $``$ 0.016 Hz reported by Cadez & Galicic (1996) was put to the test. FFT analysis on each night’s dataset binned at 1, 0.1 & 0.2 Hz were tested, but no significant signal was apparent within the 0.01 - 1 Hz window - or at greater frequencies. Consequently we must remain somewhat guarded as regards the conclusions of Cadez & Galicic (1996). The definitive discovery of a 60 second precessional signal from an isolated rotating neutron star would have profound implications for fundamental physics and thus demands conclusive and undeniable evidence for its justification. We must state that we see no such evidence to date.

## 6 Analyses - Phase Resolved Photometry

By the use of developed epoch folding algorithms, it was possible to take the deep summed time series for all three wavebands and re-construct a phase-resolved image for a specific phase region. Figure 2 shows such an ad hoc mosaic of phase-resolved images, with the phase resolution here 0.05 of phase.

It is apparent that there appears to be residual emission from the Crab pulsar during the traditional ‘off’ phase (essentially row 4 in Figure 2), as was originally noted by Peterson et al. (1978). The fact that we can obtain deep phase-resolved images which have been corrected for telescope drift and wobble with a sufficient field of view enables us to rigorously characterise the time-resolved nature of the pulsar’s emission via standard image processing techniques.

There are a number of ways to proceed with this analysis. Principally, we wish to determine 1) what the total background is in a given aperture per colour band and 2) what are the respective fluxes for phase-resolved components of emission associated with the Crab pulsar. The first step is perhaps the most important as it may then be applied to the standard derived light curves so as to correct for the background per bin and thus allow one to characterise the so-called ‘steady’ or ‘off’ phase of emission.

In Golden et al. (2000), the image processing techniques were outlined in some considerable detail for this analysis. It is sufficient to state that the IRAF daophot package was used to remove the full Crab image (and thus obtain the background per aperture per colour band in addition to the full Crab flux), and similarly for the various phase-resolved light curve components, which are localised in terms of phase in Table 1. This process involved the empirical fitting of a PSF to the Crab and making the subsequent extraction.

The ‘off’ phase here was localised in the following manner. Correcting for the estimated background above, the ‘off’ regions for the various coloured light curves at 11 $`\mu `$s resolution (3000 bins) were tested via a $`\chi ^2`$ algorithm. Here, the bin region was iteratively reduced, and tested to see if the observed flux was consistent with that expected from a non-varying source over that phase region, at the 95 % level.

As regards any nebular contribution to this photometric analysis, no significant change was noted between the estimated flux of the ‘off phase’ made via extraction using the full cycle PSF, and the flux determined by empirically fitting a PSF to the ‘off phase’ itself. This indicates that the local nebulosity does not significantly contribute to the determined fluxes in this work.

### 6.1 Phase Averaged/Integrated Colour Photometry

Actual phase-resolved photometry was done by normalizing the total Crab flux per colour band to a known integrated photometric reference, and then determining the contribution per phase region based on the fractional fluxes estimated for that particular phase-resolved component. The extinction corrected $`UBV`$ ground-based photometry of Percival et al. (1993) was used as the reference data for this purpose

Table 2 displays the tabulated fractional fluxes estimated following the photometric analysis detailed previously within each colour band - the fluxes were derived via renormalizing the fractional flux with respect to the Percival et al.’s integrated data.

In Table 3 the estimated power-law parameter $`\alpha `$ has been determined via a weighted linear least-squares fit for each individual spectral dataset, with associated errors. We have re-calculated $`\alpha `$ for the both the full range ($`UVUBV`$) & $`UBV`$ ground-based Percival et al. (1993) datasets to compare with the other power-law fits. Included in this Table are power-law fits to the FWHM and HWHM components.

### 6.2 Phase Resolved Colour Spectra - the Peaks & Bridge

Figure 3 shows as a function of photon energy the spectral forms based on the fluxes of each peak, and to the same scale and the emission from the bridge component.

It seems clear that certainly as regards the bulk flux from both peaks, that the functional forms, with $`\alpha `$ $``$ 0.07 $`\pm `$ 0.19 for peak 1 and $`\alpha `$ $``$ -0.06 $`\pm `$ 0.19 for peak 2, are to first order identical with that reported for the integrated emission, namely a flat power law ($`\alpha `$ $``$ -0.07 for the ground-based $`UBV`$, $`\alpha `$ $``$ 0.11 for the full ground-based $`UVUBV`$ Percival et al. (1993) dataset). This is not entirely surprising, as the dominant emission from the pulsar is localised in the peaks, and the fact that both are identical in form suggests a common population of electrons/magnetic field environment/Lorentz factor. Figure 4 substantiates this common source, there being no significant difference between the colour bands evident in these datasets. There appears to be, however, evidence for a somewhat steeper power-law associated with the bridge component of emission ($`\alpha `$ $``$ -0.44 $`\pm `$ 0.19).

In Figure 5 the fractional bridge/peak fluxes are shown as a function of photon energy, both suggesting similar functional forms, consolidating the common proportionality between each flux component. This suggests to first order a different population of emitting electrons / magnetic field environment / Lorentz factor in terms of the bridge component than that for the peaks. However it should be stressed that any definitive differing power-law forms are only significant at the 2$`\sigma `$ level - in this we are at the mercy of the reference data of Percival et al. (1993), where the substantial error component originates. Using other reference data of a higher photometric quality and further observations at other frequency ranges (such as in the $`R`$) would constrain these fits significantly.

### 6.3 Phase Resolved Colour Spectra - Leading & Trailing Peak Edges

In Figures 6 & 7 the HWHM flux estimates for both peaks 1 and 2 are shown - here the leading and trailing HWHM fluxes are derived for each peak. It seems clear that the leading & trailing edges for both peaks present spectral forms that are in agreement with each other (see Table 3), although we are somewhat restricted in our analysis by error considerations.

For both cases, we see suggestions for power-laws with decreasing exponents as a function of frequency, although both are consistent with an essentially flat functional form. These results are in agreement with estimates of Eikenberry & Fazio (1997), whose analysis yielded results which incline towards a decrease in flux with increasing photon energy. It is at $`JHK`$ wavebands that Eikenberry et al. (1996) observed quite profound changes in the overall spectral forms for the leading & trailing edges of both peaks - certainly in the $`UBV`$ regime, we appear to be sampling similar synchrotron emitting populations.

## 7 Analyses - Morphological Characteristics

As has been outlined previously, there has not been to date a consistent set of observational datasets that are temporally accurate (to 1 $`\mu `$s of UTC), and this has to some extent compromised previous efforts. The data taken in January 1996 provides the first such thorough sample, and we consequently first concentrate on its temporal analysis.

With an effective estimate of the background obtained for all three colour bands as previously outlined, one can proceed with a rigourous morphological analysis of the light curves themselves. This is essentially pursued by use of both least-squares fits to the Poisson ‘corrupted’ light curve profile, and the use of Monte Carlo techniques to ascertain the errors on any ideal topological fit to the data. In this way, we can hope to obtain information on a given peak’s precise maximum position in terms of phase, the extent of a plateau within the peak’s maximum region, and the peak’s FWHM and HWHM characteristics.

The basic algorithm fits third order polynomials to the leading and trailing edges of the given peak (from 30% to 90% maximum intensity), and then a sixth order polynomial across the peak cusp (from 90% to 90%). This latter fit is then used as a template in the following manner. A Monte Carlo routine uses the template as a basis from which to ‘add’ Poisson noise dependent on the template counts per specific bin. Having done this, a subsequent sixth order polynomial is fit via least-squares. From this, the phase bin with the maximum intensity is determined, and ‘noted’. This maximum phase is then used in the calculation of the FWHM and HWHM parameters, using the previously determined fits to the leading and trailing edges of a given peak.

An iterative routine then starts with the two adjacent bins to the nominated maximum bin, computes the local average per bin, and computes the $`\chi ^2`$ to test for a deviation from this local average i.e. the end of a ‘plateau’ region. The range is extended by a bin on each side, and the process is continued until significant deviations are noted at the 95% and 99% confidence intervals.

Finally, within the defined bin regions defined by these confidence ranges, a final routine starts with 2 bins, and sweeps through the ranges, noting the $`\chi ^2`$ statistic, and incrementing the bin number. The $`smallest`$ number of bin sizes that deviate from the required confidence levels are noted. The entire process is repeated from the initially derived template $`10^4`$ times. For each chosen parameter, the average is calculated, and the standard error on the $`10^4`$ samples is used to represent 1$`\sigma `$ errors on that mean.

The above outlined approach has its origin in the analytical techniques of Ransom, Eikenberry and others (Ransom et al., (1994), Eikenberry et al. (1996)) with the Monte Carlo derived estimate of a plateau region incorporated. The completed temporal analysis for the light curves in $`UBV`$ is tabulated in Table 4.

### 7.1 Peak 1 position & Peak-to-Peak Phase Difference

For each of the two peaks per colourband, as stated above, the iterative Monte Carlo routine determines the positions of the ‘best fitted’ peak maxima, and collates the difference between the two in phase, as well as recording the position in phase of the first peak. Figure 8 displays the peak-to-peak difference as a function of frequency. It would seem that there is no clear trend for the peaks phase difference to deviate within this restricted photon energy regime.

Based on the errors associated with our data, we do not see any departure from this trend reported elsewhere in the literature that the peak-to-peak phase difference is seen to contract with increasing energy (e.g. Eikenberry et al. 1997; Ramanamurthy 1994). Similar error concerns restrict any real attempt to definitively state that the $`UBV`$ main peaks arrive at differing phase intervals. The trend, if any, is that the radio-optical main peak difference in arrival phase increases with photon energy - broadly consistent with Percival et al. (1993) and other authors.

### 7.2 Full Width Half Maximum Differences

The FWHM of a pulsar’s pulsed components as a function of energy have been used for many years as a way of understanding the emission characteristics of these objects, being undoubtedly related to the geometry of the actual emission mechanism. The expected FWHM behaviour with photon energy is thus a strong function of the model under scrutiny - for example, the standard Outer Gap formalism predicts a growing FWHM with decreasing photon energy (Eikenberry et al. 1997), and early rather monolithic models envisaging emission occuring from an open cone interpret the morphologies as the collected open ‘cones’ of synchtrotron emitting population of electrons. In this latter case, a similar effect is expected. Our Monte Carlo results, as tabulated and presented in Figure 9

for both peaks substantiate this trend - although there seems to be a statistically significant greater gradiant associated with peak 2 compared with peak 1. From the previous work of Eikenberry et al. (1997) using the HSP dataset, it is impossible to make such a differentiation.

### 7.3 Half Width Half Maximum Differences

In their analysis, Eikenberry et al. (1997) found statistically significant differences between the leading and trailing half-widths for both peaks, these deviations clearly evident at the lowest photon frequencies ($`JHK`$), with the trend at our optical wavebands being an increase in measured phase extent and decreasing photon energy. Figures 10 and 11 display the results following the Monte Carlo analysis (tabulated in Table 4), and there is agreement with these earlier conclusions within the spread of errors.

The leading edges for both peaks display similar functional forms.

### 7.4 Plateau on the Main Peak

Resolution of the main peak has been discussed many times in the literature as a diagnostic as regards the basic emission mechanism for principally the Crab pulsar, as data with the required S/N is unlikely (at this stage) to come from any other object; but as outlined earlier, crude temporal resolution and not sufficiently rigorous analytical techniques have always dogged this. With the background-subtracted datasets, the same Monte Carlo code was used to ascertain the extent of non-deviating ‘best-fit’ 6th order polynomial across the cusp of the main peak with respect to the local noise (at the 95% and 99% confidence levels). The results of this approach are tabulated in Table 5, and in Figure 12 we show the estimated plateau duration as a function of photon energy.

Within the errors, we can say to first order that the cusp of the main peak is consistent with the presence of a ‘plateau’ of emission lasting some $``$ 50 $`\mu `$s in exent, which is in agreement with the previous conclusions of Komarova et al. (1996). Tantalisingly, a possible trend of decreasing plateau exent with increasing photon energy will have to await further datasets; there do appear to be indications that this may be the case from our data.

## 8 Discussion

Putting it succinctly, despite the excellent temporal resolution of our data, our analysis and subsequent conclusions are limited by both the low throughput of the imaging system and the large errors associated with the reference integrated photometry of the Crab pulsar. Nevertheless, we have resolved to an unprecendented photometric accuracy the flux fractions associated with the principal regions of the light curve, most especially that of the unpulsed component. Within errors, the datasets confirm existing trends seen in previous analyses; the intra-peak separation, peak FWHM/HWHM, and bridge flux as a function of photon energy. In this, we see nothing unusual to most model predictions. However the following results provoke further scrutiny.

### 8.1 Resolution of the Unpulsed Component of Emission

We have outlined in some considerable detail elsewhere the implications of the observed spectral form for the resolved unpulsed flux component (Golden et al. 2000). Clearly it is nonthermal, and spectrally similar to that of the bridge of emission. In our previous communication we argued that the observed unpulsed component of emission has its source in a similar electron population/magnetic field/Lorentz factor environment, but stated that an unequivocal link was not evident, and that viewing geometry was likely to play a substantial role in what we ultimately observe.

### 8.2 The temporal extent of the Main peak’s cusp

Resolution of the main peak’s cusp temporally has long been regarded as a possibly powerful diagnostic as regards determining the local physical conditions where the emission is occuring, which in the optical regime is believed to be as a result of incoherent synchrotron processes. In this analysis we have obtained rigorous estimates of the extent of this plateau, namely $``$ 55 $`\mu `$s with some indication of an increase with wavelength. It now remains to place this result in some sort of theoretical context. The models imply emission either occuring from a spatially extended region (the Outer Gap of Romani & Yadigaroglu (1995) or from a localised cone/‘trumpet’ like zone traced by the opening angle above a polar cap, and subsequently modified by the toroidal dipolar field. Thus for the former, the peaks are a consequence of many photons arriving from disparate locations but in phase due to relativistic effects, and in the latter we have a conceptually more appealing core/cone of emission. In both cases, the synchrotron emitting electron’s are expected to have a common Lorentz/local magnetic field, in order to satisfy the observed luminosity, although the intuitive concept of a ‘beam’ is problematic for the former $`ansatz`$.

For a temporal extent of 55 $`\mu `$s at the main peak, one could argue that this implies that the dominant beaming from the population of synchrotron emitting electrons from some localised area in the magnetosphere has a total angular scale of $``$ 10-1.98 radians. It follows (Sturrock et al. 1975) then that if one envisages that the total optical emission comes from these electrons with a Lorentz factor $`\gamma `$ and pitch angle $`\xi `$, then

$$\gamma \gamma _o=10^{1.98},\xi \xi _o=10^{1.98}$$
(1)
suggesting a $`low`$ Lorentz factor emission regime ($``$ $`10^2`$). If the radiation is via normal (large-angle) synchrotron emission, then the emitted spectrum peaks at the frequency

$$\nu =\frac{3}{2}\nu _G\gamma ^2\xi \mathrm{Hz}(\gamma \xi \frac{4}{3})$$
(2)
where $`\nu _G`$ is the gyrofrequency given by

$$\nu _G=10^{6.5}B\mathrm{Hz}$$
(3)
In the small angle approximation, the spectrum peaks at

$$\nu =2\nu _G\gamma \mathrm{Hz}(\gamma \xi \frac{4}{3})$$
(4)
It is easily shown that for either case,

$$\nu >2\nu _G\gamma _o\mathrm{Hz}$$
(5)
From Eikenberry et al. (1997), the optical spectrum appears to show a maximum in the vicinity of the H band ($``$ 1.8 $`\times 10^{14}`$ Hz) or $`10^{14.26}`$. Applying this with the above suggests that the local magnetic field must be

$$B3.02\times 10^5\mathrm{Gauss}$$
(6)
Assuming a dipolar $``$ $`r^3`$ law, and the canonical 10km neutron star possessing a derived surface magnetic field of 3.7 $`\times `$ $`10^{12}`$ Gauss, this corresponds to a radial distance of some 2256 km. In comparsion, the $`R_{lc}`$ is estimated to be $``$ 1523 km - thus emission must occur at some 1.5$`R_{lc}`$. This emission is not possible for the ideally orthogonally aligned rotator of Smith (1978) and others, but feasible for those pulsars whose angle $`\alpha `$ between dipole & rotational axes is $`<`$ 90o. It is clear that the assumption of restricted radial zones of emission would favour the more localised models over the spatially extended (where zonally varying local $`\gamma `$ & $`B`$ would contribute to the overall pulse profile) and that this radial estimate has many intuitive links with the emission model of Gil et al. (2000).

As a final comment, we note some evidence for the plateau’s temporal extent to grow with wavelength, errors notwithstanding. This trend, if real, could be possibly explained most easily with reference to synchrotron self-absorption. Previous discussions on this topic assumed that the observed roll-over was due to this phenomenon, yet time-resolved observations of the $`IR`$ wavebands by several authors (Middleditch et al. 1983; Eikenberry et al. 1996, 1997) do not neccesarily show evidence for a flattening of the main peak (however, Pennypacker 1981; Penny 1982). This effect should be manifested in the $`IR`$ wavebands, and it is somewhat difficult to justify at our low frequencies for a tightly localised synchrotron source - but it is possible in a radially or latitudinally extended model, where differing regions of the magnetosphere may experience the effect. Alternatively, we can define the cusp as a simple function of the local synchrotron emission processes. Here, the lower energy electrons have a wider opening angle, and thus we observe a wider plateau. Accurate numerical modelling would test this latter hypothesis.

### 8.3 Phase difference between radio & $`UBV`$ Main peak arrival

From the tabulated data (Table 4), we can say for the higher S/N data of $`B`$ & $`V`$ that the main peak pulses arrive to within 10 $`\mu `$s of each other, and that, when one considers the peak-to-peak separation, there does not seem to be any clear evidence for variation as a function of photon energy. In this we are somewhat restricted by the overall poor bin S/N for these datasets. Deeper observations would provide temporally more accurate data, but to first order, it does seem that in these wavebands, we are seeing an essentially constant general behaviour, as regards pulse arrivals. We cannot absolutely differentiate the radio/pulse lag time, due to the error on the former (some 50 $`\mu `$s), but we can place a relative lag of order $``$ -60 $`\mu `$s - that is to say, the $`UBV`$ emission trails the arrival of the radio main pulse. This might not be so difficult to understand if the radio emission had its source at some distance further in towards the neutron star surface with respect to the higher energy sources, located in the toroidally warped dipolar field, where both relativistic beaming and the local magnetic field vector might conspire to produce such a delayed arrival time (assuming localised emission). But this is not the case for the Crab, whose radio emission is believed to occur in proximity to the optical sources close to the light cylinder. Observations with the FIGARO II $`\gamma `$-ray instrument (Masnou et al. 1994) have suggested that the (0.15 - 4.0 MeV) $`\gamma `$-ray main peak $`precedes`$ that of the radio pulse by some 375 $`\pm `$ 148 $`\mu `$s, which according to the authors, points to spatially different regions of emission between the radio & $`\gamma `$-ray source regions of some 100 km in the magnetosphere. Recent RXTE observations (Rots et al. 1998) at a superior temporal resolution with respect to UTC suggest that the (5-200 keV) X-ray main peak leads the radio pulse by 264 $`\pm `$ 330 $`\mu `$s.

For our ‘best’ position, the $`V`$ band, the main peak is localised at +79.2 $`\pm `$ 10 $`\mu `$s with respect to the radio pulse - incorporating the quoted radio error enhances this to $`\pm `$ 60 $`\mu `$s. Thus to first order, the optical trails the X-ray main peak, and the X-ray trails the $`\gamma `$-ray main peaks, but their absolute location with respect to the radio is not clear, being dependent on the local epoch timing solution. The discrepancies are undoubtedly as a result of the differing quality of the different datasets (in terms of S/N), the accuracy of the radio ephemeris solution and the manner in which the main peak’s centroid is defined. From our folded light curves, it appears that the main peak in the $`UBV`$ is located $`after`$ the radio main pulse, but it is important to appreciate that there is a substantial error on the radio timing solution. What is required are simultaneous deep observations in both the optical & radio so as to guarantee as accurate a timing solution as possible. Such observations have recently been made in collaboration with colleagues from the Westerbork Radio Telescope Observatory using the TRIFFID instrument at La Palma; the data is currently under analysis. Realistic numerical modelling is possibly the only full consistent manner in which to reproduce such frequency dependent lag effects and so constrain the position(s) of emission for a standard dipolar structure under a range of viewing geometries, for both spatially extended & localised model frameworks. Such modelling studies are underway at NUI, Galway and elsewhere.

## 9 Conclusion

Definitive optical observations of the Crab and other pulsars have to date been restricted - either by poor temporal resolution both in terms of resolution and with respect to absolute UTC, and also by the limited spatial information which is crucial in distinguishing the pulsar’s emission from the background. Using the TRIFFID high speed photometer, we have obtained datasets uncompromised by these previous limitations, and in addition to confirming many of the previous conclusions of Eikenberry et al. (1997), we have determined several new results uniquely due to the technology we have implemented, including the fluxing of the unpulsed component of emission, the absolute arrival times of the three light curves with respect to one another, and the spatial extent of the plateau on the main peak. Our photometric analysis was however limited by the reference data we used to renormalize the flux ratios thus determined. It is clear that further analysis using more photometrically accurate reference data combined with further observations at extended wavebands, such as the $`UV`$ or $`R/I`$ would yield fluxes that would provide excellent leverage to the least-squares fits and thus produce more significant power-law exponents. In terms of a more rigorous temporal analysis of the light curves, deeper observations in $`UBV`$ would place limits on the extent of the main peak’s plateau, and possibly indicate if it is energy dependent. Furthermore, a greater sample of photons per pulse period may allow us to probe the possible correlation between the radio giant pulses and high energy emission mechanisms. We note that the incorporation of polarimetric technology into this nascent field of 2-d high speed photometry would yield critical data which would undoubtedly reveal more clues to puzzle over, and thus contribute towards the eventual solution to this three decade old mystery of how exactly pulsars function.

###### Acknowledgements.

ESO is thanked for the provision of their MAMA detector. We are grateful to the 6-m telescope Program Committee of the RAS for observing time allocation. We thank the engineers of SAO RAS, A. Maksimov for help in equipment preparation for the observations and the Director of SAO RAS Yu. Balega for arranging the observations. Dr Ray Butler of NUI, Galway is thanked for his assistance during the course of the photometric analysis. This work was supported by the Russian Foundation of Fundamental Research (98-02-17570), State programme ”Astronomy”, Russian Ministry of Science and Technical Politics, and the Science-Educational Centre ”Cosmion”, and funded under the INTAS programme. The support of Enterpise Ireland, the Irish Research and Development agency, is very gratefully acknowledged. Finally, we would like to thank the anonymous referee’s considerable effort and wise counsel in contributing towards the final version.


format: text

text:
# Clopper-Pearson Bounds from HEP Data CutsThis research was partially funded by the Department of Energy under contract DE-FG02-97ER41022.

## 1 Introduction

The general theory of confidence bounds (or fiducial intervals) was developed by Fisher , Neyman and Pearson . We consider a particular problem which is of interests when cuts are used to analyze high energy physics data. Typically, a neural network or some other method of performing the cuts results in probabilities (efficiencies) to tag signals more likely than background events. For instance by means of Monte Carlo (MC) simulations, these probabilities can normally be calculated. Let $`P_s`$ be the probability to tag a signal and $`P_b`$ be the likelihood to tag a background event, $`0<P_b<P_s<1`$. Out of a total number of $`N`$ data one gets in this way

$$N^Y\mathrm{tagged}\mathrm{data}.$$
(1)
It is easy to find from this the mean expectation for the signal probability. Assume that there are $`N_s`$ signals and $`N_b`$ background events in the data. Then we have

$$N=N_s+N_b\mathrm{and}N^Y=P_sN_s+P_bN_b$$

and these two equations solve for

$$p_{\mathrm{mean}}=\frac{N_s}{N}=\frac{N^YP_bN}{N(P_sP_b)}.$$
(2)
The question is, what are the implied confidence limits on the signal probability?

The special case $`P_s=1`$ and $`P_b=0`$ (sure signal detection) has been treated by Clopper and Pearson  in 1934. After briefly reviewing their approach in the next section, we derive and illustrate the general case in section 3. This is, in part, based on Ref.. In particular, the method is valid when the number of tagged data is small and returns the probability $`P_0`$ for the case that there is no signal, i.e. that the exact signal probability is $`p_{\mathrm{exact}}=0`$. This is of interest for the statistical analysis of the ongoing Higgs search at LEP . Discovery of the Higgs particle on the $`5\sigma `$ level would mean $`P_0\mathrm{0.287\hspace{0.17em}10}^6`$. Subsection 3.3 explains the use of the corresponding Fortran programs and how to download them from the web. Conclusions follow in the final section 4.

## 2 The Clopper–Pearson Confidence Limits

Let $`p`$ be the likelihood that a data point is a signal. For $`N`$ measurements the probability to observe $`k`$ signals is given by the binomial coefficient

$$b_N(k,p)=(\begin{array}{c}N\\ k\end{array})p^kq^{Nk}=\frac{N!}{k!(Nk)!}p^kq^{Nk}\mathrm{with}q=p1.$$
(3)
The probability to observe $`kN_s`$ signals is given by

$$P_{kN_s}(p)=\underset{k=N_S}{\overset{N}{}}b_N(k,p)$$
(4)
and the probability to observe $`kN_s`$ signals by

$$P_{kN_s}(p)=\underset{k=0}{\overset{N}{}}b_N(k,p).$$
(5)
For $`N=26`$ and $`N_s=10`$ the functional forms of $`P_{kN_s}(p)`$ and $`P_{kN_s}(p)`$ are depicted in figure 1.

Assume that $`N_s`$ signals are found in $`N`$ measurements and that a probability $`Q^c<0.5`$ (typical values are $`Q^c=0.16`$ or $`Q^c=0.025`$) is given. We can solve equation (4) for $`P_{kN_s}(p_{})=Q^c`$ and $`p_{}`$ is a lower bound on the true signal probability $`p_{\mathrm{exact}}`$, such that the likelihood to find $`kN_s`$ signals in $`N`$ measurement is smaller than $`Q^c`$ for every $`p_{\mathrm{exact}}<p_{}`$. Correspondingly, we can solve equation (5) for $`P_{kN_s}(p_+)=Q^c`$ and $`p_+`$ is an upper bound on the true signal probability $`p_{\mathrm{exact}}`$, such that the likelihood to find $`kN_s`$ signals in $`N`$ measurement is smaller than $`Q^c`$ for every $`p_{\mathrm{exact}}>p_+`$. Together, this combines into the Clopper–Pearson bounds: The probability to find the true signal probability in the range

$$p_{}p_{\mathrm{exact}}p_+\mathrm{is}\mathrm{larger}\mathrm{or}\mathrm{equal}\mathrm{to}P^c=12Q^c.$$
(6)
In more details the meaning of the inequality is discussed in . For $`P^c=0.68`$ ($`Q^c=0.16`$) the $`p^\pm `$ values are indicated in figure 1. Approximately, this range corresponds to the confidence of a $`1\sigma `$ error bar. Similarly the confidence range corresponding to a $`2\sigma `$ error bar, etc., can be found.

## 3 Confidence Limits from Data Cuts

We are interested in the situation where signal and background data can no longer be distinguished unambiguously. Instead, a neural network or other device yields statistical information by tagging signals with efficiency $`P_s`$ and background data with probability $`P_b`$, as discussed in the introduction.

Applying the cuts to all $`N`$ data results in $`N^Y`$ tagged data ($`0N^YN`$), composed of $`N^Y=N_s^Y+N_b^Y`$, where $`N_s^Y`$ is the number of tagged signals and $`N_b^Y`$ is the number of tagged background data. Of course, the values for $`N_s^Y`$ and $`N_b^Y`$ are not known. Our task is to determine confidence limits for the signal probability $`p`$ from the sole knowledge of $`N^Y`$. We proceed by writing down the probability density of $`N^Y`$ for given $`p`$ and, subsequently, generalizing the Clopper-Pearson method.

First, assume fixed $`N_s`$. The probability densities of $`N_s^Y`$ and $`N_b^Y`$ are binomial and thus the probability density for $`N^Y`$ is given by the convolution

$$P^Y(N^Y|N_s)=\underset{N_s^Y+N_b^Y=N^Y}{}b_{N_s}(N_s^Y,P_s)b_{N_b}(N_b^Y,P_b),N_b=NN_s.$$
(7)
Proof: For a signal event the probability to be tagged is $`P_s`$, so $`b_{N_s}(N_s^Y,P_s)`$ is the probability to tag $`N_s^Y`$ out of the $`N_s`$ signals. Similarly, the probability for a background event to become tagged is $`P_b`$ and $`b_{N_b}(N_b^Y,P_b)`$ is the probability to tag $`N_b^Y`$ of the $`N_b=NN_s`$ background events. As these two probabilities are independent, the likelihood that precisely $`N_s^Y`$ of the signals and $`N_b^Y`$ of the background events are tagged becomes the product $`b_{N_s}(N_s^Y,P_s)b_{N_b}(N_b^Y,P_b)`$. Summing over all possibilities which add up to $`N_s^Y+N_b^Y=N^Y`$ gives the result.

Summing over $`N_s`$ in (7) removes the constraint of fixed $`N_s`$ and, with $`N`$, $`p`$ fixed, the probability to tag $`k`$ events becomes

$$b_N^Y(k,p)=\underset{N_s=0}{\overset{N}{}}b_N(N_s,p)P^Y(k|N_s).$$
(8)
Fourier transformation of the convolution (7) allows for an efficient numerical calculation of the $`P(k|N_s)`$ coefficients. In analogy with equations (4) and (5) we find the probabilities to tag $`kN^Y`$ and $`kN^Y`$ events to be

$$P_{kN^Y}^Y(p)=\underset{k=N^Y}{\overset{N}{}}b_N^Y(k,p)\mathrm{and}P_{kN^Y}^Y(p)=\underset{k=0}{\overset{N^Y}{}}b_N^Y(k,p).$$
(9)
For $`N=35`$, $`N^Y=12`$, $`P_s=0.8`$ and $`P_b=0.05`$ the functions  $`P_{kN^Y}^Y(p)`$ and $`P_{kN^Y}^Y(p)`$ are depicted in figure 2. The 68% confidence range (6) is also indicated in figure 2, where the bound values $`p_\pm `$ are now defined as solutions of the equations

$$Q^c=P_{kN^Y}^Y(p_{})\mathrm{and}Q^c=P_{kN^Y}^Y(p_+).$$
(10)
The range $`[p_{},p_+]`$, obtained with $`Q^c=0.16`$, guarantees the standard one error bar confidence probability of 68% for every true signal probability $`p_{\mathrm{exact}}`$. For almost all values the actual confidence will be better. However, the bounds cannot be improved without violating the requested confidence probability for the case that $`p_{\mathrm{exact}}`$ happens to agree with either $`p_{}`$ or $`p_+`$. In the same way, bounds calculated with $`Q^c=0.023`$ ensure the standard two error bar confidence level of 95.4% or better, and so on.

### 3.1 Data Sets with Few Signals

As outlined, data sets with few signals are of of particular interest in high energy physics. Let us replace $`N^Y=12`$ of figure 2 by $`N^Y=3`$. The resulting graph is depicted in figure 3. From the $`P_{kN^Y}^Y(p)`$ curve we read off the finite probability $`P_0=0.254`$ for the likelihood that the true signal probability $`p_{\mathrm{exact}}=0`$ generates $`kN_y`$ tags. Due to this probability the lower 68% confidence bound $`p_{}`$ disappears, whereas the upper $`p_+`$ bound does still exist. In passing let us note that for the data of figure 2 we have $`P_0=\mathrm{0.69\hspace{0.17em}10}^7`$, i.e. there $`p_{\mathrm{exact}}=0`$ is ruled out on the $`5\sigma `$ level.

### 3.2 Signal Probability Distributions

To avoid frequentist objections, the cumulative signal distribution function $`F(p)`$ is, in the opinion of the author, best defined as the expectation of the researcher to find the true signal probability $`p_{\mathrm{exact}}`$ in the range $`0p_{\mathrm{exact}}p`$. Our approach allows to estimate upper and lower bounds for the cumulative signal distribution function

$$F(p)=_{\mathrm{}}^p𝑑p^{}f(p^{})\mathrm{where}f(p^{})\mathrm{is}\mathrm{the}\mathrm{signal}\mathrm{probability}\mathrm{density}.$$
(11)
Equation (10) implies

$$F_1(p)=1P_{kN^Y}^Y(p)=\underset{k=N^Y+1}{\overset{N}{}}b_N^Y(k,p)F(p)F_2(p)=P_{kN^Y}^Y(p)=\underset{k=N^Y}{\overset{N}{}}b^Y(k,p)$$
(12)
and note that $`F(0)=F_2(0)`$. Figure 4 shows the result for the same data which were used in figure 2. Any reasonable Bayesian estimates (which involves additional a-priori assumptions) should give a function $`F(p)`$ which is sandwiched between $`F_1(p)`$ and $`F_2(p)`$.

It is instructive to define peaked distribution functions  by

$$F_i^{\mathrm{peaked}}(p)=\{\begin{array}{cc}F_i(p)\mathrm{for}F_i(p)0.5\mathrm{and}\hfill & \\ 1F_i(p)\mathrm{for}F_i(p)0.5,(i=1,2).\hfill & \end{array}$$
(13)
Using the same data as in figure 4, $`F_1^{\mathrm{peaked}}(p)`$ and $`F_2^{\mathrm{peaked}}(p)`$ are depicted in figure 5. The advantages of using peaked distribution functions instead of conventional cumulative distribution functions are:

1. The ordinate becomes enlarged by a factor of two.
2. The estimated medians are located at the peaks and the probability content of the distribution is instructively displayed.
3. Bounds like those of figure 1 are easily read off.

The probability densities corresponding to the cumulative distribution functions (12) are the derivatives of the $`F_i(p)`$ with respect to $`p`$

$$f_i(p)=\frac{dF_i(p)}{dp}(i=1,2).$$
(14)
Their numerical calculation is straightforward when analytical expressions for the derivatives of the binomial coefficients in equation (8) are used. Figure 6 exhibits the results for $`f_1(p)`$ and $`f_2(p)`$ corresponding to $`F_1(p)`$ and $`F_2(p)`$ of figure 4. At $`p=0`$ the probability densities have $`\delta `$-function contributions

$$f_i(p)=F_i(0)\delta (p)+\mathrm{},(i=1,2).$$
(15)
In case of figure 6 the $`F_i(0)`$ coefficients are practically zero. However, for the probability densities corresponding to figure 3 there would be a substantial contribution: $`F_1(0)=0.096`$ and $`F_2(0)=0.254`$ in that case.

### 3.3 The Fortran Code

The Fortran code which produces the illustrations of this paper is available on the web. Start at the author’s homepage  www.hep.fsu.edu/~berg and follow the research and from there the Clopper-Pearson hyperlink. Load down all files of the Fortran Programs subdirectory into an empty directory. Any standard Fortran 77 compiler should then be able to compile the cp1.f, cp2.f and cp3.f programs. Be aware that the program files include (via include Fortran statements) some of the other files you downloaded. Running one of the programs produces a data files with the name of that program and a .d extension. Subsequently, gnuplot users can produce the graphical presentations of this paper by using the \*.plt driver files, as listed in the following.

$`\mathrm{𝚙𝚛𝚘𝚐𝚛𝚊𝚖}`$ $`\mathrm{𝚐𝚎𝚗𝚎𝚛𝚊𝚝𝚎𝚜}\mathrm{𝚏𝚒𝚕𝚎}`$ $`\mathrm{𝚞𝚜𝚎}\mathrm{𝚠𝚒𝚝𝚑}\mathrm{𝚐𝚗𝚞𝚙𝚕𝚘𝚝}\mathrm{𝚍𝚛𝚒𝚟𝚎𝚛}`$
$`\mathrm{𝚌𝚙𝟷}.𝚏`$ $`\mathrm{𝚌𝚙𝟷}.𝚍`$ $`\mathrm{𝚌𝚙𝟷}.\mathrm{𝚙𝚕𝚝}`$
$`\mathrm{𝚌𝚙𝟸}.𝚏`$ $`\mathrm{𝚌𝚙𝟸}.𝚍`$ $`\mathrm{𝚌𝚙𝟸}.\mathrm{𝚙𝚕𝚝},\mathrm{𝚌𝚙𝟺}.\mathrm{𝚙𝚕𝚝},\mathrm{𝚌𝚙𝟻}.\mathrm{𝚙𝚕𝚝},\mathrm{𝚌𝚙𝟼}.\mathrm{𝚙𝚕𝚝}`$
$`\mathrm{𝚌𝚙𝟹}.𝚏`$ $`\mathrm{𝚌𝚙𝟹}.𝚍`$ $`\mathrm{𝚌𝚙𝟹}.\mathrm{𝚙𝚕𝚝}`$
Here, the gnuplot file number corresponds to the figure number of this paper. To get encapsulated postscript files, the comment signs in front of the first two rows of each gnuplot file have to be eliminated.

## 4 Conclusions

We have calculated confidence limits, and corresponding limits of the entire cumulative distribution function, for an unknown true signal likelihood $`p_{\mathrm{exact}}`$. The only input used are the efficiencies $`P_s`$ for tagging signals, the probabilities $`P_b`$ for tagging background events, the number $`N^Y`$ of tagged data and the total number of data $`N`$. In particular, the method allows to deal with the situation where only few signals occur and yields then a finite probability for the likelihood that $`kN^Y`$ tags are observed if the true signal probability is $`p_{\mathrm{exact}}=0`$. In real life the probabilities $`P_s`$ and $`P_b`$ are most likely estimators by themselves, i.e. quantities with error bars. This causes no major problem, one just has to apply our confidence calculations to an appropriate sample and to average over the results.


text:
# Blackbody radiation in a nonextensive scenario

## I Introduction

Planck’s blackbody radiation studies constitute one of the milestones in the history of Physics. Planck’s law is satisfactorily accounted by recourse of Bose-Einstein statistics. The possible existence of small deviations from this law in the cosmic blackbody radiation has recently been discussed. They would have arisen, at the time of the matter-radiation de-coupling  as a consequence of the fact that long-range interactions seem to be associated to a nonextensive scenario . The concomitant thermostatistical treatment  is by now recognized as a new paradigm for statistical mechanical considerations. It revolves around the concept of Tsallis’ information measure $`S_q`$ , a generalization of Shannon’s one, that depends upon a real index $`q`$ and becomes identical to Shannon’s measure for the particular value $`q=1`$.

The blackbody radiation problem was discussed within the Tsallis’ nonextensive framework in . It was there conjectured that the deviations from Planck’s law detected in the pertinent data  could be attributed to a long-range gravitational influence. The study reported in  employs the so-called Curado-Tsallis unnormalized expectation values  in the limit $`q1`$. Afterwards, Lenzi  advanced an exact treatment of the problem, also within the Curado-Tsallis framework.

Nowadays, it is believed that this framework has been superseded by the so-called normalized one, advanced in , which seems to exhibit important advantages . This normalized treatment, in turn, has been considerably improved by the so-called “Optimal Lagrange Multipliers” (OLM) approach . It is then natural to revisit the problem from such a new context.

In such a vein we shall consider blackbody radiation in equilibrium within an enclosure of volume $`V`$ with the goal of ascertaining the possible $`q`$-dependence of the Planck spectrum, Steffan-Boltzmann’s law and Wien’s one. The OLM treatment seems to be indicated in view of the findings of  regarding the particular nature of the Lagrange multiplier associated to the temperature.

The paper is organized as follows: In Sect. II we present a brief OLM primer, while in Sect. III we obtain the partition function, internal energy and energy density of our problem. In Sect. IV first-order corrections in $`1q`$ are discussed and, finally, some conclusions are drawn in Sect. V.

## II Main results of the OLM formalism

For a most general quantal treatment, in a basis-independent way, consideration is required of the density operator $`\widehat{\rho }`$ that maximizes Tsallis’ entropy

$$\frac{S_q}{k}=\frac{1Tr(\widehat{\rho }^q)}{q1},$$
(1)
subject to the $`M`$ generalized expectation values $`\widehat{O}_j`$, where $`\widehat{O}_j`$ $`(j=1,\mathrm{},M)`$ denote the relevant observables.

Tsallis’ normalized probability distribution , is obtained by following the well known MaxEnt route . Instead of effecting the variational treatment of Tsallis-Mendes-Plastino (TMP) , we pursue the alternative path developed in . One maximizes Tsallis’ generalized entropy (1)  subject to the constraints

$`Tr(\widehat{\rho })`$ $`=`$ $`1`$ (2)
$`Tr\left[\widehat{\rho }^q\left(\widehat{O}_j\widehat{O}_j_q\right)\right]`$ $`=`$ $`0,`$ (3)
whose generalized expectation values

$$\widehat{O}_j_q=\frac{Tr(\widehat{\rho }^q\widehat{O}_j)}{Tr(\widehat{\rho }^q)},$$
(4)
are (assumedly) a priori known. The resulting density operator reads

$$\widehat{\rho }=\overline{Z}_q^1\left[1(1q)\underset{j}{\overset{M}{}}\lambda _j\left(\widehat{O}_j\widehat{O}_j_q\right)\right]^{\frac{1}{1q}},$$
(5)
where $`\{\lambda _j\}`$ stands for the Optimal Lagrange Multipliers’ set and $`\overline{Z}_q`$ is the partition function

$$\overline{Z}_q=Tr\left[1(1q)\underset{j}{\overset{M}{}}\lambda _j\left(\widehat{O}_j\widehat{O}_j_q\right)\right]^{\frac{1}{1q}}.$$
(6)
It is shown in  that

$$Tr(\widehat{\rho }^q)=\overline{Z}_q^{1q},$$
(7)
and that Tsallis’ entropy can be cast as

$$S_q=k\mathrm{ln}_q\overline{Z}_q,$$
(8)
with $`\mathrm{ln}_q\overline{Z}_q=(1\overline{Z}_q^{1q})/(q1)`$. These results coincide with those of TMP  in their normalized treatment. If, following , we define now

$$\mathrm{ln}_qZ_q=\mathrm{ln}_q\overline{Z}_q\overline{Z}_q^{1q}\underset{j}{}\lambda _j\widehat{O}_j_q,$$
(9)
we are straightforwardly led to

$`{\displaystyle \frac{}{\widehat{O}_j_q}}\left({\displaystyle \frac{S_q}{k}}\right)`$ $`=`$ $`\overline{Z}_q^{1q}\lambda _j`$ (10)
$`{\displaystyle \frac{}{\lambda _j}}\left(\mathrm{ln}_qZ_q\right)`$ $`=`$ $`\overline{Z}_q^{1q}\widehat{O}_j_q.`$ (11)
Equations (10) and (11) are modified Information Theory relations like those used to build up an à la Jaynes , Statistical Mechanics. The basic Legendre-structure relations can be recovered in the limit $`q1`$.

## III Blackbody radiation

### A The standard situation $`q=1`$

In order to analyze the blackbody radiation in the generalized statistical context we first consider the standard situation, $`q=1`$. We look for the equilibrium properties of the blackbody electromagnetic radiation. The appropriate thermodynamical variables are the volume $`V`$ and the temperature $`T`$ .

The Hamiltonian of the electromagnetic field, in which there are $`n_{𝐤,ϵ}`$ photons of momentum $`𝐤`$ and polarization $`ϵ`$, is given by

$$\widehat{}=\underset{_{k,ϵ}}{}\mathrm{}\omega \widehat{n}_{k,ϵ},$$
(12)
where the frequency is $`\omega =c\left|𝐤\right|`$ ($`c`$ is light’s speed) and $`n_{k,ϵ}=0,1,2,\mathrm{}`$, with no restrictions on $`\left\{n_{𝐤,ϵ}\right\}.`$ The partition function reads

$$Z_1=Tr\left(e^{\beta \widehat{}}\right).$$
(13)
For a macroscopic volume we define a density of states in the $`n`$-dimensional space $`g_n(\omega )=A_n\omega ^{n1}`$ with

$$A_n=\frac{2\tau _nV}{(4\pi c^2)^{n/2}\mathrm{\Gamma }(n/2)},$$
(14)
where $`\tau _n=n1`$ is the number of linear-independent polarizations. Eq. (13) can be written as

$$Z_1=\mathrm{exp}\left\{_0^{\mathrm{}}𝑑\omega g_n(\omega )\mathrm{ln}\left[1\mathrm{exp}\left(\beta \mathrm{}\omega \right)\right]\right\}=e^{\xi _n},$$
(15)
where

$$\xi _n=\frac{I_nA_n}{(\mathrm{}\beta )^n},$$
(16)
with

$$I_n=_0^{\mathrm{}}𝑑xx^{n1}\mathrm{ln}\left(1e^x\right)=\mathrm{\Gamma }(n)\zeta (n+1),$$
(17)
$`\zeta `$ stands for the Riemann z-function and $`\mathrm{\Gamma }`$ is the Gamma function.

### B Generalized situation $`q1`$ ($`0<q<1`$)

#### 1 Generalized Partition Function

The OLM-Tsallis generalized partition function (see Section II) is

$$\overline{Z}_q=Tr\left[1(1q)\beta \left(\widehat{H}U_q\right)\right]^{\frac{1}{1q}},$$
(18)
where $`\widehat{H}`$ is the Hamiltonian (12), $`U_q`$ is the mean energy (4), and $`\beta `$ is the optimal multiplier Lagrange. The photons’ chemical potential is, of course, taken as zero.

With the aim of calculating $`\overline{Z}_q`$ we follow the steps of . We work again in the large volume limit and regard Eq. (18) as an integral that can be evaluated using the relation

$$\frac{e^{ab}b^{1z}}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{itb}}{(a+it)^z}=\{\begin{array}{cc}\frac{1}{\mathrm{\Gamma }(z)}& \hfill \text{ for }b>0\\ 0& \hfill \text{for }b<0,\end{array}$$
(19)
with $`a>0,`$ $`Re(z)>0`$, and $`\pi /2<\mathrm{arg}(a+it)<\pi /2.`$

If we set $`a=1`$, $`b=1(1q)\beta (\widehat{H}U_q)`$ (the cut-off condition is naturally fulfilled ) and $`z=1/(1q)+1`$, the generalized partition function adopts the appearance

$$\overline{Z}_q(\beta )=\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑tK_q(t)e^{\stackrel{~}{\beta }U_q}Z_1(\stackrel{~}{\beta }),$$
(20)
with $`Z_1`$ given by Eq. (15),

$$K_q(t)=\frac{\mathrm{\Gamma }\left[(2q)/(1q)\right]\mathrm{exp}(1+it)}{2\pi (1+it)^{(2q)/(1q)}},$$
(21)
and

$$\stackrel{~}{\beta }=(1+it)(1q)\beta .$$
(22)
In order to evaluate the integral in (20) we expand the exponential term and obtain

$$\overline{Z}_q(\beta )=\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\frac{\mathrm{\Gamma }[(2q)/(1q)]}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+(1q)\beta U_q)(1+it)}}{(1+it)^{\frac{2q}{1q}+nm}},$$
(23)
where $`\xi _n`$ is given by Eq. (16).

Using again (19) with $`b=1+(1q)\beta U_q,`$ $`a=1,`$ and $`z=(2q)/(1q)+nm`$ we arrive to

$$\overline{Z}_q(\beta )=\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\frac{\mathrm{\Gamma }[(2q)/(1q)]}{\mathrm{\Gamma }[(2q)/(1q)+nm]}\left[1+(1q)\beta U_q\right]^{\frac{1}{1q}+nm},$$
(24)
and notice that an additional cut-off’s like condition must be considered $`1+(1q)\beta U_q>0`$. Otherwise, $`\overline{Z}_q(\beta )=0`$.

#### 2 The internal energy

The generalized internal energy can be cast in the fashion

$$U_q=\frac{T_1}{T_2},$$
(25)
where we have introduced the definitions

$`T_1`$ $`=`$ $`Tr\left\{\left[1(1q)\beta (\widehat{H}U_q)\right]^{\frac{q}{1q}}\widehat{H}\right\},`$ (26)
$`T_2`$ $`=`$ $`Tr\left\{\left[1(1q)\beta (\widehat{H}U_q)\right]^{\frac{q}{1q}}\right\}.`$ (27)
We first analyze $`T_2`$. The expression $`\left[1(1q)\beta (\widehat{H}U_q)\right]^{\frac{q}{1q}}`$ is evaluated by recourse to (19), where we set $`a=1,`$ $`b=1(1q)\beta (\widehat{H}U_q),`$ and $`z=1/(1q).`$ We are thus led to

$$\left[1(1q)\beta (\widehat{H}U_q)\right]^{\frac{q}{1q}}=\frac{\mathrm{\Gamma }[1/(1q)]}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+it)[1(1q)\beta (\widehat{H}U_q)]}}{(1+it)^{\frac{1}{1q}}}.$$
(28)
The trace operation yields afterwards

$$T_2=\frac{\mathrm{\Gamma }[1/(1q)]}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+it)[1+(1q)\beta U_q]}}{(1+it)^{\frac{1}{1q}}}Z_1(\stackrel{~}{\beta }),$$
(29)
with $`\stackrel{~}{\beta }`$ defined by Eq. (22). Expanding exponents and employing once again (19) leads finally to

$$T_2=\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\frac{\mathrm{\Gamma }[1/(1q)]}{\mathrm{\Gamma }[1/(1q)+nm]}\left[1+(1q)\beta U_q\right]^{\frac{q}{1q}+nm}.$$
(30)
As for $`T_2`$, we recast $`T_1`$ as

$$T_1=\frac{\mathrm{\Gamma }[1/(1q)]}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+it)[1+(1q)\beta U_q]}}{(1+it)^{\frac{1}{1q}}}Tr\left[Z_1(\stackrel{~}{\beta })\widehat{H}\right].$$
(31)
We proceed as follows: i) take advantage of the fact that

$$Tr\left[Z_1(\stackrel{~}{\beta })\widehat{H}\right]=\frac{Z_1(\stackrel{~}{\beta })}{\stackrel{~}{\beta }},$$
(32)
ii) use Eq. (15), and iii) expand again the exponential term. We find

$$T_1=\frac{n\mathrm{\Gamma }[1/(1q)]}{2\pi (1q)\beta }\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^{m+1}}{m!}\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+it)[1+(1q)\beta U_q]}}{(1+it)^{\frac{1}{1q}+n(m+1)}}.$$
(33)
The integral is evaluated as before, which yields

$$T_1=\frac{n}{(1q)\beta }\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^{m+1}}{m!}\frac{\mathrm{\Gamma }[1/(1q)]}{\mathrm{\Gamma }[1/(1q)+n(m+1)+1]}\left[1+(1q)\beta U_q\right]^{\frac{1}{1q}+n(m+1)},$$
(34)
so that the internal energy becomes

$$U_q=\frac{n\xi _n}{(1q)\beta }\frac{\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\frac{1}{\mathrm{\Gamma }[1/(1q)+n(m+1)+1]}\left[1+(1q)\beta U_q\right]^{n(m+1)+1}}{\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\frac{1}{\mathrm{\Gamma }[1/(1q)+nm]}\left[1+(1q)\beta U_q\right]^{nm}},$$
(35)
although compliance with the Tsallis’cut-off condition $`1+(1q)\beta U_q>0`$ is always to be demanded. The above is a non-linear equation to be tackled in numerical fashion. Some pertinent results are displayed in Fig. 1, for different $`q`$-values, they include the standard $`q=1`$ case that yields the Steffan Boltzmann law

$$U_1=\sigma T^4$$
(36)
with $`\sigma =5.67\times 10^8W/m^2K^4`$, the so-called Steffan-Boltzmann constant.

We notice that a Steffan Boltzmann-like law holds, if $`kT`$ is small enough, for all $`q(0<q<1)`$, with $`\sigma =\sigma (q)`$. For high $`T`$-values and $`q1`$ the power law behavior seems to persist, but the corresponding power, let us call it $`a_q`$, is no longer equal to $`4`$. Fig. 2 is a graph of $`a_q`$ versus $`q`$ for

$$U_qT^{a_q}.$$
(37)
Fig. 2 (a) shows $`a_q`$ for small $`kT`$-values. The validity of Stefan Boltzmann’s Law for a wide range of $`q`$-values is easily appreciated. Fig. 2 (b) depicts, instead, $`a_q`$ for large $`kT`$-values. In this limit, important deviations from the $`a_q=4`$ value are detected.

An intermediate temperature range is observed in Fig. 1 in which the power law behavior (of $`U_q`$ vs. $`kT`$) is violated. The corresponding transition is the more abrupt the larger $`|q1|`$. As will be seen in the next section, these violations can be attributed to a Tsallis cut-off in the energy densities.

#### 3 Energy densities

The generalized spectral energy distribution of blackbody radiation $`u_q`$ is given by

$$U_q=_0^{\mathrm{}}𝑑\omega u_q.$$
(38)
In order to obtain $`u_q`$ we work with $`U_q=T_1/T_2`$, with $`T_1`$ and $`T_2`$ given by Eqs. (26) and (27), respectively. The denominator will not suffer any change, but in the numerator we will keep the integral form of $`Z_1`$ as given in Eq. (15). Then, Eq. (32) will lead to

$$Tr\left(e^{\stackrel{~}{\beta }\widehat{H}}\widehat{H}\right)=\mathrm{}A_ne^{\xi _n}\underset{0}{\overset{\mathrm{}}{}}𝑑\omega \omega ^n\frac{e^{\stackrel{~}{\beta }\mathrm{}\omega }}{1e^{\stackrel{~}{\beta }\mathrm{}\omega }},$$
(39)
allowing us to cast $`T_1`$ in the fashion

$$T_1=\frac{\mathrm{\Gamma }[1/(1q)]}{2\pi }\mathrm{}A_n\underset{0}{\overset{\mathrm{}}{}}𝑑\omega \omega ^n\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+it)[1(1q)\beta (\mathrm{}\omega U_q)]}}{(1+it)^{\frac{1}{1q}}}\frac{e^{\xi _n}}{1e^{\stackrel{~}{\beta }\mathrm{}\omega }}.$$
(40)
Recourse to the identity

$$\frac{1}{1e^{\stackrel{~}{\beta }\mathrm{}\omega }}=\underset{s=0}{\overset{\mathrm{}}{}}\left(e^{\stackrel{~}{\beta }\mathrm{}\omega }\right)^s,$$
(41)
yields now

$$T_1=\frac{\mathrm{\Gamma }[1/(1q)]}{2\pi }\mathrm{}A_n\underset{0}{\overset{\mathrm{}}{}}𝑑\omega \omega ^n\underset{s=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑t\frac{e^{(1+it)[1(1q)\beta [\mathrm{}\omega (1+s)U_q]]}}{(1+it)^{\frac{1}{1q}+nm}}.$$
(42)
If we use Eq. (19) the integral above is easily calculated, and we obtain

$$T_1=\mathrm{\Gamma }[1/(1q)]\mathrm{}A_n\underset{0}{\overset{\mathrm{}}{}}𝑑\omega \omega ^n\underset{s=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\frac{\xi _n^m}{m!}\frac{\left[1(1q)\beta [\mathrm{}\omega (1+s)U_q]\right]^{\frac{q}{1q}+nm}}{\mathrm{\Gamma }[1/(1q)+nm]}.$$
(43)
According to (38) the ratio $`\frac{T_1}{T_2}`$ implies that

$$u_q=\mathrm{}A_n\omega ^n\frac{\underset{m=0}{\overset{\mathrm{}}{}}B_mS_m}{\underset{m=0}{\overset{\mathrm{}}{}}B_m},$$
(44)
where

$`B_m`$ $`=`$ $`{\displaystyle \frac{\xi _n^m}{m!}}{\displaystyle \frac{\left[1+(1q)\beta U_q\right]^{nm}}{\mathrm{\Gamma }[1/(1q)+nm]}},`$ (45)
$`S_m`$ $`=`$ $`{\displaystyle \underset{s=1}{\overset{\mathrm{}}{}}}\left[1{\displaystyle \frac{(1q)\beta [\mathrm{}\omega (s+1)U_q]}{1+(1q)\beta U_q}}\right]^{\frac{q}{1q}+nm}.`$ (46)
Figs. 3 and 4 are plots of $`u_q`$ vs. the frequency $`\omega `$. Notice that they quite resemble the semi-empirical Planck law for small $`|q1|`$-values. Fig. 5 clearly shows that the frequency shifts of the maxima of the above referred to curves (located at $`\omega =\omega ^{}`$) comply with the Wien law

$$\omega ^{}T,$$
(47)
for such $`q`$-values. Up to $`q=0.8`$ this nice behaviour persists. As $`q`$ departs from $`1`$ an anomalous behaviour ensues. Local maxima begin to appear in the curves $`u_q(\omega )`$ (Figs. 3 and 4). As $`T`$ grows the global maximum is shifted and the secondary ones become more pronounced, until one of them replaces the former one as the maximum maximorum. An abrupt jump in the curve $`\omega ^{}`$ vs. $`T`$ ensues as a consequence (Fig. 5).

The departures from Planck’s law as $`|q1|`$ grows, could in principle be subject to experimental verification, thus $`q`$ would be fixed.

The above referred to discontinuities in the “Wien-like” behaviour are due to the Tsallis’ cut-off operative in terms of the sums $`T_1`$ and $`T_2`$.

## IV First order corrections

We have seen in the preceding section that (44) yields, for $`q`$-values close to unity, results that quite resemble the ones given by Planck’s law. It is then reasonable to look for a perturbative expansion in $`1q`$.

Let us start with $`S_m`$ (Cf. Eq. (46)) and write it in the form

$$S_m=\underset{s=0}{\overset{\mathrm{}}{}}\mathrm{exp}\left\{\left(\frac{q}{q1}+nm\right)\mathrm{ln}\left[1(1q)\frac{\beta \mathrm{}\omega (1+s)}{1+(1q)\beta U_q}\right]\right\}.$$
(48)
A first order expansion of the logarithm yields

$$S_m\underset{s=0}{\overset{\mathrm{}}{}}\mathrm{exp}\left\{[q+nm(1q)]\frac{\beta \mathrm{}\omega (1+s)}{1+(1q)\beta U_q}\right\},$$
(49)
a power series in $`s`$ of guaranteed convergence for $`1q<1`$. Neglecting the term $`nm(1q)`$ on account of the fact that the concomitant series rapidly converges so that terms for which $`m`$ could compete with $`1/(1q)`$ are indeed negligible, we have

$$S_m\frac{1}{\mathrm{exp}\left[\frac{q\beta \mathrm{}\omega }{1+(1q)\beta U_q}\right]1},$$
(50)
and replacement in Eq. (44) yields

$$u_q\frac{\mathrm{}A_n\omega ^n}{\mathrm{exp}\left[\frac{q\beta \mathrm{}\omega }{1+(1q)\beta U_q}\right]1},$$
(51)
a first order correction for a “generalized Planck law” that yields the classical result for $`q1`$. Eq. (51) provides one then with an approximate energy density.

Using now this equation as a starting point we also get first order corrections in $`1q`$ to the Wien law. We can locate the maximum with respect to $`\omega `$ of Eq. (51) with the auxiliary definition

$$x=\frac{\omega q\beta \mathrm{}}{1(1q)\beta U_q},$$
(52)
and immediately find

$$e^x+\frac{x}{n}=1,$$
(53)
whose solution is a constant $`b`$ for each fixed value of $`n`$. For $`n=3`$ we find, for instance, $`b=2.82`$. The maximum of $`u_q(\omega )`$ for different $`T`$’s is located at distinct frequencies according to

$$\omega _i=\frac{bk}{q\mathrm{}}T_i\frac{1q}{q}\frac{U_q}{\mathrm{}},$$
(54)
a first order correction to Wien’s law.

If we integrate Eq. (51) over frequencies we find $`U_q`$ as a function of $`T`$, and a concomitant first order correction to Planck’s law

$$U_q=\mathrm{}A_n\alpha \left[\frac{1+(1q)\beta U_q}{q\beta \mathrm{}}\right]^{n+1},$$
(55)
that reduces trivially to Planck’s one in the limit $`q1`$. We call $`\alpha `$ the integral

$$\alpha =_0^{\mathrm{}}𝑑x\frac{x^n}{e^x1}=\mathrm{\Gamma }(n+1)\zeta (n+1),$$
(56)
with $`x`$ given by (52).

It may be convenient to point out that the first order relationships here obtained ignore possible cut-off problems. Anyway, in the $`0<q<1`$ range (with $`q1`$) that is relevant to the present considerations, such problems should not arise.

## V Conclusions

We have presented here exact results for the mean energy and the energy density of the blackbody problem in a nonextensive environment. A Stefan Boltzmann-like law  is obtained, for all $`0<q<1`$, although with a $`q`$-dependent SB constant $`\sigma (q)`$, for almost the whole temperature range, although for some temperatures a deviation from the power law behaviour is detected.

The energy density curves correspond to Planck’s law for $`q`$ values close to unity. For other $`q`$’s a gradual departure from the typical Planck curve is observed, which can be attributed to Tsallis’ cut-off effects. Wien’s law is also obtained for $`q`$ values close to unity.

Finally, we have presented first order perturbative results in the limit $`q1`$ for the energy density $`u_q`$ (the parameter is, of course, $`1q`$). These could be of some utility in obtaining first order nonextensive corrections to the Bose-Einstein equation.

###### Acknowledgements.

The financial support of the National Research Council (CONICET) of Argentina is gratefully acknowledged. F. Pennini acknowledges financial support from UNLP, Argentina.


text:
# Charge Conjugation Violation in Supernovae and The Neutron Shortage for R-Process Nucelosynthesis

## Abstract

Core collapse supernovae are dominated by energy transport from neutrinos. Therefore, some supernova properties could depend on symetries and features of the standard model weak interactions. The cross section for neutrino capture is larger than that for antineutrino capture by one term of order the neutrino energy over the nucleon mass. This reduces the ratio of neutrons to protons in the $`\nu `$-driven wind above a protoneutron star by approximately 20 % and may significantly hinder r-process nucleosynthesis.

Core collapse supernovae are perhaps the only present day large systems dominated by the weak interaction. They are so dense that photons and charged particles diffuse very slowly. Therefore energy transport is by neutrinos (and convection).

We beleive it may be useful to try and relate some supernova properties to the symmetries and features of the standard model weak interaction. Parity violation in a strong magnetic field could lead to an asymmetry of the explosion1 . Indeed, supernovae explode with a dipole asymmetry of order one percent in order to produce the very high ‘recoil’ velocities observed for neutron stars2 . However, calculating the expected asymmetry from P violation has proved complicated. Although explicit calculations have yielded somewhat small asymmetries3 ; 4 ; 5  it is still possible that more efficient mechanisms will be found.

In this paper we calculate effects of charge conjugation, C, violation in the Standard Model on the difference between neutrino and antineutrino interactions. In Quantum Electrodynamics C symmetry insures the cross section for $`e^{}p`$ is equal to that for $`e^+p`$ scattering (to lowest order in $`\alpha `$). In contrast, the standard model has large parity, P, and C violations (since the product CP is approxamitely conserved). Therefore the $`\overline{\nu }`$-nucleon cross sections are systematically smaller than $`\nu `$-nucleon cross sections.

However at the low $`\nu `$ energies in supernovae, time reversal symmetry limits the difference between $`\nu `$ and $`\overline{\nu }`$ cross sections. Time reversal can relate $`\nu N`$ elastic scattering and $`\overline{\nu }N`$ where the nucleon scatters from final momentum $`p_f`$ to initial momentum $`p_i`$. If the nucleon does not recoil then the $`\nu `$ and $`\overline{\nu }`$ cross sections are equal. Thus the difference between $`\nu `$ and $`\overline{\nu }`$ cross sections are expected to be of recoil order $`E/M`$ where $`E`$ is the neutrino energy and $`M`$ the nucleon mass. This ratio is relatively small in supernovae. However the coefficient multiplying $`E/M`$ involves the large weak magnetic moment of the nucleon (see below).

The standard model has larger $`\nu `$ cross sections than those for $`\overline{\nu }`$. For neutral currents, this leads to a longer mean free path for $`\overline{\nu }_x`$ compared to $`\nu _x`$ (with x=$`\mu `$ or $`\tau `$). Thus even though $`\nu _x`$ and $`\overline{\nu }_x`$ are produced in pairs, the antineutrinos escape faster leaving the star neutrino rich. The muon and tau number for the protoneutron star in a supernova could be of order $`10^{54}`$6 . Supernovae may be the only known systems with large $`\mu `$ and or $`\tau `$ number. For charged currents, the interaction difference can change the equilibrium ratio of neutrons to protons and may have important implications for nucleosynthesis. We discuss this below. To our knowledge, all previous work on nucleosynthesis in supernovae assumed equal $`\nu `$ and $`\overline{\nu }`$ interactions (aside from the n-p mass difference).

The neutrino driven wind outside of a protoneutron star is an attractive site for r-process nucleosynthesis7 . Here nuclei rapidly capture neutrons from a low density medium to produce heavy elements8 . This requires, as a bare minimum, that the initial material have more neutrons than protons. The ratio of neutrons to protons n/p in the wind depends on the rates for the two reactions:

$$\nu _e+np+e^{},$$
$`(1a)`$
$$\overline{\nu }_e+pn+e^+.$$
$`(1b)`$
The standard model cross sections for Eqs. (1a,1b) to order $`E/M`$ are,

$$\sigma =\frac{G^2\mathrm{cos}^2\theta _c}{\pi }(1+3g_a^2)E_e^2[1\gamma \frac{E}{M}\pm \delta \frac{E}{M}],$$
$`(2)`$
with $`G`$ the Fermi constant (and $`\theta _c`$ the Cabbibo angle), $`E_e=E\pm \mathrm{\Delta }`$ the energy of the charged lepton and $`\mathrm{\Delta }=1.293`$ MeV is the neutron-proton mass difference. The plus sign is for Eq. (1a) and the minus sign for Eq. (1b). We use $`g_a1.26`$. Equation (2) neglects small corrections involving the electron mass and coulomb effects while the finite nucleon size only enters at order $`(E/M)^2`$.

We refer to the $`\gamma `$ term as a recoil correction. It is the same for $`\nu `$ and $`\overline{\nu }`$, $`\gamma =(2+10g_a^2)/(1+3g_a^2)3.10`$. Finally, the $`\delta `$ term, $`\delta =4g_a(1+2F_2)/(1+3g_a^2)4.12`$, involves the interference of vector (1+2$`F_2`$) and axial ($`g_a`$) currents. This violates P, which by CP invariance also violates C. This increases the $`\nu `$ and decreases the $`\overline{\nu }`$ cross section, Note, $`F_2`$ is the isovector anomalous moment of the nucleon. (This is the weak magnetism contribution.)

The equilibrium electron fraction per baryon $`Y_e`$ (which is equal to the proton fraction assuming charge neutrality) is simply related to the rate $`\overline{\lambda }`$ for Eq. (1b) divided by the rate $`\lambda `$ for Eq. (1a).

$$Y_e=(1+\frac{\overline{\lambda }}{\lambda })^1$$
$`(3)`$
The ratio neutrons to protons is, $`\frac{n}{p}=\frac{1}{Y_e}1`$.

In ref. prl  we calculate the reaction rates by averaging Eq. (2) over neutrino spectra to get,

$$Y_e=\left(1+\frac{L_{\overline{\nu }_e}\overline{ϵ}}{L_{\nu _e}ϵ}QC\right)^1.$$
$`(4)`$
Here $`ϵ`$ ($`\overline{ϵ}`$) is the $`\nu _e`$ ($`\overline{\nu }_e`$) mean energy, $`L_{\nu _e}`$ ($`L_{\overline{\nu }_e}`$) is the $`\nu _e`$ ($`\overline{\nu }_e`$) luminosity, $`Q`$ is the correction from the reaction Q value,

$$Q=\frac{12\frac{\mathrm{\Delta }}{\overline{ϵ}}+a_0\frac{\mathrm{\Delta }^2}{\overline{ϵ}^2}}{1+2\frac{\mathrm{\Delta }}{ϵ}+a_0\frac{\mathrm{\Delta }^2}{ϵ^2}},$$
$`(5)`$
and with C violating one has a factor $`C`$,

$$C=\frac{1(\delta +\gamma )a_2\frac{\overline{ϵ}}{M}}{1+(\delta \gamma )a_2\frac{ϵ}{M}}.$$
$`(6)`$
Simply evaluating Eq. (6) for typical parameters yields $`C0.8`$. Thus, the difference between $`\nu `$ and $`\overline{\nu }`$ interactions reduces the equilibrium n/p ratio by approximately 20 %. This is an important result and will be discussed below.

Evaluating Eq. (4) for the neutrino fluxes of a Supernova simulation by Wilson 11  shows that with the $`C`$ term the neutrino driven wind starts out proton rich and ends up with about equal numbers of neutrons and protons. When C violation is included the wind is never significantlt neutron rich. It is very unlikely that successful r-process nucleosynthesis can take place in the wind of this or similar models.

With the approximately 20 % reduction in n/p from the difference between $`\nu `$ and $`\overline{\nu }`$ interactions, there appears to be very serious problems with r-process nucleosynthesis in the wind of present supernova models. In addition to the initial lack of neutrons, one has to overcome the effects of neutrino interactions during the assembly of $`\alpha `$ particles and during the r-process itself14 . These further limit the available neutrons per seed nucleus. Thus, it is unlikely that present wind models will produce a successful r-process. Of course, the wind in supernovae may not be the r-process site, although this may be unappealing (see for example8 ; 15 ). If the wind is not the site, one must look for alternative environments.

However, the effects of neutrino interactions may be very general. The only requirement is that energy transport from neutrinos plays some role in helping material out of a deep gravitational well. Given this, it is quite likely that the n/p ratio will be determined by the relative rates of Eqs. (1a,1b). Therefore differences in $`\nu `$ and $`\overline{\nu }`$ interactions may be important for just about any nucleosynthesis site that involves neutrinos.

If the $`\nu `$-driven wind is the r-process site, it is very likely, present models of the neutrino radiation in supernovae are incomplete. The high values of $`Y_e`$ make it almost impossible to have a successful r-process by only changing matter properties, such as the entropy. The neutrino fluxes will (almost assuredly) need to be changed.

Changes in the astrophysics used in the simulations or new neutrino physics such as neutrino oscillations17  could change $`\overline{ϵ}`$, $`ϵ`$ and or the luminosities and lead to a more neutron rich wind. The oscillations of more energetic $`\overline{\nu }_x`$ with $`\overline{\nu }_e`$ could increase $`\overline{ϵ}`$. However, we have some information on the $`\overline{\nu }_e`$ spectrum from SN1987a18 . Thus one can not increase $`\overline{ϵ}`$ without limit. Indeed if anything, the Kamiokande data suggest a lower $`\overline{ϵ}`$. Any model which tries to solve r-process nucleosynthesis problems by increasing $`\overline{ϵ}`$ should first check consistency with SN1987a observations19 . Alternative modifications could include oscillations of $`\nu _e`$ to a sterile neutrino or a lowering of $`ϵ`$. (However, we know of no model which lowers $`ϵ`$.) Whatever the modification of the neutrino fluxes, one will still need to include the differences between $`\nu `$ and $`\overline{\nu }`$ interactions in order to accurately calculate n/p.

In conclusion, supernovae are one of the few large systems dominated by energy transport from weakly interacting neutrinos. Therefore, some supernova properties may depend on symmetries and features of the standard model weak interactions. The cross secton for neutrino capture is larger than that for antineutrino capture by a term of order the neutrino energy over the nucleon mass. This difference between neutrino and antineutrino interactions reduces the ratio of neutrons to protons in the $`\nu `$-driven wind above a protoneutron star by approximately 20 % and may significantly hinder r-process nucleosynthesis.

This work was supported in part by DOE grant: DE-FG02-87ER40365.