format: text

text:
# Search for cosmic ray sources using muons detected by the MACRO experiment

## 1 Introduction

The interest in the search for cosmic ray point sources identified from a measurable flux of underground muons has historical motivations mainly because of the “CygX-3 saga”. Cygnus X-3 (CygX-3) is a galactic binary system well studied in all types of electromagnetic radiation, most notably in the X-rays. At $`\gamma `$-ray energies, $`E10^{15}`$ eV, the Kiel Extensive Air Shower array (EAS) initially reported an excess of events from CygX-3 correlated with its 4.8 hour binary period. This observation appeared to confirm previous results at lower energies . Subsequently, several groups operating EAS array experiments and atmospheric Čerenkov telescopes confirmed the signals with different statistical significance and at different energy thresholds . These results prompted a new generation of more sensitive experiments using new techniques. The new experiments reported clear evidence of TeV $`\gamma `$-ray emissions from many galactic ($`i.e.`$, Crab and SN1006 ) and the extragalactic AGN-like ($`i.e.`$, Mrk421 and Mrk501 ) sources. Finally, in 1990-1991 the new generation detectors CASA-MIA and CYGNUS put stringent upper limits to the $`\gamma `$-ray signal from CygX-3, which excluded early observations.

Soon after the initial EAS detection from CygX-3, two underground experiments (Soudan 1  and NUSEX ) reported excesses in TeV muons pointing back to the CygX-3 direction when correlating event arrival time with the known orbital period of the source. These detections suggested new physics beyond the Standard Model since muons are the products of the decay of charged pions and kaons, which are produced by primary Cosmic Ray (CR) interactions with atmospheric nuclei. Charged CR nuclei cannot propagate directly from CygX-3 across 10 kpc to the Earth and there is only a very low probability of UHE $`\gamma `$-rays producing TeV muons underground. Investigations by other underground detectors , however failed to detect any significant muon excess from that source or any other.

The MACRO experiment, which ended its operations life in December 2000, was a large area underground detector able to reconstruct muon arrival directions to very high accuracy . The apparatus started operations in 1989 and was completed in April 1994. A preliminary investigation using a limited data sample of $`1.8\times 10^6`$ muons collected by its first and second supermodules was published in 1993 . Since then, MACRO has increased the muon sample by a factor of 30, its pointing capability has been accurately determined and its direction reconstruction capability carefully studied.

In this paper we present the final results on the search for a muon excess with respect to the evaluated background by surveying the sky in declination from -15 to 90 and from several candidate UHE gamma ray emitters. We also present a more sensitive search for CygX-3 and HerX-1 using their known periodicities. Finally we searched for bursting behavior of Mrk421 and Mrk501.

## 2 The MACRO detector

MACRO was a multipurpose underground detector located in the Gran Sasso underground Laboratory (LNGS-INFN). It was designed to search for rare events in the cosmic radiation and its sensitivity was optimized to detect supermassive magnetic monopoles. The detector  had a modular structure and dimensions of $`76.5\times 12\times 9.3`$ m3 with a total acceptance for an isotropic flux of about 104 m2sr. The rock overburden had a minimal depth of 3150 m.w.e and an average of 3700 m.w.e.

The minimal energy of surface muons that trigger the MACRO apparatus was about 1.3 TeV. The rock coverage was very irregular and the actual slant depth under which a source is observable was taken into account when evaluating the muon flux at the surface. The full detector acceptance for downgoing muons with zenith angles $`72^{}`$ was about 3100 m2sr.

The detector has worked with different configurations starting with 1 supermodule at the beginning and finally with the full configuration consisting of 6 supermodules and the upper half called the “Attico”. The active detection elements are planes of streamer tubes for tracking and liquid scintillation counters for fast timing. The lower half of the detector is filled with streamer tube planes alternating with trays of crushed rock. The upper part is hollow and contains the electronics racks and work areas. There are 10 horizontal streamer tube planes in the bottom half of the detector, and 4 planes on the top, all with wires and 27 strips, providing stereo readout of the detector hits. Six vertical planes of streamer tubes cover each side of the detector.

## 3 Data selection

The data used in this analysis were collected in the period May 1989- December 2000. The total number of recorded muons is approximately 60 million.

In order to optimize the quality of the tracking reconstruction and to eliminate periods when the detector was malfunctioning, we applied several cuts on a run-by-run basis. These run cuts include a check on the efficiency of the streamer tubes system whose data were used for the track reconstruction, and a check on the event rate. The streamer tube efficiency was obtained using the sub-sample of tracks crossing all the 10 lower streamer tube planes; runs having an average efficiency smaller than 90$`\%`$ for the wire view and 85$`\%`$ for the strip view were discarded. As the detector configuration changed during the data taking, the average value of the counting rate was computed for each run. Those with abnormal rates, $`i.e.`$, runs having muon counting rate that deviated by more than $`\pm 3\sigma `$ from the average, were cut.

In addition to the described run cuts, we applied event cuts by excluding:

* events having zenith angles larger than 72; the cut is due to the large uncertainties in the rock depth crossed by these events;
* events whose reconstruction in one of the two streamer tube stereo views is missing;
* muons in a bundle with multiplicity larger than 2; this cut is necessary because the reconstructed direction in space for high multiplicity events is not completely reliable;
* events that had no arrival time given by the Universal Time Clock (UTC);
* events that do not cross at least 3 planes of the lower part of the apparatus; this cut excludes low energy muons  coming from possible secondary interactions and it improves the angular resolution. The minimal muon energy to satisfy this cut is 1 GeV.

These cuts reduce the sample to 49.9 million well-reconstructed single and double muon events during 74,073 hours of livetime. Table 1 shows, for different data taking periods, the detector configuration, the number of events which survived the cuts and the effective livetime.

## 4 Background estimation

Several methods were suggested to evaluate the background when looking for an excess of counts from a fixed direction in the sky . Many experiments simply consider, as an estimation of the background, the average value of the counts of the surrounding sky bins with respect to a selected one . An alternative way is to average the counts of all the bins at the same declination $`\delta `$, except the chosen one in a particular value of the right ascension (RA). The assumption that cosmic rays arrive uniformly from any direction in the sky is implicit in this method. For an underground detector both methods require that the unequal distribution of the overburden and the dead-time of the apparatus be taken carefully into account .

In our analysis, we adopted a different approach. Assuming that the arrival directions and the arrival time of underground muons are uncorrelated variables, we constructed sets of simulated events by randomly coupling the times and the directions of each event in local coordinates. A mandatory requirement for this method is a good accuracy in the measurement of the arrival direction and of the recorded time. For this reason, we excluded from the analysis all the runs and the events with errors in the readout time, as described in the previous section. We determined that the optimum number of background events to be associated with each real event, minimizing the computer processing time, is 25. A total number of $`1.3\times 10^8`$ simulated muons were generated as background. Figure 1 presents the distribution of observed events $`vs.`$ the right ascension RA (a) and the sine of the declination sin$`\delta `$ (b) together with the simulated background. The fluctuations in the simulated distributions are small compared with those in the data. The unevenness in RA is due to the effect of the dead times. The sin$`\delta `$ distribution reflects the shape of the mountain and also the exposure of MACRO (see ).

## 5 Best sky bin definition

It is usual for astronomical telescopes to define a point spread function (PSF) that in most cases is a simple bi-dimensional normal distribution. The PSF can be used to define the optimal source bin, $`i.e.`$, the half-angle $`\theta `$ of a cone centered on a source giving the maximum signal over background ratio $`\frac{S}{\sqrt{B}}`$. For a normal PSF distribution with variance $`\sigma `$ it is shown in  that the best value of the half-angle of a cone is $`\theta =1.58\sigma `$. This cone contains, on average, 72$`\%`$ of the total number of events.

The MACRO PSF was defined using the double muon sample. A double muon event (as in general multiple muon events) is produced by an interaction of a primary nucleus at the top of the atmosphere; muons are expected to arrive at the earth surface practically parallel. When reaching an underground detector, the reconstructed spatial directions of these two muons differ as a consequence of the independent scattering of each muon. We can therefore estimate the MACRO PSF by using as a variable the distribution of the angular separation of the reconstructed muon directions $`\theta `$ for double muon events. This value must be divided by a $`\sqrt{2}`$ factor to account for the two independent scatterings .

We found  that 50$`\%`$ of double muons events are contained in a cone of 0.5 half-angle, and that 72$`\%`$ of events are contained in a cone of 1.05 half-angle. Initially, this result suggests that the best value for the MACRO search bin (assuming that the $`\theta `$ is normally distributed) is an half-angle cone of 1.05/$`\sqrt{2}`$. To produce the bidimensional ($`RA,\delta `$) distribution representing our PSF (shown in Fig. 2) we used the differences between the reconstructed muon directions in RA and declination $`\delta `$, divided by $`\sqrt{2}`$. The main contribution to the differences between the reconstructed directions is the multiple Coulomb scattering (MCS) of muons in the rock. The PSF shown in Fig. 2 has a non-gaussian shape, with more events for small displacements and longer tails. The latter were produced by stochastic interaction process in the muon energy loss.

To obtain the value of the best search cone angle, we performed a dedicated Monte Carlo. A source was simulated in a particular ($`RA,\delta `$) cell, and the events were extracted according to the PSF of Fig. 2. The background events were generated according to a flat density distribution. We calculated the ratio $`\frac{S}{\sqrt{B}}`$ as a function of the search cone angle, as shown in Fig. 3. The maximum of the $`\frac{S}{\sqrt{B}}`$ ratio is reached for a search cone angle of $``$ 0.4-0.5, which thus will represent the best choice value for MACRO.

As an independent check of the validity of the simulation, we studied the signal due to the Moon shadowing effect (in this case, a lack of events) measured by our detector . The cumulative number of missing events obscured by the Moon’s disk, as function of the angular distance from the Moon center, shows maximum evidence for a deficit using a cone of 0.45 half-angle.

From all these indications the best value for the search cone angle was chosen to 0.5. However, in the analyses with very low bin contents such as those used in the search for flaring activities, we use an enlarged search window due to poor statistics.

## 6 Search for DC sources

In the search for a steady excess of muons from any direction of the sky (DC point-like sources search) we performed an all-sky survey without a priori assumptions. We divided the sky into 37176 bins of equal solid angle ($`\mathrm{\Delta }`$$`\mathrm{\Omega }`$= 2.3$`\times `$ $`10^4`$ sr; $`\mathrm{\Delta }`$RA = $`1^{}`$, $`\mathrm{\Delta }`$sin$`\delta `$ = 0.013). These bins have the same $`\mathrm{\Delta }`$$`\mathrm{\Omega }`$ as a narrow cone of half-angle $`0.5^{}`$. We examined the sky bins looking for significant deviations from the simulated background. The deviation was defined as $`\sigma =\frac{(n_{obs}n_{exp})}{\sqrt{n_{exp}}}`$, where $`n_{obs}`$ is the number of events observed in each bin and $`n_{exp}`$ is the expected number of background events in that bin.

We used three different sky grids, each displaced by half bin width in RA, in sin($`\delta `$) and in both coordinates, to take into account the possibility that a source can be located at the edge of one of the bins. Fig. 4 shows the distributions of the deviations for the four sky-maps with the best-fit Gaussian function superposed. The positive deviations in the first map are reported in the bidimensional map in RA and sin$`\delta `$ of Fig. 5. No particular pattern or clustering of positive deviations is observed. The line indicating the galactic plane is also shown. We conclude that there is no evidence for a steady source emitting muons in our data.

We calculated the 95$`\%`$ C.L muon flux for all sky bins using the formula :

$$\mathrm{J}_\mu ^{stdy}(95\%)\frac{n_\mu (95\%)}{KA_{eff}T_{exp}}\mathrm{cm}^2\mathrm{s}^1$$
(1)
where:

1. $`n_\mu `$(95%) is the upper limit for the number of muons in the bin at 95%C.L., where $`n_{obs}`$ and $`n_{exp}`$ are respectively the number of observed and expected events in that bin. $`n_\mu `$(95%) was calculated, according to , as the value for which:

$$\frac{2}{\sqrt{\pi }}_{n_\mu (95\%)}^{\mathrm{}}\frac{e^{(n_\mu \overline{n_\mu })^2/2\sigma ^2}}{\sqrt{2}\sigma }\mathrm{d}n_\mu =0.05$$
(2)
with $`\overline{n_\mu }`$=$`n_{obs}n_{exp}`$ and $`\sigma ^2`$=$`n_{exp}`$.
2. $`A_{eff}(i)`$ is the average effective detector area for every bin. It was computed by averaging the projected area seen by each muon and taking into account the geometrical and the tracking reconstruction efficiencies:

$$\mathrm{A}_{eff}(i)=\frac{1}{n_{obs}(i)}\underset{j=1}{\overset{n_{obs}(i)}{}}A(RA_i,\delta _i)$$
(3)
where $`RA_i`$ and $`\delta _i`$ are the muon arrival right ascension and declination.
3. $`T_{exp}`$ is the exposure time computed as the time in which each bin is visible. A bin is not visible in our apparatus, if its zenith angle is larger than 72.
4. K is a correction factor which takes into account the fraction of muons within the bin dimensions for an hypothetical source placed at the bin center. We calculated K=0.5 for a search cone of 0.5 half-angle and K=0.72 for 1.0 half-angle.

The distribution of the 95% C.L. upper limit $`J_\mu ^{stdy}`$(95%) for all the 37176 bins in the sky region accessible for MACRO has an average value of $`2.3\times 10^{13}\mathrm{cm}^2\mathrm{s}^1`$. For almost all bins, the upper limit ranges between $`1\times 10^{13}J_\mu ^{stdy}(95\%)4\times 10^{13}\mathrm{cm}^2\mathrm{s}^1`$.

Finally we also investigated some selected interesting point sources (in Table 2) identified by surface telescopes and EAS arrays . None of the selected sources in the list exhibits a significant deviation from the background.

## 7 Searches for modulated signals from CygX-3 and HerX-1

Periodic UHE gamma ray sources are attractive from the observational point of view. When it is known that the emission from the source is modulated with a certain periodicity, the signal-to-noise ratio is improved by $`\sqrt{N}`$  by folding the event arrival time modulo the source period into N bins. This folding procedure was employed to analyze the data from the direction of CygX-3 and HerX-1 .

We used the quadratic ephemeris reported in Table 3 and obtained from the fit to the observed X-ray light curve , . The phase diagrams for events coming from a cone of 1, centered on the position of the two sources, are shown in Fig. 6. The expected number of background muons in each phase bin is also shown with the dashed lines.

The largest positive deviation (1.8 $`\sigma `$) in the CygX-3 cycle is in the phase bin 0.1-0.2. The largest positive deviation for HerX-1 (1.7 $`\sigma `$) is in the phase bin 0.7-0.8. Using the values of the largest fluctuation we computed the 95$`\%`$ C.L. upper limits to the modulated muon flux using equation 1, and the results are reported in Table 3. In Fig. 7 our computed upper limit for the modulated emission from CygX-3 is compared with those of previous underground experiments. Our level is the lowest reached by such detectors. The problem of understanding past positive observations still remain unsolved.

## 8 Search for flaring activity from Mkn421 and Mkn501

During the 1995-1999, several Čerenkov telescopes reported observations of flaring activities up to 10-20 TeV from the two celestial objects Mkn421 and Mkn501. These observations prompted us to search for bursting muon signals in a 1 half-cone, around the position of these sources of U.H.E. photons. We used two different methods.

In the first method we studied, similarly to , the accumulation rate of events coming from each of the two sources, adding day-by-day the differences between the measured number of events and the calculated background. Figures 8 and 9 show the cumulative excesses as a function of date (Julian days) since MACRO starting data taking. The cumulative excess presents fluctuations, but never becomes significant.

In the second method we assumed (as in ) that the background has a Poissonian distribution. If $`n_{exp}`$ is the expected daily background from the direction of a source, then the probability to observe a random fluctuation of the background large as the observed $`n_{obs}`$ events in a day, is given by:

$$P=1\underset{n=0}{\overset{n_{obs}1}{}}\frac{\alpha ^n}{(1+\alpha )^{n+n_{exp}+1}}\frac{(n+n_{exp})!}{n!\times n_{exp}!}$$
(4)
where $`\alpha `$ is the ratio of the ON-source time to the OFF-source time. Because we extract 25 simulated background events for each real one, $`\alpha `$ is 0.04. For a Poissonian background, the cumulative frequency distribution of P is expected to be a power law with index -1, in logarithmic scale.

Figures 10 and 11 show the value of -$`\mathrm{log}P`$ evaluated day-by-day for Mkn421 and Mkn501, respectively. In Fig. 10 the date of the largest fluctuations for Mkn421 with respect to the background are also indicated.

To verify the probability of such a positive fluctuation from Mkn421, we computed the quantity -$`\mathrm{log}P`$ for a set of selected bins having an exposure similar to the two Mkn objects, each monitored for about 3600 days. Figure 12 shows the log-log plot of the cumulative frequency distribution for all these sky bins. Since this cumulative distribution has a slope close to -1, as expected in the case of no source detected, we can use the Poissonian statistics to compute the expectation to see large fluctuations.

We set as an attention level a probability of $`10^3`$ for a fluctuation. We found 4 days with probability larger than $`10^3`$ for Mkn421, and none for Mkn501. We observed the Mkn421 source for 3600 days: assuming as attention level for any source the probability value of $`10^3`$, the expected number of random fluctuations with probability lower than $`10^3`$ is $`3600(1\times 10^3)=3.6`$. Therefore the Poissonian probability to observe 4 random fluctuation in the same period, with an average of 3.6, is about 20$`\%`$. This probability value therefore excluded a positive observation of a burst from this source.

## 9 Conclusions

Since February 1989, the MACRO detector collected 49.9 million well-constructed muons. Using this sample, we searched for muon excesses above background from all visible sky directions and from known astrophysical sources. No significant excesses were found from the all-sky survey. We computed the 95% confidence level upper limit $`J_\mu ^{stdy}(95\%)`$ for a steady muon flux for all the 37176 sky cells; the average value of $`J_\mu ^{stdy}(95\%)`$ is equal to 2.3$`\times `$$`10^{13}`$ $`\mathrm{cm}^2\mathrm{s}^1`$. We analyzed the muons coming from the direction of Cyg X-3 and Her X-1, searching for a modulated emission, with a negative result. The search for a steady or episodic emission coming from Mkn421 and Mkn501 was made with two different methods. We found no muon excess above the estimated background.

## 10 Acknowledgments

We gratefully acknowledge the support of the director and of the staff of the Laboratori Nazionali del Gran Sasso and the invaluable assistance of the technical staff of the Institutions participating in the experiment. We thank the Istituto Nazionale di Fisica Nucleare (INFN), the U.S. Department of Energy and the U.S. National Science Foundation for their generous support of the MACRO experiment. We thank INFN, ICTP (Trieste), WorldLab and NATO for providing fellowships and grants (FAI) for non Italian citizens.

This work is dedicated to the memory of Lynn Miller.


text:
# Comment on “Entropy and Wigner function”

## Abstract

In a recent paper \[Phys. Rev. E 62, 4556 (2000)\] Manfredi and Feix proposed an alternative definition of quantum entropy based on Wigner phase-space distribution functions and discussed its properties. They proposed also some simple rules to construct positive-definite Wigner functions which are crucial in their discussion. We show, however, that the main results obtained and their discussion are incorrect.

Entropy is without any doubt one of the most important physical concepts, not only in thermodynamics and statistical mechanics, but also in various areas of quantum theory. It remains to be also one of the most mysterious quantities, especially outside the traditional phenomenological thermodynamics.

In a recent paper PRE\_62\_4665  Manfredi and Feix proposed an alternative definition of quantum entropy with quadratic dependence on the Wigner distribution function $`W(x,p)`$:

$$S_2=1(2\pi \mathrm{})^D_\mathrm{\Gamma }𝑑x𝑑p(W(x,p))^2$$
(1)
The main advantage of this expression is that it is an exact Weyl transform of the linear entropy:

$$S_2=1\mathrm{T}r(\widehat{\rho }^2)$$
(2)
and therefore all properties of the linear entropy could be easily extended to Eq.(1).

For a quantum state per se, i.e., without any measurements being involved, represented through a density operator $`\widehat{\rho }`$, the (dimensionless) quantum entropy may be defined after von Neumann GN\_1927\_279  as:

$$S=\mathrm{T}r(\widehat{\rho }\mathrm{ln}\widehat{\rho })$$
(3)
The presence of the $`\mathrm{ln}\widehat{\rho }`$ term precludes an easy translation of Eq.(3) to, e.g., the Weyl-Wigner-Moyal phase-space picture of quantum mechanics.

The classical limit of the von Neumann entropy is given by the Wehrl entropy RepMathPhys\_16\_353 ; JMP\_25\_1507 :

$$S_{cl}=_\mathrm{\Gamma }dpdq\rho _H(p,q;\kappa ))\mathrm{ln}\rho _H(p,q;\kappa )$$
(4)
where $`\rho _H`$ denotes the respective Husimi phase-space distribution function parametrized with some arbitrary positive constant $`\kappa `$. The Wehrl entropy is possibly the closest quantity to the classical Gibbs entropy (Eq.(4) in PRE\_62\_4665 ), but they are obviously very different from each other. Therefore, the discussion of the classical limit of the von Neumann entropy presented in Sec. I of PRE\_62\_4665 , although appealing, must be considered as incorrect.

The linear entropy Eq.(2), as well as the term $`\mathrm{T}r(\widehat{\rho }^2)`$ alone, is obviously a measure of the purity of the involved quantum state, therefore the opinion (PRE\_62\_4665 , p. 4666) that $`W^2𝑑x𝑑p`$ “has no obvious physical meaning” seems to be surprising.

Discussing the properties of the quantum entropy, the authors have found that the Wigner distribution function which maximizes the entropy is defined as follows (PRE\_62\_4665 , p. 4668): “ $`W=\mathrm{c}onst=\mathrm{\Omega }^1`$ within a phase space region of volume (area) equal to $`\mathrm{\Omega }`$, and $`W=0`$ elsewhere”, i.e., $`W`$ is simply a characteristic function of this phase space area normalized to unity. It could be shown that such a phase-space function, and in general any compactly supported phase-space function cannot be regarded as a valid Wigner distribution function (cf. JJW\_PLA133\_459 ; DL\_PhysicaA182\_643 ; JFAA4\_723  and references therein).

In Sec. IV (PRE\_62\_4665 , p. 4669) the authors discuss extensively the Gaussian smoothing of Wigner distribution functions. A claim is made that the result of Gaussian smoothing of any Wigner distribution function “… is still an admissible Wigner function”. The admissibility is then defined in exactly the same way as the necessary and sufficient conditions for a phase-space function to be a Wigner distribution function (Ref. 4 in PRE\_62\_4665 , see also PRA\_34\_1 ). In other words, any admissible phase-space function have to be a genuine Wigner distribution function and vice versa.

It is well known (cf., e.g., PRep\_259\_147 ) that the Gaussian smoothed Wigner phase-space distribution functions are Husimi functions, which are associated with the (generalized) antinormal operator ordering. Also, Husimi functions does not give proper marginal densities in the general case. In contrast, the Weyl-Wigner representation employs a different, symmetrized operator ordering, and the Wigner distribution functions give always proper marginal densities both in position and momentum spaces. Therefore, the claim that Gaussian smoothing of any Wigner distribution function gives always a non-negative Wigner distribution function cannot be true as inconsistent with well known results, and the offered proof must be flawed.

Indeed, the key point in the authors’ reasoning is the substitution of $`W(x,p)`$, corresponding to $`\psi (x)`$, by $`W_1(x,p)=W(x,p)`$ (corresponding to $`\psi (x)`$) made in Eq.(33) of PRE\_62\_4665 . This substitution is not legitimate in the general case, because only $`W(x,p)`$ obeying $`W(x,p)=W(x,p)`$ could be allowed here.

To demonstrate that explicitly, let us repeat here orderly the steps performed in PRE\_62\_4665 , starting from the l.h.s of Eq. (32) with $`\overline{W}(x,p)`$ taken from Eq. (30):

$`{\displaystyle 𝑑x𝑑p\overline{W}(x,p)F(x,p)}`$ (5)
$`={\displaystyle 𝑑x𝑑pF(x,p)𝑑x^{}𝑑p^{}W(x^{},p^{})K(xx^{},pp^{})}`$
$`={\displaystyle 𝑑x𝑑pF(x,p)𝑑x^{\prime \prime }𝑑p^{\prime \prime }W(xx^{\prime \prime },pp^{\prime \prime })K(x^{\prime \prime },p^{\prime \prime })}`$
$`={\displaystyle 𝑑x^{\prime \prime }𝑑p^{\prime \prime }K(x^{\prime \prime },p^{\prime \prime })𝑑x𝑑pW(xx^{\prime \prime },pp^{\prime \prime })F(x,p)}`$
The last expression Eq. (5) obtained above is equivalent to

$$𝑑x^{\prime \prime }𝑑p^{\prime \prime }K(x^{\prime \prime },p^{\prime \prime })𝑑x𝑑pW(x^{\prime \prime }x,p^{\prime \prime }p)F(x,p)$$
(6)
containing the desired convolution $`WF(x^{\prime \prime },p^{\prime \prime })`$, if and only if:

$$W(x,p)W(x,p)$$
(7)
The suggested further (PRE\_62\_4665 , p. 4670) way: “to construct a phase-space distribution function which is both positive and admissible \[…\] by smoothing an arbitrary (but positive) function of phase space variables, again with a Gaussian kernel. …” is also incorrect. This should be not a surprise even at a first glance, because an arbitrary, but non-negative phase-space function will be almost always not a sufficiently regular one. A simple counterexample may be furnished here by the function:

$$f_a(x,p)=\mathrm{exp}\left\{a(x^2+p^2)\right\},a>a_0>0$$
(8)
where the parameter $`a`$ is sufficiently large (the threshold value $`a_0`$ depends on the smoothing Gaussian). This everywhere positive phase-space function cannot be smoothed with a Gaussian, because the respective integrals are divergent.

Finally, the procedure of taking the classical limit of Eq.(1) by “keeping the Wigner function fixed, and letting Planck’s constant go to zero”, discussed in Sec. VI (PRE\_62\_4665 , p. 4673) cannot be so simple, due to the explicit dependence of the Wigner distribution function on the Planck’s constant. An excellent and detailed discussion of these and related problems may be found in PTRSLA\_287\_237 . However, for pure states we have:

$$_\mathrm{\Gamma }𝑑x𝑑pW(x,p)^2=(2\pi \mathrm{})^D$$
(9)
which in combination with Eq.(1) gives always $`S_2=0`$ rather than $`S_2=1`$ (indicated as the “unpleasant property” of $`S_2`$ in PRE\_62\_4665 ).

###### Acknowledgements.

It is a pleasure to thank Professor Jens Peder Dahl for his extremely warm hospitality. The support received from the Danish National Research Council is also acknowledged.
Addendum

This Comment has been rejected by the PRE Editorial Board as not “\[…\] helpful in straightening the matter \[…\]” (quoted from the Editorial Board Member report).

Whatever it is supposed to mean, let me only stress here once again the most important issues raised in my Comment:

* Gaussian-smoothed Wigner distribution functions are in general not Wigner distribution functions themselves.
* Arbitrary positive functions of phase-space variables cannot be in general converted to Wigner distribution functions via Gaussian smoothing.

These issues are crucial for the criticized paper, where quantum entropy is defined as a state functional in the Weyl-Wigner-Moyal picture (cf. Eq.(8) in PRE\_62\_4665 ). The usage of Gaussian-smoothed Wigner distribution functions instead of ordinary Wigner distribution functions means here a representation switch to the Husimi picture, with all consequences for the derived results. For example, the inequality $`S_2[\overline{W}]S_2[W]`$ (cf. Eq.(44) in PRE\_62\_4665 ) does not indicate that smoothing generally increases the $`S_2`$ entropy, because $`S_2[\overline{W}]`$ and $`S_2[W]`$ are well-defined in distinct phase-space representations.


text:
# Light bullets and optical collapse in vacuum

## 1 Introduction

According to quantum electrodynamics (QED), the non-classical phenomenon of photon–photon scattering can take place due to the exchange of virtual electron–positron pairs. This is a second order effect (in terms of the fine structure constant $`\alpha e^2/4\pi \epsilon _0\mathrm{}c1/137`$), which in standard notation can be formulated in terms of the Euler–Heisenberg Lagrangian density

$$_{\mathrm{EH}}=\epsilon _0+\xi (4^2+7𝒢^2),$$
(1)
where $`\xi 2\alpha ^2\epsilon _0^2\mathrm{}^3/45m_e^4c^5`$, $`\frac{1}{2}(E^2c^2B^2)`$, $`𝒢c𝐄𝐁`$, $`e`$ is the electron charge, $`c`$ the velocity of light, $`2\pi \mathrm{}`$ the Planck constant and $`m_e`$ the electron mass. The latter terms in expression (1) represent the effects of vacuum polarization and magnetization. We note that $`=𝒢=0`$ in the limit of parallel propagating waves. It is therefore necessary to use other waves in order to obtain an effect from the QED corrections. Furthermore, it has recently been shown that QED vacuum nonlinearities can lead to self-focusing of laser beams , and that such effects might be experimentally viable within a decade. In Ref.  it was found that guiding structures can be useful for the purpose of detecting QED effects, and a suggestion was presented for an experimental setup that gives a measurable signal.

In the present work we thus consider waves guided by two parallel conducting planes, in order to study the possibility of using nonlinear QED effects to balance diffraction and dispersion in directions parallel to the planes. We find that our system is governed by a $`1+2`$ dimensional nonlinear Schrödinger equation (NLSE). Using approximate variational methods, which have proven to be useful in similar situations (see e.g. ), we find that the nonlinearities can counteract diffraction and dispersion to form 2-dimensional light bullet solutions, i.e., solitary solutions to the NLSE which preserve the envelopes of the fields. The formation of light bullets was in other contexts (optical media) first studied by Silberberg .

## 2 Light bullet solutions

In a medium with polarization $`𝐏`$ and magnetization $`𝐌`$, the general wave equations for $`𝐄`$ and $`𝐁`$ are

$$\frac{1}{c^2}\frac{^2𝐄}{t^2}^2𝐄=\mu _0\left[\frac{^2𝐏}{t^2}+c^2(𝐏)+\frac{}{t}(\times 𝐌)\right]$$
(2)
and

$$\frac{1}{c^2}\frac{^2𝐁}{t^2}^2𝐁=\mu _0\left[\times (\times 𝐌)+\frac{}{t}(\times 𝐏)\right].$$
(3)
Furthermore, the effective polarization and magnetization in vacuum due to photon–photon scattering induced by the exchange of virtual electron–positron pairs are given by (see, e.g., Ref. )

$$𝐏=2\xi \left[2(E^2c^2B^2)𝐄+7c^2(𝐄𝐁)𝐁\right]$$

and

$$𝐌=2c^2\xi \left[2(E^2c^2B^2)𝐁+7(𝐄𝐁)𝐄\right].$$

We consider propagation between two parallel conducting planes with spacing $`x_0`$ (i.e., the region $`0xx_0`$ is vacuum surrounded by the plates that, as a starting point, are assumed to be perfectly conducting). We assume that only one $`\mathrm{TE}_{n0}`$ mode ($`n=1,2,\mathrm{}`$) is present. To linear order, this gives the fields

$`B_z`$ $`=`$ $`{\displaystyle \frac{\pi }{x_0}}\stackrel{~}{A}\mathrm{cos}\left({\displaystyle \frac{n\pi x}{x_0}}\right)\mathrm{exp}[\mathrm{i}(kz\omega t)]+\mathrm{c}.\mathrm{c}.,`$ (4a)
$`B_x`$ $`=`$ $`\mathrm{i}k\stackrel{~}{A}\mathrm{sin}\left({\displaystyle \frac{n\pi x}{x_0}}\right)\mathrm{exp}[\mathrm{i}(kz\omega t)]+\mathrm{c}.\mathrm{c}.,`$ (4b)
$`E_y`$ $`=`$ $`\mathrm{i}\omega \stackrel{~}{A}\mathrm{sin}\left({\displaystyle \frac{n\pi x}{x_0}}\right)\mathrm{exp}[\mathrm{i}(kz\omega t)]+\mathrm{c}.\mathrm{c}.,`$ (4c)
together with $`0\omega ^2/c^2k^2n^2\pi ^2/x_0^2`$. Here c.c. stands for complex conjugate.

It should be noted that the beams propagate with a 90 degree angle with respect to each other for the selfconsistent solutions obtained in a similar situation in Ref. . This would correspond to the relation $`k=n\pi /x_0`$ in our case. However, we have no such limitation in our derivation. The result of a 90 degree beam angle of Ref.  comes from requiring an explicit symmetry between the wave equations for $`𝐄`$ and $`𝐁`$. A similar requirement is not needed in our case, since the $`\mathrm{TE}_{n0}`$-modes introduced by the conducting planes do not have a complete electric-magnetic symmetry as, for example, the electric fields vanish close to the planes. What one still might worry about is whether or not both Faraday’s and Ampere’s law, which are equivalent to the two wave equations for $`𝐄`$ and $`𝐁`$, are properly solved in our case. However, the equation that will be derived for the vector potential is equivalent to Ampere’s law, which is equivalent to the wave equation for the electric field, and using the vector potential we therefore automatically solve Faraday’s law. It thus follows that our solution for the vector potential corresponds to electromagnetic fields that are consistent with the two wave equations for $`𝐄`$ and $`𝐁`$.

We have expressed the fields in terms of the vector potential amplitude $`\stackrel{~}{A}`$, given by $`𝐀=(0,A,0)`$, where

$$A=\stackrel{~}{A}\mathrm{sin}\left(\frac{n\pi x}{x_0}\right)\mathrm{exp}[\mathrm{i}(kz\omega t)]+\mathrm{c}.\mathrm{c}.$$
(5)
using the radiation gauge ($`\varphi =0`$). A nonlinear dispersion relation can be derived by inserting the linear expression for the fields into the right hand side of Eq. (3), taking the $`z`$-component of the resulting equation, and separating into orthogonal trigonometric functions. From that equation the coefficients in the NLSE can be found. However, here we will present an alternative derivation, starting directly from the Euler-Heisenberg Lagrangian (1), which is more elegant and gives the same result. Naturally we must first express $`_{\mathrm{EH}}`$ in terms of the electromagnetic potential, since it is the basis for the original variational principle. We take $`\stackrel{~}{A}=\stackrel{~}{A}(t,z)`$, and assume $`\stackrel{~}{A}`$ to be weakly modulated ($`|\stackrel{~}{A}/t||\omega \stackrel{~}{A}|`$, $`|\stackrel{~}{A}/z||k\stackrel{~}{A}|)`$, but we omit the slow $`y`$-dependence, since such a dependence makes the derivation technically more complicated without adding extra understanding. A diffraction term in the $`y`$-direction is lost by this procedure, but this effect is trivially added afterwards. To lowest order, we thus omit the nonlinear terms and the slow derivatives in $``$. Noting that $`=(1/2)[(A/t)^2c^2(A)^2]`$ and $`𝒢=0`$, we find that the non-oscillating part (i.e. the part that does not vanish after integration) is

$$_0/\epsilon _0=\frac{1}{2}(\omega ^2k^2c^2)|\stackrel{~}{A}|^2\mathrm{sin}^2\left(\frac{n\pi x}{x_0}\right)\frac{n^2\pi ^2c^2}{2x_0^2}|\stackrel{~}{A}|^2\mathrm{cos}^2\left(\frac{n\pi x}{x_0}\right).$$
(6)
After performing the integration over the region between the plates we find that the lowest order part of the action vanishes identically as the dispersion relation is here considered to be satisfied. Going to the next order of approximation in the Lagrangian we include first order slow derivatives, which yields $`_1/\epsilon _0=_0/\epsilon _0+\mathrm{i}\omega [(A/t)A^{}(A^{}/t)A]\mathrm{i}kc^2[(A/z)A^{}(A^{}/z)A]`$. After variation this leads to an equation where the envelope moves with the group velocity. The next order and final approximation includes second order slow derivatives and the $`^2`$ term in the Lagrangian. After performing the $`x`$-integration, dropping $`_0`$ (since the action is identically zero because of the dispersion relation) and eliminating the second order slow time-derivatives using $`^2/t^2v_g^2^2/z^2`$, the final expression for the Lagrangian is

$`_2`$ $`=`$ $`\mathrm{i}\omega \epsilon _0\left({\displaystyle \frac{\stackrel{~}{A}}{t}}\stackrel{~}{A}^{}{\displaystyle \frac{\stackrel{~}{A}^{}}{t}}\stackrel{~}{A}\right)\mathrm{i}kc^2\epsilon _0\left({\displaystyle \frac{\stackrel{~}{A}}{z}}\stackrel{~}{A}^{}{\displaystyle \frac{\stackrel{~}{A}^{}}{z}}\stackrel{~}{A}\right)`$ (7)
$`+(c^2v_g^2)\epsilon _0\left|{\displaystyle \frac{\stackrel{~}{A}}{z}}\right|^2+{\displaystyle \frac{3n^4c^4\pi ^4\xi }{x_0^4}}|\stackrel{~}{A}|^4,`$
where $``$ stands for integration over $`x`$. Variation with respect to $`\stackrel{~}{A}^{}`$ leads to

$$2\mathrm{i}\omega \left(\frac{}{t}+v_g\frac{}{z}\right)A+c^2\frac{^2A}{y^2}+\omega v_g^{}\frac{^2A}{z^2}+\frac{6n^4c^4\pi ^4\xi }{\epsilon _0x_0^4}|A|^2A=0,$$
(8)
where $`v_g`$ is the group velocity and $`v_g^{}`$ the group dispersion that follows from the linear dispersion relation. We have also added the diffraction term in the $`y`$-direction corresponding to the full amplitude dependence $`A=\stackrel{~}{A}(t,y,z)`$. Changing to a system moving with the group velocity while rescaling the coordinates and the amplitude according to

$$\begin{array}{cc}\tau =\omega t/2,\hfill & \upsilon =\omega y/c,\hfill \\ \xi =\sqrt{\omega /v_g^{}}(zv_gt),\hfill & a=\sqrt{\frac{6n^4c^4\pi ^4\xi }{\omega \epsilon _0x_0^4}}A,\hfill \end{array}$$
(9)
we obtain

$$\mathrm{i}\frac{a}{\tau }+\frac{^2a}{\upsilon ^2}+\frac{^2a}{\xi ^2}+|a|^2a=0.$$
(10)
which corresponds to the rescaled Lagrangian density

$$=\frac{\mathrm{i}}{2}\left(a^{}\frac{a}{\tau }a\frac{a^{}}{\tau }\right)\left|\frac{a}{\upsilon }\right|^2\left|\frac{a}{\xi }\right|^2+\frac{1}{2}|a|^4,$$
(11)
where the diffraction term in the $`y`$-direction is also included.

From now on we will look for cylindrically symmetric solutions of Eq. (10), i.e., $`a=a(t,\rho )`$, where $`\rho ^2=\upsilon ^2+\xi ^2`$. Equation (10) can then be written

$$\mathrm{i}\frac{a}{\tau }+\frac{1}{\rho }\frac{}{\rho }\left(\rho \frac{a}{\rho }\right)+|a|^2a=0,$$
(12)
while the Lagrangian density (11) takes the form

$$=\frac{\mathrm{i}}{2}\left(a^{}\frac{a}{\tau }a\frac{a^{}}{\tau }\right)\left|\frac{a}{\rho }\right|^2+\frac{1}{2}|a|^4,$$
(13)
with the action given by $`𝒜=\rho 𝑑\rho 𝑑\tau `$. Equation (12) is a 2-dimensional radially symmetric NLSE. Although exact solutions of this equation are not available, numerical and approximate analysis give a clear picture of the solutions. In particular, it has been shown that Eq. (12) allows a stationary solution where the diffractive/dispersive spreading of the pulse in the coordinate $`\rho `$ is balanced by the focusing effect created by the nonlinearity. An accurate analytical approximation of the dynamics of the pulse-like solutions of Eq. (12) can be obtained using direct variational methods involving the Lagrangian $``$ given by expression (13) and subsequent Rayleigh–Ritz optimization based on suitably chosen trial functions (see e.g.  and references therein). A convenient trial function for the present problem is

$$a_T(\tau ,\rho )=F(\tau )\mathrm{sech}\left[\frac{\rho }{f(\tau )}\right]\mathrm{exp}\left[\mathrm{i}b(\tau )\rho ^2\right]$$
(14)
which involves a complex amplitude $`F(\tau )`$, a pulse width $`f(\tau )`$, and a quadratic phase function modeling the phase front curvature. Inserting this ansatz into the variational integral and integrating over $`\rho `$, a reduced variational problem is obtained for the parameter functions, $`F(\tau )`$, $`F^{}(\tau )`$, $`f(\tau )`$ and $`b(\tau )`$. The subsequent reduced Euler–Lagrange equations can be rearranged to give $`F`$, $`F^{}`$ and $`b`$ as explicit functions of the width $`f(\tau )`$, which satisfies the equation

$$\frac{d^2f}{d\tau ^2}=\gamma \left(1\frac{I}{I_c}\right)\frac{1}{f^3}$$
(15)
where $`\gamma =4(\mathrm{ln}4+1)/(27\zeta (3))0.29`$, $`I(\tau )=f^2(\tau )|F(\tau )|^2=f^2(0)|F(0)|^2=I_0`$, and $`I_c=(2\mathrm{ln}2+1)/(4\mathrm{ln}21)1.35`$.111Here $`\zeta (p)`$ is the Riemann zeta function. Obviously, stationary solutions exist when the pulse power satisfies $`I_0=I_c`$. Furthermore, the full solution of Eq. (15) is

$$f(\tau )=f(0)\sqrt{1+\frac{\gamma }{f^2(0)}\left(1\frac{I_0}{I_c}\right)\tau ^2},$$
(16)
which shows that the stationary solution is unstable and either collapses to zero width in a finite time when $`I_0>I_c`$, or diffracts monotonously towards infinite width when $`I_0<I_c`$.

## 3 Discussion and conclusion

From Sec. 2, it is clear that the most interesting alternative is $`I_0>I_c`$, in which case the QED vacuum nonlinearities play a crucial role, leading to a collapse. In dimensional units the inequality $`I_0>I_c`$ roughly leads to $`E_{\mathrm{init}}^2k^2r_{\mathrm{init}}^2>E_{\mathrm{crit}}^2`$, where $`E_{\mathrm{init}}`$ and $`r_{\mathrm{init}}`$ are the initial electric field and beam radius, respectively, and the critical QED electric field is defined by $`E_{\mathrm{crit}}^2=\epsilon _0/\xi `$. Naturally this unbounded self-focusing will eventually be saturated by some kind of higher order nonlinear mechanism. However, this will not occur before an electric field level of order $`EE_{\mathrm{crit}}`$ is reached, which is the field strength when both our perturbative nonlinear calculation scheme and our starting expression, the Euler–Heisenberg Lagrangian, breaks down. For such extreme energy densities, higher order Feynman diagrams must be included in the QED description, and possibly the corresponding physical effects may then counteract the collapse scenario, resulting in a saturated beam radius $`r_{\mathrm{sat}}r_{\mathrm{init}}E_{\mathrm{init}}/E_{\mathrm{crit}}`$. The possibility to reach such extreme intensity levels is of very much interest. However, it is clear that setting up the conditions necessary for vacuum light bullets to be formed, is a technological challenge.

Presently the electric fields that can be supported by the walls before field emission takes place are of the order $`10^8`$ V/m . For such field strengths we note from the above estimates that we must have $`kr_{\mathrm{init}}10^{11}`$, in which case either the initial beam radius becomes unrealistically large (of the order of km), or the wavelength to short (i.e. in the short uv or x-ray regime) for the walls to be conducting. Two things should be stressed, however. Firstly it should be noted that the experimentally possible field level before field emission takes place have increased significantly in recent years . Secondly, in our geometry there is no normal component of the electric field at the wall surfaces, which may allow for a higher (central) field strength than normal, relaxing the necessary value on $`kr_{\mathrm{init}}`$. A slightly different but related technological question is whether the huge beam powers needed for vacuum self-focusing can be reached. This issue was discussed in Ref. , where it was concluded that it may occur within the next 10–15 years, given the current rate of technological improvement.

A very interesting question from a principal point of view, would be whether it is possible to have fully three-dimensional QED-structures which do not require any guiding support. However, it seems that this issue cannot be addressed within a perturbational approach, and further research is thus required.

In this paper we have shown that the photon–photon scattering, due to the exchange of virtual electron–positron pairs, that effectively gives rise to self-interaction terms in Maxwell’s equations, can cause optical collapse in vacuum if the right conditions are met. More specifically, high intensity electromagnetic waves guided by two parallel conducting planes can form light bullets, which can collapse if the intensity of the beams are high enough. It is possible that higher order corrections to the Euler–Heisenberg Lagrangian can counteract such a collapse scenario, thus leading to stable optical vacuum structures. Further research is also necessary to find out more about how the optimum configurations can be changed, both in the present paper and in previous work (e.g. ).


text:
# 1 Introduction

## 1 Introduction

Radio–loud active galactic nuclei, galactic sources showing apparent superluminal motions (GS), and gamma–ray bursts (GRBs) are the three classes of sources where there is evidence for large quantity of matter moving very close to the speed of light. It is natural to wonder if these sources work in a similar way, namely if the machine able to produce, collimate and accelerate matter to relativistic velocities is the same. Furthermore one would like to know if in all three systems the radiation we see comes from the transformation of bulk kinetic (or magnetic) into random energy and, if so, if the mechanisms able to do that are similar. If this turns out to be true, then the galactic superluminal sources could be used as a sort of laboratory for the study of radio–loud AGNs, both because their vicinity (i.e. high fluxes) allows more detailed studies, but also because the ratio in black hole masses is probably $`10^8`$: if all timescales scale with the black hole mass, a millisecond in the life of GRS 1915+105 is equivalent to one day in the life of a blazars (see the recent review by Blandford 2002).

GRBs, on the other hand, are the most powerful events occurring after the Big Bang. But this does not mean that we are witnessing a single tremendous explosion. Indeed, GRBs have durations lasting up to hundreds of seconds, and recent modeling associates this time with the duration of the central engine powering them. Since a stellar mass black hole is involved, we have that a GRB of 100 s duration lasts for $``$a million of dynamical times. Put in this way, the “explosion” is instead a quasi steady–state process!

Thus the comparison among the three systems and their central engines is certainly legitimate, hopefully leading to a deeper understanding of how the formation, launch and acceleration of jet work.

## 2 The bulk Lorentz factor

Blazars — The best evidences for bulk relativistic motion in blazars come from VLBI observations of knots of radio emission moving superluminally. Apparent speeds up to $`\beta _{\mathrm{app}}=30h^1`$ are measured (Jorstad et al. 2001, see also these proceedings), at least in those blazars that have been discovered by EGRET to be powerful $`\gamma `$–ray emitters. With $`H_0=65`$ Mpc km-1 s-1, we have objects moving with $`\beta _{\mathrm{app}}40`$, which indicates even larger Lorentz factors $`\mathrm{\Gamma }`$ (note that $`\mathrm{\Gamma }`$ must always be larger than $`\beta _{\mathrm{app}}`$). These values are larger than the “fiducial” $`\mathrm{\Gamma }10`$ we are used to, and can help solving the problem of the large radio brightness temperature of intraday variable sources (see Wagner & Witzel 1995 for a review), once corrected for interstellar scintillation (sources scintillate if their angular dimensions are small, hence the limit on the brightness temperature; Kedziora–Chudczer et al. 2001). Indications of a large degree of beaming come also from spectral fitting, especially of low powerful, TeV emitting, blazars where beaming factors $`\delta >20`$ are derived (Tavecchio et al. 1998). There are also indications that the jet could be structured, with a fast “spine” surrounded by a slower “layer” (Laing 1993, Chiaberge et al. 2001), explaining for example why the non–thermal radiation of the core of radio–galaxies is not as faint as predicted if the plasma is moving at the same large speed and if they are observed at large angles with respect to the jet axis.

Galactic superluminals — Besides SS 433 ($`\beta =0.26`$), GRS 1915+105 and GRO J1655–40 (both of them have $`\beta 0.9`$ and $`\mathrm{\Gamma }`$2.5) there are other recently discovered sources showing radio jets with moving features (see the reviews by Mirabel & Rodriguez 1999 and by Fender 2001). Among these, SAX J1819.3–2525 seems particularly interesting, since a new determination of its distance put this source at more than 8 kpc, much further away than previously thought (Orosz et al. 2001). If one assumes that the jets were ejected during one of the bright X–ray flares one then finds $`\mathrm{\Gamma }10`$.

Gamma–Ray Bursts — We do not have direct measurement of $`\mathrm{\Gamma }`$ during the prompt emission, so we must rely on theory to estimate it: the “fiducial” value is $`\mathrm{\Gamma }`$100–300. It cannot be much less than 100 to explain the very fast (millisecond) variability observed in the $`\gamma `$–ray light curve of their prompt emission. It cannot be much larger than that since otherwise we would have that the source (the fireball) becomes transparent to Thomson scattering before the acceleration phase ends, letting the internal radiation escape with a quasi black body spectrum (which is not observed, see Piran 1999 for a review). A similar constraint ($`\mathrm{\Gamma }>100`$) comes also from associating the emission above 100 MeV seen in a few cases with the prompt emission (see Fishman & Meegan 1995 for a review), thus imposing that the source must be transparent for the $`\gamma \gamma e^+e^{}`$ process. Similarly to blazars, there is the possibility that also the “jets” of GRBs are structured, with the bulk Lorentz factor decreasing for increasing angles from the jet axis (Rossi, Lazzati & Rees, 2002).

## 3 The power of jets

Blazars — Without any doubt, the knowledge of the power and the energy carried by jets is the prime parameter to start any modeling. Despite this, the power of blazar (and radiogalaxy) jets is poorly known. This is due to two main reasons: 1) the radiation produced by relativistic jets is highly beamed, and the flux we receive is therefore highly enhanced (in blazars) or dimmed (in radiogalaxies) by relativistic effects; 2) we still do not know for sure the matter content of jets, i.e. we are still debating whether if most of the matter, is made by electron–positron pairs or by a normal electron–proton plasma. On the other hand, we know that the radio lobes of radio–galaxies and blazars are a sort of calorimeter, since the cooling time are long on these scales, and minimum energy considerations allow to estimate lower limits on the total energetics. Yet again, even in this case we do not know either the contribution of protons to the total energy or the emitting plasma filling factor (see however recent results due to Chandra observations, Fabian et al. 2002). Bearing in mind these limitations, the total energy of a radio lobe, divided by its lifetime (estimated from spectral aging or from advance motion) allows to estimate the average power required by the radio lobes to emit and expand. This is the average power a jet must have. This estimate has been done, among others, by Rawling & Saunders (1991): they find an average power ranging from $`10^{43}`$$`10^{44}`$ erg s-1 for FR I radiogalaxies to $`10^{46}`$$`10^{47}`$ erg s-1 for FR II radiogalaxies and radio–loud quasars.

One can also calculate the power carried by the jet by inferring its density through modeling the observed SED and requiring that the jet carries at least the particles and the magnetic field necessary to make the radiation we see. This has been done on the pc scale by Celotti & Fabian (1993), on sub–pc scale (the $`\gamma `$–ray emitting zone) by Celotti & Ghisellini (2002, see also Ghisellini & Celotti 2002), and on the hundreds of kpc scale (the X–ray jets seen by Chandra) by Celotti, Ghisellini & Chiaberge (2001) and Tavecchio et al. (2000). These studies suggest large values of the power transported by the jet and require the presence of a dynamically dominating proton component (see also arguments by Sikora & Madejski 2000).

Galactic superluminals — The jet of Galactic superluminal sources is not a steady feature, and the conditions for the sporadic launching of jets in GRS 1915+105 and GRO J1655–40 are still under debate. But during the major flare events in GRS 1915+105 we can calculate the jet power by knowing the bulk Lorentz factor, the viewing angle, and assuming that the jet is at least carrying the particles and magnetic field that account for the synchrotron emission we see. Therefore, besides the bulk motion of the emitting particles, one has to account also for some of the jet power being in the form of a Poynting vector. We therefore have a powerful tool to find the minimum jet power, which corresponds to rough equipartition between the magnetic field and the particle density in the comoving frame. In this way Gliozzi, Bodo & Ghisellini (1999) estimated a minimum jet power in GRS 1915+105 ranging between $`10^{39}`$ erg s-1 (pure electron–positron jet) and $`10^{40}`$ erg s-1 (proton–electron jet).

Gamma–Ray Bursts — Are GRB collimated? If their emission were isotropic, then some of them would emit more than $`10^{54}`$ erg only in $`\gamma `$–rays. This is equivalent to one solar mass completely transformed into energy. One should not be particularly impressed by that, since the rotational extractable energy of black hole is 29 per cent of its mass, and therefore reaches $`5\times 10^{54}`$ erg for a 10 $`M_{}`$ black hole. On the other hand we do have some indications from afterglows that most of the emission is collimated in a cone of semiaperture $`\theta _\mathrm{j}`$, and that therefore the true radiated energy is a factor $`(1\mathrm{cos}\theta _\mathrm{j})`$ less than the isotropic equivalent (for two jets). The precise amount of collimation is uncertain (Frail et al. 2001, Ghisellini et al. 2002), and we do not know yet the value for the efficiency in transforming bulk energy into radiation. Most energy is however contained in the first prompt phase, and not in the much longer afterglow. Note that in the case of GRBs the word “jet” is somewhat misleading, since the radial extension of the “jet” is of order of hundreds of light seconds. This is the reason why Piran proposed the term ”flying pancakes”. Bearing the above uncertainties in mind, the energetics of GRBs is in the range $`10^{51}`$$`10^{53}`$ ergs.

### 3.1 Comparison with accretion disk power

Blazars — The luminosity produced by the accretion disk in blazars $`L_\mathrm{d}`$ can be derived “directly” only when the blue bump is visible, i.e. when it is not completely swamped by the beamed non–thermal continuum or by the optical emission of the host galaxy. In many (but not all) cases we can however estimate it using the emission lines, both narrow (e.g. Rawlings & Saunders 1991) and broad (Celotti, Padovani & Ghisellini 1997). In completely featureless BL Lacs, instead, we can derive only upper limits for the luminosity of the disk. The most remarkable fact of these studies is that there is a rough equality between the average power of the jet and the power in the accretion disk, once an average covering factor of $`10^2`$ and $`10^1`$ is assumed for the narrow and broad line clouds, respectively.

Galactic superluminals — Here the disk luminosity is emitted at relatively soft X–ray energies, where photoelectric absorption can be an issue (possibly making SS 433 not so bright in X–rays). For GRS 1915+105 the X–ray luminosity, during flares, can be $`10^{39}`$ erg s-1 (see e.g. Mirabel & Rodriguez 1999). The recent determination of the mass of its black hole (i.e. 14$`\pm 4`$ solar masses, Greiner 2001) makes the accretion disk of GRS 1915+105 Eddington limited. The ratio $`L_\mathrm{j}/L_\mathrm{d}`$ of the jet power to accretion luminosity is then of the order of 1–10 during major ejection events.

Gamma–Ray Bursts — In GRBs the power due to accretion is completely unobservable. We must rely on theory and models to infer it. There is a growing consensus that long GBRs (i.e. GRBs whose prompt emission lasts for more than 2 seconds: these are the only ones for which the good locations allowed the follow–up observations in other bands) are associated with the collapse of stars more massive than ordinary supernovae (and therefore called hypernovae), or a two–step collapse forming first a netron star and then a black hole (the SupraNova scenario of Vietri & Stella 1998). In any case the final resulting compact object is a fast spinning black hole surrounded by a very dense torus of mass $`M_\mathrm{T}`$0.1–0.2$`M_{}`$. It is possible that the duration of the GRB prompt emission is associated with the accretion time. If this is the case, the accretion rate is obviously huge even if the radiation that can escape is a tiny fraction, since the accreting material is completely opaque (scattering optical depths of order $`10^{12}`$ or so, see Table 1).

## 4 Outflowing vs accreting mass rates

For powerful blazars and Galactic Superluminals we can write:

$$L_\mathrm{d}=\eta \dot{M}_{\mathrm{in}}c^2$$
(1)
$$L_\mathrm{j}=\mathrm{\Gamma }\dot{M}_{\mathrm{out}}c^2$$
(2)
and therefore

$$\mathrm{Blazars},\mathrm{GS}\frac{\dot{M}_{\mathrm{out}}}{\dot{M}_{\mathrm{in}}}=\frac{\eta }{\mathrm{\Gamma }}\frac{L_\mathrm{j}}{L_\mathrm{d}}=\mathrm{\hspace{0.17em}10}^2\frac{\eta _1}{\mathrm{\Gamma }_1}\frac{L_\mathrm{j}}{L_\mathrm{d}}$$
(3)
where $`\eta `$ is the efficiency in converting the accreted mass into energy, $`\dot{M}_{\mathrm{in}}`$ and $`\dot{M}_{\mathrm{out}}`$ are the accretion and the mass outflowing rates, respectively. We use the notation $`Q=10^xQ_x`$.

As mentioned, in GRBs we cannot estimate the mass accretion rate directly. Let us assume that the accretion process lasts for the duration of the burst, and that the total accreted matter corresponds to the mass of the torus surrounding the black hole, $`M_\mathrm{T}=0.1M_{\mathrm{T},1}M_{}`$. In this case

$$\dot{M}_{\mathrm{in}}=\mathrm{\hspace{0.17em}2}\times 10^{32}\frac{M_{\mathrm{T},1}}{t_{\mathrm{burst}}}\mathrm{g}\mathrm{s}^1$$
(4)
yielding

$$\mathrm{GRB}\frac{\dot{M}_{\mathrm{out}}}{\dot{M}_{\mathrm{in}}}=\frac{L_\mathrm{j}t_{\mathrm{burst}}}{\mathrm{\Gamma }M_\mathrm{T}c^2}=\mathrm{\hspace{0.17em}5}\times 10^3\frac{E_{\mathrm{burst},52}}{\mathrm{\Gamma }_2M_{\mathrm{T},1}}.$$
(5)
We then obtain the remarkable result that (at least in this model) the ratio of outflowing and accreting mass is roughly the same in blazars, Galactic Superluminals and Gamma Ray bursts. If true, this result could explain why GRBs fireballs have that particular barion loading, which allows them to be at the same time light enough to be relativistic (contrary to supernovae), but heavy enough to convert all internal energy into bulk motion (requiring the fireball not to become transparent too early for scattering). Alternatively, one can assume that there is a typical (and fixed) barion loading in jets, thus deriving the resulting bulk Lorentz factor.

## 5 The central engine

There is some consensus (although certainly not unanimous) that the powering mechanism of jets is the extraction of rotational energy of a spinning black hole, through the Blandford & Znajek (1977) process. As an order of magnitude estimate, the power which can be extracted by this mechanism is (Blandford & Znajek 1977; Rees et al. 1982)

$`\mathrm{Blazars}`$ $`L_{BZ}`$ $`\mathrm{\hspace{0.17em}6}\times 10^{46}\left({\displaystyle \frac{a}{m}}\right)^2M_9^2B_4^2\mathrm{erg}\mathrm{s}^1`$
$`\mathrm{GRBs}`$ $`L_{BZ}`$ $`\mathrm{\hspace{0.17em}6}\times 10^{50}\left({\displaystyle \frac{a}{m}}\right)^2M_1^2B_{14}^2\mathrm{erg}\mathrm{s}^1`$
$`\mathrm{GS}`$ $`L_{BZ}`$ $`\mathrm{\hspace{0.17em}6}\times 10^{38}\left({\displaystyle \frac{a}{m}}\right)^2M_1^2B_8^2\mathrm{erg}\mathrm{s}^1`$ (6)
where $`(a/m)`$ is the specific black hole angular momentum ($`1`$ for maximally rotating black holes).

In the following we will compare the $`B`$–values required for the Blandford & Znajek mechanism to work with the gravitational energy density of the matter at $`R_s`$. To do that, we will assume that the accretion process, in blazars and in Galactic Superluminals, is converting in radiation a fraction $`\eta =0.1\eta _1`$ of the total accreted mass–energy: $`L_\mathrm{d}=\eta \dot{M}_{\mathrm{in}}c^2`$. This allows to estimate $`\dot{M}_{\mathrm{in}}`$ from the blue bump luminosity, or, when this is not visible, through the luminosity of the emission lines, (assuming they reprocess a fraction of $`L_\mathrm{d}`$ of the same order of the covering factor of the emission line clouds, i.e. 10% for the broad and 1% for the narrow emission lines, respectively).

For GRBs, instead, the mass accretion rate can be derived assuming, as a working hypothesis, that their duration is controlled by the time needed to accrete the mass of the dense torus surrounding the spinning black hole. In this way $`\dot{M}_{\mathrm{in}}M_\mathrm{T}/t_{\mathrm{burst}}`$.

To find the corresponding density $`n_\mathrm{p}`$ of the matter we use

$$\dot{M}_{\mathrm{in}}=\mathrm{\hspace{0.17em}2}\pi Rh\beta _\mathrm{r}n_\mathrm{p}m_\mathrm{p}c$$
(7)
where $`\beta _\mathrm{r}c`$ is the radial inflow velocity, and $`h`$ is the height of the disk at the radius $`R`$. The gravitational energy density $``$, close to the Schwartzschild radius is

$$n_\mathrm{p}m_\mathrm{p}c^2=\frac{\dot{M}_{\mathrm{in}}c}{2\pi (Rh/R_\mathrm{s}^2)R_\mathrm{s}^2\beta _\mathrm{r}}$$
(8)
Equipartition of the magnetic energy density \[$`B_{\mathrm{eq}}=(8\pi )^{1/2}`$\] gives

$`\mathrm{Blazars},\mathrm{GS}B_{\mathrm{eq}}`$ $`=`$ $`\left[{\displaystyle \frac{4L_\mathrm{d}}{\eta c(Rh/R_\mathrm{s}^2)R_\mathrm{s}^2\beta _\mathrm{r}}}\right]^{1/2}`$
$`\mathrm{GRBs}B_{\mathrm{eq}}`$ $`=`$ $`\left[{\displaystyle \frac{4M_\mathrm{T}c}{t_{\mathrm{burst}}(Rh/R_\mathrm{s}^2)R_\mathrm{s}^2\beta _\mathrm{r}}}\right]^{1/2}`$ (9)
the fiducial values of $`B_{\mathrm{eq}}`$ are listed in Table 1 for the three classes of sources. One can see that the derived values are of the order required by the Blandford–Znajek process to provide the right power.

This can be taken as promising, since the found values match the required ones: note that it requires that the magnetic field energy density is of order of the gravitational energy density, which, depending on the accretion conditions, can be much larger than either the gas pressure of the accreting matter or the radiation energy density produced by accretion. This latter possibility appears to occur for GRBs.

In Table 1 we also list the value of the scattering optical depth, calculated in the vertical direction:

$$\tau _\mathrm{T}\sigma _\mathrm{T}n_\mathrm{p}h.$$
(10)
We can see that in blazars it is of order unity, a few in (flaring) Galactic superluminals and huge in GRBs.

## 6 Internal shocks

In 1978, Rees proposed that the jet of M87 could be powered by collisions among different part of the jet itself, moving at different speeds. When colliding, they would produce shocks, giving rise to the non–thermal radiation we see. Although born in the AGN field, this idea of internal shocks grew up more robustly in the GRB field, to become the “paradigm” to explain their prompt emission (see e.g. Rees & Mészáros 1994). Faster and later shells can then catch up slower earlier ones, dissipating part of their bulk kinetic energy into radiation. However, all shells are and remain relativistic: after the collision the merged shells move with a bulk Lorentz factor which is intermediate between the two initial ones. It is then clear that this mechanism has a limited efficiency, because only a small fraction of the bulk energy can be converted into radiation (unless the contrast between the two initial Lorentz factor is huge, see Beloborodov 2000, but also Ghisellini 2002).

In blazars, on the other hand, we require a small efficiency, since most of the bulk energy has to go un–dissipated to the outer radio–lobes. This model is very promising, since it can explain some basic properties of blazars:

* About 10 per cent of the power in bulk motion is dissipated into radiation, while 90 per cent is available to power the radio lobes.
* Most of the dissipation occurs at a few hundreds of Schwarzschild radii, far enough from the black hole and accretion disk to avoid relevant $`\gamma `$$`\gamma `$ $``$ $`e^\pm `$ absorption of $`\gamma `$–rays, but still within (in powerful objects) the broad line region, which provides seed photons for the Compton scattering process (leading to a dominant high energy emission).
* The shells that have already collided once can collide again, at larger distances from the black hole with reduced efficiency, since the $`\mathrm{\Gamma }`$–contrast of the colliding shells decreases. This can explain why the luminosity of jets decreases (but it does not vanish) with distance. Fig. 1 shows the efficiency (i.e. the fraction of the luminosity emitted to the jet carried power) as a function of distance for a simulation we have done for Mkn 421 (Guetta et al. 2002).
* Variability is a built–in feature of the model. See Zhang et al. (2002, see also this volume), for an example of blazar variability that might be explained by internal shocks.
* There is the possibility (not yet studied in detail) of a link between the flares at optical and $`\gamma `$–ray energies and the flares in the radio–mm band (see Spada et al. 2001 for an example of correlation between the mm and the optical–$`\gamma `$–ray fluxes in 3C 279).

A detail study via numerical simulations of the predictions of this model in blazars has been carried out by Spada et al. (2001) for powerful objects (like 3C 279). We are now studying less powerful objects (like Mkn 421 and BL Lac itself, Guetta et al. 2002). As discussed by Kaiser, Sunyaev & Spruit (2000), internal shocks could also work for Galactic Superluminals, although detailed simulations have yet to be done.

### 6.1 Internal shocks and the blazar sequence

One of the important characteristics of the internal shock model for blazars is that the injection of energy in each collision is finite in time. The generation of relativistic particles lasts for the collision time, which is of the order of the time needed for one shell to overtake the other, i.e. $`t_{\mathrm{inj}}^{}\mathrm{\Delta }R^{}/c`$ in the comoving frame. The energy distribution of the relativistic particles responsible for the emission, in this case, never reaches a steady state, and the maximum emitted flux corresponds to the end of the injection (if the energy distribution of particles, $`N(\gamma )`$, is steeper than $`N(\gamma )\gamma ^2`$). This can explain why there is a sequence in the SED of blazars: in low power sources lacking emission lines the cooling time is long, and only the most energetic particles can cool in $`t_{\mathrm{inj}}`$. Suppose to inject a power law of relativistic particles, between $`\gamma _{\mathrm{min}}`$ and $`\gamma _{\mathrm{max}}`$, with slope $`s`$ \[i.e. $`Q(\gamma )\gamma ^s`$\]. After the injection phase (i.e. after $`t_{\mathrm{inj}}`$), only those particles for which $`t_{\mathrm{cool}}(\gamma )<t_{\mathrm{inj}}`$ have cooled, producing a break in the particle distribution at $`\gamma _\mathrm{c}>\gamma _{\mathrm{min}}`$. Below $`\gamma _\mathrm{c}`$ the particle distribution retains its original slope \[$`N(\gamma )\gamma ^s`$\], while above $`\gamma _\mathrm{c}`$ the distribution has a slope steeper by one $`N(\gamma )\gamma ^{(s+1)}`$. If $`2<s<3`$ the peaks of the synchrotron and the inverse Compton spectra are produced by electrons at $`\gamma _\mathrm{c}`$. For powerful sources, on the other hand, electrons of almost all energies can cool in a time $`t_{\mathrm{inj}}`$, due to the larger magnetic field in their jet, but especially due to the contribution of the broad line photons which enhances the inverse Compton cooling rate. In this case the particle distribution $`N(\gamma )\gamma ^{(s+1)}`$ above $`\gamma _{\mathrm{min}}`$ and $`N(\gamma )\gamma ^2`$ below.

Therefore sources will have different peak locations, according to the amount of cooling suffered during the time of the shell–shell collision. Calling $`\gamma _{\mathrm{peak}}`$ the random Lorentz factor of those electrons emitting most of the radiation (i.e. emitting at the peaks of the SED), we find a correlation between $`\gamma _{\mathrm{peak}}`$ and the amount of radiative plus magnetic energy density $`U`$ (as measured in the comoving frame) which has two branches: in powerful objects $`\gamma _{\mathrm{peak}}=\gamma _{\mathrm{min}}U^{0.5}`$, while in weak lineless BL Lacs we have $`\gamma _{\mathrm{peak}}=\gamma _\mathrm{c}U^1`$ (Ghisellini, Celotti & Costamante, 2002).

## 7 Conclusions

We have rather compelling indications that relativistic jets in different systems are produced via similar processes, and that also the way in which they dissipate part of their kinetic energy into radiation may be similar. The matter outflowing in relativistic jets may be of the order of 1 per cent of the matter inflowing into the black hole. If powered by spinning black holes, the value of the magnetic field close to the horizon must be of the order of the gravitational energy density of the accreting matter, which can be much more of the gas or radiation pressure.

There is now growing consensus that the massive black hole in radio–loud AGNs exceeds $`10^8M_{}`$ (see e.g. Lacy et al. 2001; Ghisellini & Celotti 2001). If maximally spinning, its rotational energy is more than enough to account for the energy contained in the radio–lobes of radio–loud sources. There may be a problem, however, in the mechanisms proposed so far to extract this rotational energy. They must be very efficient. Indeed, jets can transport more power than that radiated by accretion disks: this is obviously true in GRBs, but it is also the case for lineless BL Lac objects and Galactic Superluminal sources during major radio flares. In fact, as already known at low accretion rates, the accretion disk luminosity could be lowered by the de–coupling of proton and electrons (Rees et al. 1982, Narayan, Garcia & McClintock 1997), while, at high accretion rates, the optical depth of the accretion disk is large enough to trap photons (Begelman 1979). Jets seem not to be bound by these limitations.

## Acknowledgements

We are grateful to ASI and the Italian MIUR (AC) for support.