text:
# Linear Programming helps solving large multi-unit combinatorial auctions

## 1 Background, Introduction and Previous Work

The last decade has seen extensive denationalizations and a fundamental change in economic patterns brought about by the use of the Internet as a world-wide market place. Both phenomena promoted the use of auctions to a central place among economic mechanisms. The study of complex auctions is now a most active area of research, drawing from the fields of Mechanism Design in Economics and from the Theory of Algorithms in Computer Science. Combinatorial auctions in which a large number of items are sold and in which bidders may express preferences over bundles, i.e., subsets, of items have drawn a lot of attention recently .

We consider here the “winner determination” problem, i.e., the algorithmic problem of finding the allocation that maximizes the declared social welfare. We consider both (single-unit) combinatorial auctions and multi-unit combinatorial auctions in which a number of identical items of each type are sold. In , this problem was shown to be NP-hard. In , it is shown that it is even hard to find a non-trivial approximation to the optimal allocation. Those results consider worst-case complexity and do not necessarily imply that, in practice, one cannot find the optimal solution of real-life large combinatorial auctions. Some of the papers above proposed winner-determination algorithms and tried to show they perform well even on large auctions. Others proposed mechanisms based on non-optimal allocations. A number of researchers proposed sets of problems on which to test candidate algorithms: most notably Leyton-Brown, Shoham and Tennenholtz , Leyton-Brown, Pearson and Shoham  and de Vries and Vohra .

We implemented a Branch-and-Bound search program along the lines of  for finding the optimal allocation in multi-unit combinatorial auctions. This program makes heavy use of a linear programming routine to bound from above the value of extensions of partial allocations. Common wisdom in the field has it that linear programming is too time-consuming to be used in practice, and one should use lighter, if if less accurate, methods to obtain upper-bounds. We report here results that go against this common wisdom: pruning is of the essence and one should not hesitate to devote all the time necessary to compute tight upper-bounds in order to prune aggressively. Such conclusions hold for all the distributions described above. We present, in Section 2.3, an important way of saving on those expensive calls to the LP routine. Section 2 presents the algorithm we used and the remainder of the paper discusses our experimental results.

## 2 Branch and Bound

### 2.1 Winner Determination in Multi-Unit Combinatorial Auctions

We assume $`n`$ different commodities and $`k_j`$ identical items of commodity $`j`$ for $`j=1,\mathrm{}n`$. The problem we are trying to solve is that of determining the optimal allocation of the items. A set of bids is given: each bid requests a number of units (possibly zero) from each commodity and offers a price $`p`$ for the whole set. A subset of the set of all bids is conflict-free if, for each commodity, the sum of the units requested does not surpass the number of units for sale. The problem is to find a conflict-free subset that maximizes the sum of the prices proposed.

### 2.2 Our Algorithm

The algorithm we experimented with is very similar to the one described in . We assume that none of the bids submitted requests more units than available, i.e., each singleton set of bids is conflict-free. Our algorithm consists of an initialization phase and a main (recursive) routine.

In the initialization phase, we use a host of fast heuristics to find a good allocation. In practice, this initialization phase is performed very fast and gives an allocation that is not much less than the optimal one. Many times, one, in fact, gets the optimal allocation, but one, obviously, does not know it is optimal. The heuristics used are all members of the greedy family: one chooses a bid $`b`$, grants it, subtracts the quantities requested by $`b`$ from the stock of available units, eliminates all bids that cannot be granted anymore, i.e., those that are not conflict-free, chooses a bid, and so on. Each heuristic provides a feasible allocation. One keeps the best one. The different heuristics we propose to use differ only in the way they choose the next bid to be granted. The best allocation found is kept.

The main (recursive) routine is given:

1. a partial allocation, i.e., a set of bids already granted,
2. a stock of units available,
3. a set of bids not yet granted, each of them requesting only a number of units not larger than the number of units in stock, and
4. it uses a global variable describing the best allocation found so far.

The problem is solved by calling the main routine with an empty partial allocation, the full stock of available units, the set of all bids, and after having initialized the best allocation to the best allocation found in the initialization phase.

The main routine is:

1. Stop if the third argument is empty, i.e., no more bids waiting, return,
2. Update if the value of the partial allocation is larger than that of the best solution found so far, update this best solution,
3. Bound bound from above the value of the optimal allocation of the units left in stock to the bids that are the third argument of the routine. This is done by solving the Linear Programming problem that is the fractional relaxation of the Integer Programming problem describing the optimal allocation,
4. Prune checks whether the search can be stopped, i.e., the subtree rooted at the current node can be pruned. The condition for pruning is that the sum of the values of the partial allocation and the solution of the Linear Programming problem just found be not greater than the value of the best allocation found so far. If the condition holds, return. If the condition does not hold, go on to the next step,
5. Choose choose a bid from the third argument,
6. Left Call call recursively the main routine with a partial allocation augmented by the bid just chosen, the stock of available units diminished by subtracting the units requested by the bid chosen and the list of bids diminished by subtracting the chosen bid and all bids that request more units than available in (the diminished) stock,
7. Right Call call recursively the main routine with the original partial allocation, the original stock and the list of bids without the chosen bid.

The Left Call solves the auction obtained after the chosen bid has been granted. The Right Call solves the auction obtained after the chosen bid has been denied.

The only step under-specified in this routine is Choose: we shall discuss at length how the next bid to branch on should be chosen in Section 6. The only time-consuming step in the main routine is the Linear Programming step in Bound. We have a way to save on those expensive calls to the LP routine.

### 2.3 Saving on calls to LP

Notice that, in Right Call, the problem to be solved is very similar to the original one: same partial allocation, same stock. This recursive call will trigger a call to the LP routine (call it LP2) solving a LP problem that is very similar to the one just solved (call it LP1). The only difference is that the bid just chosen has disappeared from the list of bids considered by LP2. If this bid did not enter the optimal (fractional) solution of LP1, i.e., if its coefficient in the solution of LP1 was equal to zero, then, the solution of LP2 is the same as that of LP1 and therefore one can save a call to the LP routine. One can easily see that, in this case, no pruning can take place. Figure 1 shows that the saving can be considerable.

The auctions are taken from . On the x-axis: the number of bids considered, on the y-axis: under log scale, the number of nodes visited or the number of calls to the LP routine.

## 3 Experimental Parameters

We ran all our experiments on a Pentium Pro 930 MHz with 512K of memory, running Windows 2000. Our programs were coded in Matlab (R12). We used the linprog facility of the Optimization ToolBox for Linear Programming: this is a primal-dual algorithm. All our results are averages over fifteen auctions drawn from a distribution to be specified for each graph. We did experience some memory problems (probably due to some remaining bugs in Matlab Release 12), but found Matlab extremely easy to use and its linprog facility excellent. The remainder of the paper presents, in graphic form, the results of our experiments. Some of those graphs are difficult to decipher in black and white: we apologize and will produce more readable graphs for the conference paper.

## 4 Experimental Results

### 4.1 Global Picture

Figure 2 describes the quality of the optimal solution, the lower bound obtained in the initialization phase and two upper bounds for multi-unit auctions. The upper curve (extended norm bound) describes the easy to compute upper bound proposed in . The second (from the top) curve describes the upper bound given by Linear Programming. The latter can be proved to be always at most equal to the former. On the x-axis: the number of bids considered (from 250 to 750) and the y-axis the average value (in monetary terms, numbers are not important) of the allocations. This graph shows two things. First, that the lower bound computed in the initialization phase and the upper bound computed by LP are very close. The optimal solution is, obviously, in-between. The larger the auction, the better the convergence of those three values. The second thing is that the easy-to-compute upper-bound of  is not even close to the LP bound and seems to get (relatively) worse as the number of bids increases.

The main claim of this paper is that using the time-consuming LP routine is preferable to the use of a lighter but easier to compute upper bound in Bound. Figure 3 compares the running times of the algorithm above that uses LP and the running time of the same algorithm using the easy-to-compute upper bound proposed by Leyton-Brown and al.

The distribution of auctions is taken from , for auctions of different sizes. On the x-axis: the number of bids. On the y-axis: the running time in seconds (log scale). By using LP, one gains two orders of magnitude for large auctions. The lighter upper bound is better for very small numbers of bidders.

## 5 Initialization Phase

In the initialization phase, in order to find a good solution to be used for pruning, we propose to use many greedy heuristics and keep the best solution found. Figure 4 describes the quality of 16 different greedy heuristics on auctions of varying size (100 to 1000 bids, always 10 goods) taken from the distribution of .

On the y-axis: the percentage of the value of the optimal solution reached by the approximation: 100 means perfect approximation, 50 means half of the optimal solution. One expects that criteria that pick up more promising bids first give better results. This is what happens. One sees three main groups. The middle group includes two criteria: picking at random and picking following the (in fact random) order with which bids were given. The lower (worst) group of criteria comprises criteria that choose unattractive bids first. The higher group contains all criteria we thought attractive. It contains both unnormalized and normalized criteria (see ). The square-root criterion:

$$r(a)=\frac{p}{\sqrt{_{j=1}^nq_j}},$$
(1)
that has been proved to be theoretically optimal, gives the best approximation. The average price per unit gives markedly inferior results. Two criteria based on Linear Programming give pretty good results. The first one solves the relaxed LP problem for the original auction and orders the bids according to their coefficient in the optimal solution, in descending order. The second, improved one, is an adaptative version of the first: it solves a LP problem before choosing each bid (after having eliminated all bids that cannot be granted in full). More on those criteria is found in Section 6.

## 6 Choosing the bid on which to branch

The main factor in determining the efficiency of our Branch-and-Bound algorithm is the criterion used in choosing the next bid on which to branch in step Choose. Figure 5 compares different such criteria.

The graph compares six different criteria. The auctions were taken from . The x-axis shows the number of bids and the y-axis, with a log scale, shows the average number of the nodes visited. The square root criterion of Equation 1 seems best. Nevertheless the LP coefficient criterion and the LP improved criterion described in Section 5 did very well and even seem better than the square root criterion for large numbers of bids. Figure 6 extends the previous graph to larger auctions, for the best criteria only (for the other criteria we could not finish the search). The criteria based on LP are definitely best.

## 7 Run Time

We now describe running times for auctions of different sizes, under different distributions. Figure 7 uses the distribution proposed in .

Auctions have been generated for 10, 12, 14 and 16 goods. On the x-axis: the number of bids. On the y-axis: running times in seconds (log scale). The curves are sub-linear on a logarithmic scale and exhibit the sub-exponential running time of our algorithm. This graph improves substantially on the similar graph found in , showing, again, the advantage of using LP. This graph also shows faster running times for large auctions than the ones presented in . Figure 8 presents similar results for a fixed number of bids (125), for different numbers of goods.

The curve is very clearly sub-linear. Contrary to CAMUS, our BB algorithm exhibits a sub-exponential dependence on the number of goods.

Figures 9, 10, 11 and 12 describe the running times of our algorithm over four distributions described by de Vries and Vohra . The auctions generated there are single-unit combinatorial auctions.

Figure 13 concerns auctions drawn according to the CATS multipaths distribution of the Test Suite in .

Notice that our algorithm solves quite easily even auctions with up to 20000 bids. It seems, then, that the CATS multipaths distribution generates problems that are much easier, on average, than those generated by the distribution previously proposed in . This claim is supported by 14.

For large numbers of bids, the CATS multipaths distribution indeed generates problems that are easier to solve. There is a simple explanation as to why our BB algorithm solves easily CATS auctions with very large numbers of bids: in fact, one call to the LP routine is typically enough. What happens is that the initialization phase finds an allocation of value $`v`$ and the first call to the LP routine finds a fractional solution of same value $`v`$. This means that the LP fractional problem has an integer solution and the initialization phase found an optimal solution. Our BB algorithm prunes the first node. Figure 15 shows the running times of our algorithm for the distributions of single unit combinatorial auctions suggested by CATS.

## 8 Conclusions

Our first conclusion is that, even though winner determination is theoretically intractable, it can be performed in practice for auctions of a few hundreds goods among a few thousands bids, without a super-computer and without low-level programming optimization. The suggestion of  to consider approximately-efficient mechanisms, may have been premature. The challenge is now to tackle auctions of a few thousands goods (such as the FCC auction). But for such situations, the conceptual framework we have worked with in this paper is most probably not adequate: bidders cannot be expected to describe their preferences by an explicit list of bids. Suitable languages for expressing preferences must be found and algorithmic problems studied in such a new framework.

Our second conclusion is that the difficulty in solving combinatorial auctions stems more from the number of items for sale than from the number of bids. If the number of bids is large, with respect to the number of items, Linear Programming often finds the optimal integer solution in one call.

Our third conclusion is that the method of choosing the next bid on which to branch is extremely important in determining the efficiency of the search. More work is needed to understand what makes a good choice and why the square-root criterion of Equation 1 and the LP adaptative criterion described in Section 6 are so good.


format: text

text:
# Gravitational radiation from collisions at the speed of light: a massless particle falling into a Schwarzschild black hole

## Abstract

We compute spectra, waveforms, angular distribution and total gravitational energy of the gravitiational radiation emitted during the radial infall of a massless particle into a Schwarzschild black hole. Our fully relativistic approach shows that (i) less than 50% of the total energy radiated to infinity is carried by quadrupole waves, (ii) the spectra is flat, and (iii) the zero frequency limit agrees extremely well with a prediction by Smarr. This process may be looked at as the limiting case of collisions between massive particles traveling at nearly the speed of light, by identifying the energy $`E`$ of the massless particle with $`m_0\gamma `$, $`m_0`$ being the mass of the test particle and $`\gamma `$ the Lorentz boost parameter. We comment on the implications for the two black hole collision at nearly the speed of light process, where we obtain a 13.3% wave generation efficiency, and compare our results with previous results by D’Eath and Payne.

preprint: DF/IST-02.2002

The study of gravitational wave emission by astrophysical objects has been for the last decades one of the most fascinating topics in General Relativity. This enthusiasm is of course partly due to the possibility of detecting gravitational waves by projects such as GEO600 geo , LIGO ligo  or VIRGO virgo , already operating. Since gravity couples very weakly to matter, one needs to have powerful sources of in order to hope for the detection of the gravitational waves. Of the candidate sources, black holes stand out naturally, as they provide huge warehouses of energy, a fraction of which may be converted into gravitational waves, by processes such as collisions between two black holes. As it often happens, the most interesting processes are the most difficult to handle, and events such as black hole-black hole collisions are no exception. An efficient description of such events requires the use of the full non linear Einstein’s equations, which only begin to be manageable by numerical methods, and state-of-the-art computing. In recent years we have witnessed serious progress in this field anninos , and we are now able to evolve numerically the collision of two black holes, provided their initial separation is not much larger than a few Schwarzschild radius. At the same time these numerical results have been supplemented with results from first and second order perturbation theory gleiser , which simultaneously served as guidance into the numerical codes. The agreement between the two methods is not only reassuring, but it is also in fact impressive that a linearization of Einstein’s equations yield such good results (as Smarr smarr1  puts it, “the agreement is so remarkable that something deep must be at work”). In connection with this kind of events, the use of perturbation methods goes back as far as 1970, when Zerilli zerilli1  and Davis et al davis  first computed the gravitational energy radiated away during the infall from rest at infinity of a small test particle of mass $`m_0`$, into a Schwarzschild black hole with mass $`M`$. Later, Ruffini ruffini  generalized these results to allow for an initial velocity of the test particle (this problem has recently been the subject of further study lousto , in order to investigate the question of choosing appropriate initial data for black hole collisions). Soon after, one began to realize that the limit $`m_0M`$ describing the collision of two black holes did predict reasonable results, still within perturbation theory, thereby making perturbation theory an inexpensive tool to study important phenomena.

In this paper we shall extend the results of Davis et al davis  and Ruffini ruffini  by considering a massless test particle falling in from infinity through a radial geodesic. This process describes the collision of an infalling test particle in the limit that the initial velocity goes to the speed of light, thereby extending the range of Ruffini’s results into larger Lorentz boost parameters $`\gamma `$s. If then one continues to rely on the agreement between perturbation theory and the fully numerical outputs, these results presumably describe the collision of two black holes near the speed of light, these events have been extensively studied through matching techniques by D’Eath (a good review can be found in his book D'Eath ), and have also been studied in smarr2 . The extension is straightforward, the mathematics involved are quite standard, but the process has never been studied. Again supposing that these results hold for the head on collision of two black holes travelling towards each other at the speed of light, we have a very simple and tractable problem which can serve as a guide and supplement the results obtained by Smarr and by D’Eath. Another strong motivation for this work comes from the possibility of black hole formation in TeV-scale gravity dimopoulos . Previous estimates on how this process develops, in particular the final mass of the black hole formed by the collision of relativistic particles have relied heavily upon the computations of D’Eath and Payne payne . A fresher look at the problem is therefore recommended, and a comparation between our results with results obtained years ago D'Eath ; smarr2  and with recent results eardley ; solodukhin  are in order.

Our fully relativistic results show an impressive agreement with results by Smarr smarr2  for collisions of massive particles near the speed of light, namely a flat spectrum, with a zero frequency limit (ZFL) $`(\frac{dE}{d\omega })_{\omega =0}=0.4244m_0^2\gamma ^2`$. We also show that Smarr underestimated the total energy radiated to infinity, which we estimate to be $`\mathrm{\Delta }E=0.26m_0^2\gamma ^2/M`$, with $`M`$ the black hole mass. The quadrupole part of the perturbation carries less than 50% of this energy. When applied to the head on collision of two black holes moving at the speed of light, we obtain an efficiency for gravitational wave generation of 13%, quite close to D’Eath and Payne’s result of 16% D'Eath ; payne .

Since the mathematical formalism for this problem has been thoroughly exploited over the years, we will just outline the procedure. Treating the massless particle as a perturbation, we write the metric functions for this spacetime, black hole + infalling particle, as

$$g_{ab}(x^\nu )=g_{ab}^{(0)}(x^\nu )+h_{ab}(x^\nu ),$$
(1)
where the metric $`g_{ab}^{(0)}(x^\nu )`$ is the background metric, (given by some known solution of Einstein’s equations), which we now specialize to the Schwarzschild metric

$$ds^2=f(r)dt^2+\frac{dr^2}{f(r)}+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),$$
(2)
where $`f(r)=12M/r`$. Also, $`h_{ab}(x^\nu )`$ is a small perturbation, induced by the massless test particle, which is described by the stress energy tensor

$$T^{\mu \nu }=\frac{p_0}{(g)^{1/2}}𝑑\lambda \delta ^4(xz(\lambda ))\dot{z}^\mu \dot{z}^\nu .$$
(3)
Here, $`z^\nu `$ is the trajectory of the particle along the world-line, parametrized by an affine parameter $`\lambda `$ (the proper time in the case of a massive particle), and $`p_0`$ is the momentum of the particle. To proceed, we decompose Einstein’s equations $`G_{ab}=8\pi T_{ab}`$ in tensorial spherical harmonics and specialize to the Regge-Wheeler regge  gauge. For our case, in which the particle falls straight in, only even parity perturbations survive. Finally, following Zerilli’s zerilli2  prescription, we arrive at a wavefunction (a function of the time $`t`$ and radial $`r`$ coordinates only) whose evolution can be followed by the wave equation

$$\frac{^2\stackrel{~}{𝐙}(\omega ,r)}{r_{}^2}+\left[\omega ^2V(r)\right]\stackrel{~}{𝐙}(\omega ,r)=(12M/r)S,$$
(4)
Here, the $`l`$dependent potential $`V`$ is given by

$$V(r)=\frac{f(r)\left[2\sigma ^2(\sigma +1)r^3+6\sigma ^2r^2M+18\sigma rM^2+18M^3\right]}{r^3(3M+\sigma r)^2},$$
(5)
where $`\sigma =\frac{(l1)(l+2)}{2}`$ and the tortoise coordinate $`r_{}`$ is defined as $`\frac{r}{r_{}}=f(r)`$. The passage from the time variable $`t`$ to the frequency $`\omega `$ has been achieved through a Fourier transform, $`\stackrel{~}{𝐙}(\omega ,r)=\frac{1}{(2\pi )^{1/2}}_{\mathrm{}}^{\mathrm{}}e^{i\omega t}Z(t,r)𝑑t`$. The source $`S`$ depends entirely on the stress energy tensor of the particle and on whether or not it is massive. The difference between massive and massless particles lies on the geodesics they follow. The radial geodesics for massive particles are:

$$\frac{dT}{dr}=\frac{E}{f(r)(E^21+2M/r)^{1/2}};\frac{dt}{d\tau }=\frac{E}{f(r)},$$
(6)
where $`E`$ is a conserved energy parameter: For example, if the particle has velocity $`v_{\mathrm{}}`$ at infinity then $`E=\frac{1}{(1v_{\mathrm{}}^2)^{1/2}}\gamma `$. On the other hand, the radial geodesics for massless particles are described by

$$\frac{dT}{dr}=\frac{1}{f(r)};\frac{dt}{d\tau }=\frac{ϵ_0}{f(r)},$$
(7)
where again $`ϵ_0`$ is a conserved energy parameter, which in relativistic units is simply $`p_0`$. We shall however keep $`ϵ_0`$ for future use, to see more directly the connection between massless particles and massive ones traveling close to the speed of light. One can see that, on putting $`p_0m_0`$, $`ϵ_0\gamma `$ and $`\gamma \mathrm{}`$, the radial null geodesics reduce to radial timelike geodesics, so that all the results we shall obtain in this paper can be carried over to the case of ultrarelativistic (massive) test particles falling into a Schwarzschild black hole. For massless particles the source term $`S`$ is

$$S=\frac{4ip_0e^{i\omega r_{}}ϵ_0(4l+2)^{1/2}\sigma }{w(3M+\sigma r)^2}.$$
(8)
To get the energy spectra, we use

$$\frac{dE}{d\omega }=\frac{1}{32\pi }\frac{(l+2)!}{(l2)!}\omega ^2|\stackrel{~}{𝐙}(\omega ,r)|^2,$$
(9)
and to reconstructe the wavefunction $`Z(t,r)`$ one uses the inverse Fourier transform

$$Z(t,r)=\frac{1}{(2\pi )^{1/2}}_{\mathrm{}}^{\mathrm{}}e^{i\omega t}\stackrel{~}{𝐙}(\omega ,r)𝑑\omega .$$
(10)
Now we have to find $`\stackrel{~}{𝐙}(\omega ,r)`$ from the differential equation (4). This is accomplished by a Green’s function technique. Imposing the usual boundary conditions, i.e., only ingoing waves at the horizon and outgoing waves at infinity, we get that, near infinity,

$$\stackrel{~}{𝐙}=\frac{1}{W}_{r_+}^{\mathrm{}}z_LS𝑑r.$$
(11)
Here, $`z_L`$ is a homogeneous solution of (4) which asymptotically behaves as

$`z_Le^{i\omega r_{}},r_{}\mathrm{}`$ (12)
$`z_LB(\omega )e^{i\omega r_{}}+C(\omega )e^{i\omega r_{}},r_{}+\mathrm{}.`$ (13)
$`W`$ is the wronskian of the homogeneous solutions of (4). These solutios are, $`z_L`$ which has just been defined, and $`z_R`$ which behaves as $`z_Re^{i\omega r_{}},r_{}+\mathrm{}`$. From this follows that $`W=2i\omega C(\omega )`$. We find $`C(\omega )`$ by solving (4) with the right hand side set to zero, and with the starting condition $`z_L=e^{i\omega r_{}}`$ imposed at a large negative value of $`r_{}`$. For computational purposes good accuracy is hard to achieve with the form (13), so we used an asymptotic solution one order higher in $`1/(\omega r)`$:

$`z_L=B(\omega )(1+{\displaystyle \frac{i(\sigma +1)}{\omega r}})e^{i\omega r_{}}+`$
$`+C(\omega )(1{\displaystyle \frac{i(\sigma +1)}{\omega r}})e^{i\omega r_{}},r_{}+\mathrm{}.`$ (14)
In the numerical work, we chose to adopt $`r`$ as the independent variable, thereby avoiding the numerical inversion of $`r_{}(r)`$. A fourth order Runge-Kutta routine started the integration of $`z_L`$ near the horizon, at $`r=r_i=2M+2Mϵ`$, with tipically $`ϵ=10^5`$. It then integrated out to large values of $`r`$, where one matches $`z_L`$ extracted numerically with the asymptotic solution (14), in order to find $`C(\omega )`$. To find $`Z(t,r)`$ the integral in (10) is done by Simpson’s rule. For both routines Richardson extrapolation is used. The results for the wavefunction $`Z(t,r)`$ as a function of the retarded time $`utr_{}`$ are shown in Figure 1, for the three lowest radiatable multipoles, $`l=2,3`$ and $`4`$.

As expected from the work of Ruffini ruffini  and Ferrari and Ruffini ferrari , the wavefunction is not zero at very early times, reflecting the fact that the particle begins to fall with non zero velocity. At late times, the $`l=2`$ (for example) signal is dominated by quasinormal ringing with frequency $`\omega 0.35/M`$, the lowest quasinormal frequency for this spacetime chandra . The energy spectra is shown in Figure 2, for the four lowest radiatable multipoles. First, as expected from Smarr’s work smarr2 , the spectra is flat, up to a certain critical value of the frequency, after which it rapidly decreases to zero. This ($`l`$-dependent) critical frequency is well approximated, for each $`l`$-pole, by the fundamental quasinormal frequency. In Table 1, we list the zero frequency limit (ZFL) for the first ten lowest radiatable multipoles.

For high values of the angular quantum number $`l`$, a good fit to our numerical data is

$$\left(\frac{dE_l}{d\omega }\right)_{\omega =0}=\frac{2.25}{l^3}ϵ_0^2p_0^2$$
(15)
We therefore estimate the zero ZFL as

$`\left({\displaystyle \frac{dE_l}{d\omega }}\right)_{\omega =0}=\left[{\displaystyle \underset{l=2}{\overset{l=11}{}}}\left({\displaystyle \frac{dE_l}{d\omega }}\right)_{\omega =0}\right]+{\displaystyle \frac{1}{2}}{\displaystyle \frac{2.25}{12^2}}ϵ_0^2p_0^2`$
$`=0.4244ϵ_0^2p_0^2.`$ (16)
To calculate the total energy radiated to infinity, we proceed as follows: as we said, the spectra goes as $`2.25/l^3`$ as long as $`\omega <\omega _{lQN}`$, where $`\omega _{lQN}`$ is the lowest quasinormal frequency for that $`l`$-pole. For $`\omega >\omega _{lQN}`$, $`dE/d\omega 0`$ (In fact, our numerical data shows that $`dE/d\omega e^{27\alpha \omega M}`$, with $`\alpha `$ a factor of order unity, for $`\omega >\omega _{lQN}`$). Now, from the work of Ferrari and Mashhoon ferrari2  and Schutz and Will schutz , one knows that for large $`l`$, $`\omega _{lQN}\frac{l+1/2}{3^{3/2}M}`$. Therefore, for large $`l`$ the energy radiated to infinity in each multipole is

$$\mathrm{\Delta }E_l=\frac{2.25(l+1/2)}{3^{3/2}l^3}\frac{ϵ_0^2p_0^2}{M},$$
(17)
and an estimate to the total energy radiated is then

$$\mathrm{\Delta }E=\underset{l}{}\mathrm{\Delta }E_l=0.262\frac{ϵ_0^2p_0^2}{M}$$
(18)
Let us now make the bridge between these results and previous results on collisions between massive particles at nearly the speed of light ruffini ; D'Eath ; smarr2 . As we mentioned, putting $`p_0m_0`$ and $`ϵ_0\gamma `$ does the trick. So for ultrarelativistic test particles with mass $`m_0`$ falling into a Schwarzschild black hole, one should have $`(dE/d\omega )_{\omega =0}=0.4244m_0^2\gamma ^2`$ and $`\mathrm{\Delta }E=0.262\frac{m_0^2\gamma ^2}{M}`$. Smarr smarr2  obtains

$`\left({\displaystyle \frac{dE}{d\omega }}\right)_{\omega =0}^{\mathrm{Smarr}}={\displaystyle \frac{8}{6\pi }}m_0^2\gamma ^20.4244m_0^2\gamma ^2,`$
$`\mathrm{\Delta }E_{\mathrm{Smarr}}=0.2{\displaystyle \frac{m_0^2\gamma ^2}{M}}.`$ (19)
So Smarr’s result for the ZFL is in excellent agreement with ours, while his result for the total energy radiated is seen to be an underestimate. As we know from the work of Davis et al davis  for a particle falling from infinity with $`v_{\mathrm{}}=0`$ most of the radiation is carried by the $`l=2`$ mode. Not so here, in fact in our case less than 50% is carried in the quadrupole mode (we obtain $`\mathrm{\Delta }E_{l=2}=0.1\frac{ϵ_0^2p_0^2}{M}`$, $`\mathrm{\Delta }E_{l=3}=0.0447\frac{ϵ_0^2p_0^2}{M}`$ ). This is reflected in the angular distribution of the radiated energy (power per solid angle)

$$\frac{dE}{d\mathrm{\Omega }}=\mathrm{\Delta }E_l\frac{(l2)!}{(l+2)!}\left[2\frac{^2}{\theta ^2}Y_{l0}+l(l+1)Y_{l0}\right]^2,$$
(20)
which we plot in Figure 3. Compare with figure 5 of smarr2 .

If, as Smarr, we continue to assume that something deep is at work, and that these results can be carried over to the case of two equal mass black holes flying towards each other at (close to) the speed of light, we obtain a wave generation efficiency of 13 %. This is in close agreement with results by D’Eath and Payne D'Eath ; payne , who obtain a 16% efficiency (Smarr’s results cannot be trusted in this regime, as shown by Payne payne2 ). Now, D’Eath and Payne’s results were achieved by cutting an infinite series for the news function at the second term, so one has to take those results with some care. However, the agreement we find between ours and their results lead us to believe that once again perturbation theory has a much more wider realm of validity. To our knowledge, this is the first alternative to D’Eath and Payne’s method of computing the energy release in such events.

This work was partially funded by Fundação para a Ciência e Tecnologia (FCT) through project PESO/PRO/2000/4014. V.C. also acknowledges finantial support from FCT through PRAXIS XXI programme. J. P. S. L. thanks Observatório Nacional do Rio de Janeiro for hospitality.


text:
# Potential Energy Landscape Equation of State

## Abstract

Depth, number, and shape of the basins of the potential energy landscape are the key ingredients of the inherent structure thermodynamic formalism introduced by Stillinger and Weber \[F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1982)\]. Within this formalism, an equation of state based only on the volume dependence of these landscape properties is derived. Vibrational and configurational contributions to pressure are sorted out in a transparent way. Predictions are successfully compared with data from extensive molecular dynamics simulations of a simple model for the fragile liquid orthoterphenyl.

Recent years have seen a resurgence in studies devoted to modeling the thermodynamics of supercooled liquids thermo4 ; thermo3 ; angell ; thermo2 ; wales . Such studies aim to elucidate the physics of the liquid-glass transition, to develop a thermodynamic description of out of equilibrium systems and to provide keys for a deeper understanding of the dynamics of supercooled states debenedetti01 . Numerical studies, boosted by increased computational resources which now allow simulations to track the slowing down of the dynamics over more than 6 decades in time, are providing quantitative estimates for the free energy of simple model systems lennard\_jones\_pes ; heuer97 ; scala2000 ; sastry2001 . The availability of such data provides stringent tests of the theoretical predictions scala2000 ; sastry2001 ; ivannature  and helps in the understanding of basic mechanisms associated with the behavior of thermodynamic and dynamic quantities close to the glass transition.

Among the thermodynamic formalisms amenable to numerical investigation, a central role is played by the Inherent Structure (IS) formalism introduced by Stillinger and Weber sw . Properties of the potential energy landscape (PEL), such as depth, number and shape of the basins of the potential energy surface are calculated and used in the evaluation of the liquid free energy scala2000 ; sastry2001 ; ivannature ; mossa02  In the IS formalism, the system free energy is expressed as a sum of an entropic contribution, accounting for the number of the available basins, and a vibrational contribution, expressing the free energy of the system when constrained in one of the basins sw .

Important progress has been made after the discovery that — for models of fragile liquids — the number $`\mathrm{\Omega }(e_{IS})`$ of distinct basins of depth $`e_{IS}`$ in a system of $`N`$ atoms or molecules is well described by a Gaussian distribution sastry2001 ; heuer00

$$\mathrm{\Omega }(e_{IS})=e^{\alpha N}\frac{e^{(e_{IS}E_o)^2/2\sigma ^2}}{(2\pi \sigma ^2)^{1/2}}.$$
(1)
Here the amplitude $`e^{\alpha N}`$ accounts for the total number of basins. Numerical studies of models for fragile liquids have also shown that the basin free energy can be written as the depth $`e_{IS}`$ plus a vibrational contribution which, in the harmonic approximation, has the well known form

$$F_{vib}(e_{IS},T)=k_BT\underset{i=1}{\overset{M}{}}\mathrm{ln}(\beta \mathrm{}\omega _i(e_{IS})),$$
(2)
where $`\omega _i(e_{IS})`$ is the $`i`$-th normal mode frequency ($`i=1\mathrm{}M`$) and $`\beta =1/k_BT`$. The $`M`$ normal mode frequencies define the shape of the basin. If relevant, anharmonic corrections can also be accounted for ivannature ; mossa02 . The quantity $`_{i=1}^M\mathrm{ln}(\omega _i(e_{IS})/\omega _o)`$ (where $`\omega _o`$ is the frequency unit) is found to depend linearly on the basin depth sastry2001 , i.e. can be written, in terms of two parameters $`a`$ and $`b`$, as

$$\underset{i=1}{\overset{M}{}}\mathrm{ln}(\omega _i(e_{IS})/\omega _o)=a+be_{IS}.$$
(3)
Hence, the vibrational free energy can be written as

$$F_{vib}(e_{IS},T)=F_{vib}(E_o,T)+k_BTb(e_{IS}E_o).$$
(4)
Within the two assumptions of Eq. (1) — Gaussian distribution of basin depths — and Eq. (4) — linear dependence of the basin free energy on $`e_{IS}`$ — an exact evaluation of the partition function can be carried out. The corresponding Helmholtz free energy is given by sastry2001

$`F(T)`$ $`=`$ $`TS_{conf}(T)+e_{IS}(T)`$
$`+`$ $`F_{vib}(E_o,T)+k_BTb(e_{IS}(T)E_o),`$
where

$$e_{IS}(T)=(E_ob\sigma ^2)\beta \sigma ^2=e_{\mathrm{}}\beta \sigma ^2,$$
(6)
and

$`S_{conf}(T)/k_B`$ $`=`$ $`\alpha N(e_{IS}(T)E_o)^2/2\sigma ^2`$ (7)
$`=`$ $`S_{\mathrm{}}/k_Bb\beta \sigma ^2\beta ^2\sigma ^2/2.`$
In the above equations, $`e_{\mathrm{}}`$ and $`S_{\mathrm{}}`$ are defined as the value of $`e_{IS}`$ and $`S_{conf}`$ at infinite $`T`$. Eqs. (Potential Energy Landscape Equation of State)-(7) show that, along constant volume $`V`$ paths, the behavior of the thermodynamic quantities is controlled by the values of the PEL properties, as given by $`\alpha `$, $`\sigma `$, $`E_0`$ (from Eq. (1)) and by $`a`$ and $`b`$ (from Eq. (3)).

In this Letter we study the volume dependence of Eq. (Potential Energy Landscape Equation of State) to provide an expression for the equation of state (EOS) based completely on landscape properties (PEL-EOS) debenedettieos . This study provides a significant insight into the understanding of the pressure $`P`$ and opens the way for detailed comparisons between experimental measurements (usually performed at constant $`P`$) and theoretical approaches based on the IS formalism. It may also help in developing an IS-based thermodynamic description of out-of-equilibrium (glass) states and a theoretical definition of the concepts of fictive $`P`$ and $`T`$aging .

In thermodynamics, $`P`$ is defined as the (negative) $`V`$ derivative of the Helmholtz free energy. Hence $`P`$ if fully determined by the $`V`$ dependence of the landscape properties $`\alpha `$, $`\sigma `$, $`E_0`$, $`a`$ and $`b`$. Eq. (Potential Energy Landscape Equation of State) shows that $`P`$ can be split into three main contributions: a configurational one, $`P_{conf}`$ — related to the change in the number of available basins with $`V`$; an $`e_{IS}`$ one, $`P_{e_{IS}}`$ — related to the change in basin depth with $`V`$; and a vibrational one, $`P_{vib}`$ — related to the change in the shape of the explored basin with $`V`$. The $`T`$ dependence of each contribution can then be studied independently. The explicit expressions for these contributions are:

$`P_{conf}(T,V)=T{\displaystyle \frac{}{V}}S_{\mathrm{}}{\displaystyle \frac{}{V}}(b\sigma ^2){\displaystyle \frac{1}{T}}{\displaystyle \frac{}{V}}(\sigma ^2/2k_B)`$ (8)
$$P_{e_{IS}}(T,V)=\frac{}{V}e_{\mathrm{}}+\frac{1}{T}\frac{}{V}(\sigma ^2/k_B)$$
(9)
$$P_{vib}(T,V)=k_BT\frac{}{V}(a+be_{\mathrm{}})+\frac{}{V}(b\sigma ^2).$$
(10)
By grouping together all the contributions with the same $`T`$ dependence, $`P`$ can be expressed in term of $`V`$ derivatives of only three combinations of the five PEL parametersanomali

$$P(T,V)=P_{const}+TP___T+T^1P_{_{_{1/T}}},$$
(11)
where $`P_{const}=e_{\mathrm{}}/V`$, $`P___T=S_{\mathrm{}}/Vk_B\frac{}{V}(a+be_{\mathrm{}})`$ and $`P_{_{_{1/T}}}=\frac{}{V}(\sigma ^2/2k_B)`$.

The present formalism also provides an expression for the so-called inherent structure equation of state (IS-EOS), $`P_{IS}(T_{IS},V)`$ debenedettieos ; laviolette ; roberts ; sastry00 , i.e. the relation between the pressure and volume of the typical IS and the temperature $`T_{IS}`$ of the equilibrium ensemble from which configurations were extracted. Indeed, the constant $`V`$ steepest descent minimization procedure, which realizes the search for the closest local minimum starting from an equilibrium configuration, suppresses all the vibrational components (hence $`P_{vib}=0`$) but keeps $`P_{conf}`$ and $`P_{e_{IS}}`$ frozen at their “equilibrium” initial value. As a result $`P_{IS}`$, a purely mechanical quantity, can be expressed as

$`P_{IS}(T_{IS},V)=P_{conf}+P_{e_{is}}=`$
$`{\displaystyle \frac{}{V}}E_o+T_{IS}{\displaystyle \frac{}{V}}S_{\mathrm{}}+{\displaystyle \frac{1}{k_BT_{IS}}}{\displaystyle \frac{}{V}}(\sigma ^2/2).`$ (12)
The present approach also predicts the behavior of $`P`$ when an IS configuration is heated from $`T=0`$ at constant $`V`$. While the system remains in the same basin, $`P`$ is given by $`P_{IS}(T_{IS},V)+P_{vib}(T,T_{IS},V)`$ where the only (linear) $`T`$-dependent part arises from

$$P_{vib}(T,T_{IS},V)=k_BT\frac{}{V}(a+be_{IS}(T_{IS})).$$
(13)
We now apply the present derivation to the simple Lewis and Wanström (LW) model for the fragile molecular liquid orthoterphenyl (oTP) mossa02 ; lewis . The LW model is a rigid three-site model, with intermolecular site-site interactions described by the Lennard Jones (LJ) potential lewis . Its simplicity allows one to reach simulation times of the order of $`\mu s`$ mossa02 . In this model, as in the LJ case, the anharmonic contributions are negligible, $`e_{IS}(T)`$ goes as $`1/T`$, and $`_{i=1}^M\mathrm{ln}(\omega _i(e_{IS})/\omega _o)`$ is linear in $`e_{IS}`$ mossa02 . Hence this model is a perfect candidate for testing the validity of the PEL-EOS proposed here.

We use the excellent data base of state points presented in Ref. mossa02  i) to calculate the $`V`$ dependence of the PEL parameters; ii) to derive the EOS for the oTP model, and iii) to compare it with the “exact” EOS calculated using the virial expression, as commonly implemented in molecular dynamics (MD) codes.

Fig. 1 shows, for five constant $`V`$ paths, the simultaneous fit of the $`T`$ dependence of $`S_{conf}`$, the basin curvatures, and of $`e_{IS}(T)`$, according to Eqs. (7), (3), and (6). The possibility of fitting simultaneously, with the same values of $`\alpha `$, $`\sigma `$, $`E_0`$, $`a`$ and $`b`$, the quantities $`e_{IS}(T)`$, $`_{i=1}^M\mathrm{ln}(\omega _i(e_{IS})/\omega _o)`$ and $`S_{conf}(T)`$, supports the validity of the two main assumptions, i.e. Eq. (1) and Eq. (4).

The $`V`$ dependence of the three combinations of fit parameters, $`e_{\mathrm{}}`$, $`S_{\mathrm{}}k_B(a+be_{\mathrm{}})`$ and $`\sigma ^2/2k_B`$ is shown in Fig.2. $`P(V,T)`$ can be calculated from the $`V`$ derivative of these quantities according to Eq. (11) and compared with the MD data. Such comparison is reported in Fig.3. The very good agreement between the two set of data at all $`V`$’s confirms that a very accurate EOS based on PES properties has been derived for this oTP model.

The $`T`$ dependence of the individual contributions can be evaluated according to Eqs. (8)-(10). Fig. 4 shows — along a constant $`V`$ and a constant $`P`$ path — $`P_{vib}`$ and $`P_{conf}+P_{e_{IS}}`$. We note that at constant $`V`$, both $`P`$ components are increasing with $`T`$, while at constant $`P`$, $`P_{conf}+P_{e_{IS}}`$ decreases on heating to compensate for the increase of $`P_{vib}`$.

The PES-EOS allows us also to contrast the isobaric and isochoric $`T`$ dependence of $`S_{conf}`$, $`e_{IS}(T)`$ and $`_{i=1}^M\mathrm{ln}(\omega _i/\omega _o)`$. Such comparison is reported in Fig. 1. Here we note the faster decrease of $`S_{conf}(T)`$ and $`e_{IS}(T)`$ along constant $`P`$ paths as well as the different trend in the change of basin shapenota .

As a further check of the quality of the calculated EOS for the oTP model, Fig. 5 compares the MD (virial) and PEL (Eq. (12)) IS-EOS. Since for the oTP model the term linear in $`T`$ in Eq. (12) is negligible lanave02 , $`P_{IS}`$ follows, to a good extent, a $`1/T_{IS}`$ law. The quality of the comparison supports the interpretation of $`P_{IS}`$ as the $`V`$ derivative of the depth and number of basins sampled in the corresponding thermodynamic equilibrium state.

Finally, Fig. 3 compares MD data (open symbols) and PEL (dashed lines) predictions in a run where $`T`$ is increased starting from the $`T=0`$ IS configuration, as previously discussed. The entire simulation is much shorter than the time needed to change basin, so that only the vibrational degrees of freedom are thermalized. Also in this case, the PEL expression accounts for the observed linear increase of $`P`$.

The present EOS, based exclusively on PEL properties, can and will be used to address important issues in the thermodynamics of supercooled liquids footnotebivariata , such as the study of the intrinsic limit of stability of the liquid state debenedettieos ; roberts ; sastry00  and the IS-based thermodynamic description of aging aging .


text:
# On the Mass-Period Correlation of the Extrasolar Planets

## 1 INTRODUCTION

The mass distribution of the extrasolar planets was recognized to be a key feature of the growing new population since the first few detections (e.g., Basri & Marcy 1997; Mayor, Queloz & Udry 1998; Mazeh, Goldberg & Latham 1998; Heacox 1999; Mazeh 1999; Stepinski & Black 2000). Recent studies showed that the mass distribution is probably flat or slightly decreasing in log M (Jorissen, Mayor & Udry 2001; Zucker & Mazeh 2001; Tabachnik & Tremaine 2002; Lineweaver & Grether 2002), and has a distinct cutoff at about 10 Jupiter masses(=$`M_J`$). This paper focuses on a possible dependence between the extrasolar planets masses and their orbital periods.

Every study of the mass-period relation has to take into account the strong observational selection effect that prohibits the detection of low-mass–long-period planets, because of their small radial-velocity amplitudes. Tabachnik & Tremaine (2002) studied the mass and the period distributions simultaneously, assuming they were two independent power-law distributions. They found that the uncertainties of the exponents of the two variables are highly correlated, but attributed their findings to the observational selection effect. This paper shows that the mass-period correlation found in the sample of known extrasolar planets cannot be attributed solely to the observational selection effect. We show that there is an additional real dependency between the mass and the period of the extrasolar planets, manifested as a significant paucity of high-mass–short-period planets. Since such planets are the easiest to detect, this paucity is probably not the result of any selection effect.

In Section 2 we present our analysis of the mass-period correlation of the whole sample of known planets. Section 3 shows that the small subsample of planets in binary stellar systems may have different, opposite, correlation. Section 4 shortly discusses our findings in terms of some theories for the existence of giant planets close to their host stars, the migration model in particular. A preliminary version of this work was presented in Mazeh & Zucker (2002).

## 2 Analysis

Figure 1 presents the minimum masses of all known extrasolar planets as a function of their orbital periods. The data were taken from the web-site of the California Planet Search Team111 http://exoplanets.org/planet\_table.txt, and were updated as of December 2001. We chose to plot the two axes with logarithmic scales, because the frequency of planets, up to 10 $`M_J`$, is almost flat in log M (Jorissen, Mayor & Udry 2001; Zucker & Mazeh 2001; Tabachnik & Tremaine 2002), as well as in log P (Heacox 1999; Stepinski & Black 2001; Mazeh & Zucker 2002; but see a somewhat different approach by Tabachnik & Tremaine 2002).

We concentrated in this analysis on a trapezoidal area in the minimum-mass–period parameter space bounded by dashed lines in the figure. The upper boundary corresponds to the 10 $`M_J`$ cutoff in planets masses (e.g., Jorissen, Mayor & Udry 2001; Zucker & Mazeh 2001). Although Zucker & Mazeh (2001) suggested a probable small higher-mass tail beyond the 10 $`M_J`$ line, the distribution is flat up to, probably, 10 $`M_J`$. We plot the five objects above the line for completeness. The two vertical lines represent the minimum and maximum orbital periods found in the sample.

The ascending line at the bottom of the figure corresponds to a constant radial-velocity amplitude, $`K`$, of 25 m s-1. The detection rate of the present planet-search projects below this line is low. This is easily seen in the figure, which includes only five planets below this border line. Again, we plotted these five planets only for the sake of completeness. We assume that planets detected above that line have all been reported. We now proceed to analyze the 66 planets inside the trapezoid, assuming a constant detection rate over its area.

A close examination of Figure 1 reveals a paucity of planets at the high-mass–short-period corner of the trapezoid. Only three planets appear at that corner. This is certainly not a selection effect, because planets at that part of the diagram have the largest radial-velocity amplitudes, and therefore are the easiest to detect.

It is not clear yet what is the shape of the area in which we find low frequency of planets. It might have, for example, a rectangular shape bordered by $`\mathrm{log}P=1.6`$ and $`\mathrm{log}(M_2\mathrm{sin}i)=0.3`$, or could have a wedge shape, bordered by the line from $`(\mathrm{log}P,\mathrm{log}(M_2\mathrm{sin}i))=(0.46,0.2)`$ to $`(1.5,1)`$. In any case, it seems that there are enough planets in the trapezoid to render this paucity significant.

To estimate quantitatively the statistical significance of the high-mass–short-period paucity seen in the figure we first consider the mass-period correlation coefficient of the sample of planets in the trapezoid. The resulting value was 0.661. We claim that this high value means there is a real correlation in the planets population — higher than the one induced by the selection effect. In terms of statistical hypotheses testing, we have to reject the null hypothesis that there is no mass-period correlation in the planets population, and the correlation we find in the sample comes solely from the wedge-shape of the area removed by the selection effect.

To assess the statistical significance of the null hypothesis rejection we used Monte-Carlo simulations in which we created an artificial sample, randomly drawn out of a two-dimensional uniform distribution in log-mass and log-period, between the period limits of the trapezoid, and between 0.175 and 10 $`M_J`$. To simulate the selection effect we have discarded every planet whose implied radial velocity was too small and drawn another one instead, until we had in hand 66 planets. The process was repeated $`10^6`$ times, calculating the correlation for each simulated sample. Figure 2 shows the histogram of the simulated correlation coefficients together with the value derived for the original sample. The simulated correlation values distributed around a mean value of 0.336 with a standard deviation of 0.104. Only 208 simulations yielded a value larger than 0.661. Thus, we can conclude that the null hypothesis of two uniform uncorrelated random distributions of log-mass and log-period can be rejected with a 99.98% confidence level. Of course, in the currently small sample of planets, even a single additional detection of a planet with extremely short period and large mass can alter the results significantly.

The seminal work of Tabachnik & Tremaine (2002) has used maximum likelihood calculation to study the mass and period distributions of the extrasolar planets. Assuming the two distributions have power-law shape and are mutually independent, they derived a positive power of $`0.26\pm 0.06`$ for the period distribution and a negative one, $`0.12\pm 0.10`$, for the mass distribution. Such a distribution creates a deficiency of high-mass–short-period planets that might have produced the effect we report here. In order to check whether this is the case, we repeated the simulation process with Tabachnik and Tremaine’s distribution. Indeed, the correlation values were somewhat higher than in the flat distribution case, with an average value of 0.403 and standard deviation 0.109. However, here again only 3823 simulations out of $`10^6`$ yielded a value larger than 0.661. We can, therefore, reject this hypothesis with a 99.62% confidence level. We also checked the extremely unlikely case where the true values of both exponents are 2-$`\sigma `$ away from the ones derived by Tabachnik & Tremaine (2002). Two independent distributions, with $`0.38`$ for the exponent of the period distribution and $`0.32`$ for that of the mass distribution, are still rejected at a 98% confidence level. The somewhat lower significance of the last rejection indicates that we can not reject all possible hypotheses where the mass and the period are uncorrelated, and one might come up with a specific distribution that will reproduce the reported effect, without a real correlation between the two variables.

## 3 Planets in Binary Stars

Having established the significance of the paucity of the short-period massive planets, we can examine Figure 1 and try to see what distinguishes the few planets in the high-mass–short-period corner. The three planets that appear to be somewhat isolated in that corner of the figure are $`\tau `$ Boo b, HD 195019 b and Gl 86 b — all of them are planets found in wide stellar binaries (Hale 1994; Fischer et al. 1999; Els et al. 2001). This raises the possibility that planets in binary stellar systems have a different mass-period distribution. If this is true, it can tell us about the possible effect the binarity of the parent star might have on the formation (Boss 1998; Nelson 2000) and orbital evolution of their planets (e.g., Mazeh, Krymolowski & Rosenfeld 1997; Holman, Touma & Tremaine 1997; Innanen et al. 1997; Holman & Wiegert 1999).

In order to check this hypothesis we assembled a subsample of all the planet-hosting stars that we can safely tag as binaries. Using the WDS catalog, we decided, somewhat arbitrarily, to consider only those binaries whose projected separation is smaller than 1000 AU, assuming wider binaries would not have influenced the formation and evolution of their planets. We also discarded binaries that have angular separations larger than 10″, whose secondaries were observed to be fainter than 12th magnitude, and had no other evidence for a physical association (mainly common proper motion). We added to the final list Gl 86 which was announced after the most recent publication of the WDS. We were left with nine binary stars: HD 142, Gl 86, HD 19994, $`ϵ`$ Eri, $`\tau `$ Boo, HD 178911 B, 16 Cyg, HD 195019, HD 217107. Figure 3 shows the two separate subsamples of planets. $`ϵ`$ Eri lies outside our nominal trapezoid and thus only eight planets constitute our binary subsample.

The difference between the two subsamples in Figure 3 is striking. Apart from the enhanced paucity of planets in the upper left corner of the figure in the single-star population, the binary population shows an opposite trend. A negative-slope straight line, with a slope of $``$0.15, can be fitted to this small subsample. This is opposed to the slope of the single-star planets, for which we fit a straight line with a slope of $`+`$0.33. As we showed in the previous section, only part of this positive slope is due to the specific shape of the trapezoid considered, and the population of planets does show a positive slope because of the short-period–high-mass paucity. The remarked difference is also reflected by the negative mass-period correlation of $``$0.459 for the binary planets.

The difference between the two populations is manifested also in the higher value of the correlation coefficient of the subsample of 58 single-star planets — 0.783. In order to test the significance of the difference between the sample with and without the binary planets, we ran Monte-Carlo simulations again. Each iteration consisted of removing a random set of 8 planets from the original sample of 66 and re-calculating the correlation. The results of $`10^6`$ iterations are depicted in Figure 4. The simulated values average at 0.661, the original value of the parent sample, and have a standard deviation of 0.028. Only 17 out of $`10^6`$ simulations yielded a value higher than 0.783, thus implying a significance of 99.998% to the higher correlation of the single-star planets.

## 4 Discussion

We have presented evidence for a substantial deficit of massive planets with short orbital periods. We have shown that the correlation seen in the sample of known planets is higher than the one predicted by published independent power-law models and the selection effect, although we need more points to corroborate our finding. We regard this deficit as a refinement of the initial surprising discoveries of giant planets in close orbits (e.g., Mayor & Queloz 1995). We have shown that only planets with masses below about 2 $`M_J`$ can frequently be found orbiting single stars with periods shorter than about 40 days. This may serve to refine or update the models that were devised to explain the existence of 51 Peg-like planets.

We can divide the models for the existence of planets with short orbital periods into three broad categories. The most favored model is that of planetary migration (e.g., Lin, Bodenheimer, & Richardson 1996). In this model a planet is formed by core accretion at a distance of the order of 5 AU or further (but see Bodenheimer, Hubickyj & Lissauer 2000 for a somewhat different approach), and then is pushed toward the parent star by interaction with the accretion disk. Alternative models are migration by interaction with other planets (e.g., Weidenschilling & Marzari 1996) or planetesimals (Murray et al. 1998). A completely different approach assumes planet formation by disk instability (e.g., Boss 1997a). Obviously, if the instability ends up as a planet far away from the parent star, one needs a migration mechanism to account for the close-in planets (Boss 1997b). In what follows we will try to comment on the implications of our findings on each of the three categories of models.

Two effects within the migration scenario can contribute to the paucity of massive planets with small orbits:

(I) Massive planets open a gap in the disc, and consequently slow their migration rate substantially (e.g., Ward 1997; Trilling et al. 1998; Nelson et al. 2000). We expect, therefore, to find the more massive planets at distances closer to their formation sites.

(II) Trilling et al. (1998) pointed out that when planets get too close to their parent stars they loose substantial fraction of their mass through Roche-lobe overflow. In the specific parameters presented by Trilling et al., planets above 3.4 $`M_J`$ do not migrate significantly from their formation site. Most of the planets with initial masses below 3.4 $`M_J`$ loose substantial fraction of their mass through Roche-lobe overflow. In fact, planets with initial masses below 3.36 $`M_J`$ loose most of their mass, and are left with masses smaller than 0.4 $`M_{}`$.

The two effects, discussed already by theoretical studies of the migration model, contribute to the paucity of very massive planets with small orbits. The second effect contributes to the mass-radius correlation for short distances, on the order of tenths of an AU, whereas the first effect dominates the correlation for larger distances. Actually, Trilling et al. (1998) already published a mass-period diagram with a series of models they ran (see their Figure 7), in which the effect we find here can be clearly seen.

However, the first effect can be canceled out if the planet mass depends on the disk mass. This is so because the size of the gap depends on the ratio between the planet mass and the disk mass (e.g., Trilling et al. 1998). Suppose, for example, that the planet mass scales with the disk mass. Then more massive planets could move inwards before they opened a gap as much as less massive planets do. Therefore, to account for the correlation we see in the data we have to assume that planets are formed with masses that do not scale with the disk mass. In other words, planets with different masses can be formed in disks with similar masses.

Interaction with other planet(s) was suggested mainly to explain the high eccentricities observed for some of the known extrasolar planets (e.g., Weidenschilling & Marzari 1996, Rasio & Ford 1996; Ford, Havlickova & Rasio 2001). A few models include an accompanying disk that absorbs the angular momentum necessary to enable the migration (e.g., Murray, Paskowitz & Holman 2001).

Suppose the interaction with the other, as yet undetected, planet is the dominant mechanism for the migration of the known planet. Such a scenario can account for the observed mass-period correlation if the mass of the undetected planet is independent of the mass of the known planet. In such a case, the ability of the unseen planet to push the known planet to smaller radii is limited only to small-mass planets, consistent with the mass-period correlation we see in the data. If, on the other hand, the mass ratio between the planets is similar in all planetary systems, the migration caused by the planet-planet interaction should not be limited to small-mass planets, contrary to our findings.

The same argument applies to migration driven by interaction with planetesimals. Murray et al. (1998) comment that “if the mass of the planetesimal disk interior to the planet is of order of the planet mass, the planet can migrate nearly to the surface of the star”. In other words, the migration depends on the mass ratio between the planet and the planetesimal disk. Therefore, to account for the observed mass-period correlation we have to assume that the planet mass does not scale with the mass of the planetesimal disk.

Planet formation via disk instability (Boss 1997a; 1998a; 2000; 2001) could also account for the mass-period correlation, in principle, if we assume an in-situ formation for this model. We can speculate that the mass of the forming planet depends on the mass available in the disk at the vicinity of the instability. Suppose, for example, that the mass associated with the instability is some fraction of the mass within the radius of the instability center. This can lead to more massive planets forming at large distances, hence the effect we have detected. However, one still needs to establish by detailed numerical simulations that this model can work at small distances from the star. If, on the other hand, the disk instability can work only at large distances, and thus still requires a mechanism for migration, then dependence of the planet mass on the disk mass might make it difficult to account for the paucity we detected.

Two theoretical studies have considered the implication of binarity of a star on the formation of planets around one of its components. One report on numerical study (Boss 1998b) indicated that the presence of a stellar companion can induce a rapid instability even for disks that are stable otherwise. A somewhat more recent work (Nelson 2000) suggested an opposite effect, claiming that planets are unlikely to form in certain binary systems. Both works dealt only with a companion at a distance of 40–50 AU, and more comprehensive studies are needed to explore the implication of a companion on planet formation. In any case, both works indicate that the planets in binary systems might not have the same mass-period distribution, consistent with our findings.

To summarize, it seems as if almost all models for the existence of close-in giant planets can account for the mass-period correlation of the single-star planets, although the correlation seems a more natural outcome of the model that assumes migration driven by a disk-planet interaction. In any case, this correlation can put, in principle, some constraints on the different models. The number of planets known today is only marginally enough to characterize the details of the short-period–massive-planets deficit, apart from establishing its existence. Nevertheless, we already can conclude that for all models we need planet masses that do not scale with the mass of the disk/planetesimals/other planets. More data can illuminate the finer details of this phenomenon and help to better tune the theories for close-in giant planets.

We wish to thank S. Tremaine for important comments on an earlier version of the manuscript. We are indebted to the referee, M. Holman, for his careful reading of the paper, his useful comments and advise, and for pointing out an error in one of our simulations. This research has made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory, and was supported by the Israeli Science Foundation (grant no. 40/00).


text:
# 1 Introduction

## 1 Introduction

Our understanding of Gamma-Ray Bursts (GRBs) was revolutionized during the last ten years. According to the generally accepted Fireball model (see e.g. ) the gamma-rays and their subsequent multiwavelength afterglow are produced when an ultra-relativistic flow is slowed down. All current observations, from prompt emission to late afterglow, from $`\gamma `$-rays via X-ray optical and IR to radio, are consistent with this model.

Observations of GRB host galaxies revealed the association of (long) GRBs (there is no information on the positions of short GRBs) with star forming Galaxies  and that these (long) GRBs follow the star formation rate . There is evidence (still inconclusive) that some long bursts are associated with Supernovae . These observations suggest that the progenitors of GRBs are massive stars. There are some reasonable ideas how does a collapsing star  or a merging binary  produce the required $`10^{51}`$ergs. But it is not clear how does the GRB’s “inner engine” accelerates and collimated the relativistic flow. Today this is the most interesting (and most difficult) open question concerning GRBs. This is also the most relevant question for this conference where the possibility that transitions to strange stars power GRBs has been discussed .

We summarize, here, the known constrains on the “inner engines” of GRBs. These constrains arise mostly from the temporal structure of GRBs. Several years ago Sari and Piran  have shown that variable GRBs can be generated only by internal interactions111This interaction is usually considered as a collisionless shock. However the exact nature of the interaction is unimportant for most of our arguments. within the flow. To produce internal shocks the central engine must produce a long and variable wind. This leads to a powerful NO GO theorem: Variable GRBs cannot be produced from a single explosion. This NO GO theorem rules out explosive GRB models that produce the relativistic flow in a single explosion. Kobayashi et al  have shown that the observed internal shocks light curve reflects almost directly the temporal activity of the inner engine. This is the best direct evidence on what is happening at the center of the GRB.

We review the arguments leading to these conclusions. We also discuss new observational results  and a new theoretical toy model  that explains these observations within the internal shocks paradigm. This toy model suggests that the “inner engine” is producing a variable Lorentz factor wind by modulating the mass ejection of a roughly constant energy flow. The other alternative of modulating the energy of a constant mass flow is ruled out. We conclude by summarizing the various constrains on the “inner engines”. We leave to the reader the task of examining the implications of these results to his/hers favorite model222See  for a discussion of the implications for possible accretion models.

## 2 Energetics and Beaming

The most important factor in any model is the total amount of energy that it releases. Redshift measurements have lead to alarming estimates of more than $`10^{54}`$ergs in some bursts . When factoring in the efficiency the requirements exceed a solar rest mass energy.

However, these early estimates assumed isotropic emission. Jet breaks in afterglow light curves lead to estimates of the beaming factors. When those are taken into account we find a “modest” practically constant energy release of $`10^{51}`$ergs . This lower energy budget allows for many possible models. At the same time it introduces an additional requirement on the central engine. It has to collimate the relativistic flow to narrow beaming angles (at times or order of $`1^o2^o`$).

We cannot estimate directly the total energy released by the “inner engine”. However, here are two possible estimates: $`E_\gamma `$, the energy released as $`\gamma `$-rays and $`E_K`$, the kinetic energy during the adiabatic afterglow phase. Remarkably both energies are comparable. This last observation implies that the conversion efficiency of the initial relativistic kinetic energy to $`\gamma `$-rays must be very high.

## 3 Time Scales In GRBs - Observations

Most GRBs are highly variable. Fig 1 depicts the light curve of a typical variable GRB (GRB920627). The variability time scale, $`\delta t`$, is determined by the width of the peaks. $`\delta t`$ is much shorter (in some cases by a more than a factor of 100) then $`T`$, the duration of the burst. Variability on a time scale of milliseconds is seen in some long bursts . However, not all bursts are variable. We stress that our discussion applies to variable bursts and it is not applicable to the small subset of smooth ones.

A comparison of the pulse width distribution and the pulse separation, $`\mathrm{\Delta }t`$ , distribution, reveals an excess of long intervals . These long intervals can be classified as quiescent periods , relatively long periods of several dozen seconds with no activity. When excluding the quiescent periods we  find that both distributions are lognormal with a comparable parameters: The average pulse interval, $`\overline{\mathrm{\Delta }}t=1.3sec`$ is larger by a factor 1.3 then the average pulse width $`\overline{\delta }t=1sec`$. One also finds that the pulse widths are correlated with the preceding interval .

The results described so far are for long bursts. the variability of short ($`T<2`$sec) bursts is more difficult to analyze. The duration of these bursts is closer to the limiting resolution of the detectors. Still we find that most ($`66\%`$) short bursts are variable with $`\delta t/T<0.1`$ .

## 4 Time Scales In GRBs - Theory

Consider a spherical relativistic emitting shell with a radius $`R`$, a width $`\mathrm{\Delta }`$ and a Lorentz factor $`\mathrm{\Gamma }`$. This can be a whole spherical shell or a spherical like section of a jet whose opening angle $`\theta `$ is larger than $`\mathrm{\Gamma }^1`$. Because of relativistic beaming an observer would observe radiation only from a region of angular size $`\mathrm{\Gamma }^1`$. Photons emitted by matter moving directly towards the observer (point A in Fig. 3) will arrive first. Photons emitted by matter moving at an angle $`\mathrm{\Gamma }^1`$ (point D in Fig. 3) would arrive after $`t_{ang}=R/2c\mathrm{\Gamma }^2`$. This is also the time of arrival of photons emitted by matter moving directly towards the observer but emitted at $`2R`$ (point C in Fig. 3). Thus, $`t_{rad}t_{ang}`$ . This coincidence is the first part of the NO GO theorem that rules out single explosions a sources of GRBs.

At a given point particles are continuously accelerated and emit radiation as long as the shell with a width $`\mathrm{\Delta }`$ is crossing this point. The photons emitted at the front of this shell will reach the observer a time $`T_\mathrm{\Delta }=\mathrm{\Delta }/c`$ before those emitted from the rear (point B in Fig. 3). In fact photons are emitted slightly longer as it takes some time for the accelerated photons to cool. For most reasonable parameters the cooling time is much shorter from the other time scales  and we ignore it hereafter.

Light curves are divided to two classes according to the ratio between $`T_\mathrm{\Delta }`$ and $`t_{ang}t_{rad}`$. The emission from different angular points smoothes the signal on a time scale $`t_{ang}`$. If $`T_\mathrm{\Delta }t_{ang}t_{rad}`$ the resulting burst will be smooth with a width $`t_{ang}t_{rad}`$. The second part of the NO GO theorem follows from the hydrodynamics of external shocks. Sari and Piran  have shown that for external shocks $`\mathrm{\Delta }/cR/c\mathrm{\Gamma }^2t_{rad}t_{ang}`$. External shocks can produce only smooth bursts!

A necessary condition for the production of a variable light curve is that $`T_\mathrm{\Delta }=\mathrm{\Delta }/c>t_{ang}`$. This can be easily satisfied within internal shocks (see Fig 4). Consider an “inner engine” emitting a relativistic wind active over a time $`T_\mathrm{\Delta }=\mathrm{\Delta }/c`$ ($`\mathrm{\Delta }`$ is the overall width of the flow in the observer frame). The source is variable on a scale $`L/c`$. The internal shocks will take place at $`R_sL\mathrm{\Gamma }^2`$. At this place the angular time and the radial time satisfy: $`t_{ang}t_{rad}L/c`$. Internal shocks continue as long as the source is active, thus the overall observed duration $`T=T_\mathrm{\Delta }`$ reflects the time that the “inner engine” is active. Note that now $`t_{ang}L/c<T_\mathrm{\Delta }`$ is trivially satisfied. The observed variability time scale in the light curve, $`\delta t`$, reflects the variability of the source $`L/c`$. While the overall duration of the burst reflects the overall duration of the activity of the “inner engine”.

Numerical simulations  have shown that not only the time scales are preserved but the source’s temporal behaviour is reproduced on an almost one to one basis in the observed light curve. We will return to this point in section 6 in which we describe a simple toy model that explains this result.

## 5 Caveats and Complications

Clearly the way to get around the NO GO theorem is if $`t_{ang}<t_{rad}`$. In this case one can identify $`t_{rad}`$ with the duration of the burst and $`t_{ang}`$ as the variability time scale. The observed variability would require in this case that: $`t_{ang}/t_{rad}=\delta t/T`$.

One can imagine an inhomogeneous external medium which is clumpy on a scale $`dR/\mathrm{\Gamma }`$ (see Fig 5). Consider such a clump located at an angle $`\theta \mathrm{\Gamma }^1`$ to the direction of motion of the matter towards the observer. The resulting angular time, which is the difference in arrival time between the first and the last photons emitted from this region would be:$`d/c\mathrm{\Gamma }`$. Now $`t_{ang}<t_{rad}`$ and it seems that one can get around the NO GO theorem.

Sari and Piran  have shown that such a configuration would be extremely inefficient. This third component of the NO GO theorem rules out this caveat. The observations limit the size of the clumps to $`dc\mathrm{\Gamma }\delta t`$ and the location of the shock to $`RcT\mathrm{\Gamma }^2`$. The number of clumps within the observed angular cone with an opening angle $`\mathrm{\Gamma }^1`$ equals the number of pulses which is of the order $`T/\delta t`$. The covering factor of the clumps can be directly estimated in terms of the observed parameters by multiplying the number of clumps $`T/\delta t`$ times their area $`d^2=(\delta t\mathrm{\Gamma })^2`$ and dividing by the cross section of the cone $`(R/\mathrm{\Gamma })^2`$. The resulting covering factor equals $`\delta t/T1`$. The efficiency of conversion of kinetic energy to $`\gamma `$-rays in this scenario is smaller than this covering factor which for a typical variable burst could be smaller than $`10^2`$.

We turn now to several attempts to find a way around this result. We will not discuss here the feasibility of the suggested models (namely is it likely that the surrounding matter will be clumpy on the needed length scale , or can an inner engine eject “bullets”  or “cannon balls”  with an angular width of $`10^2`$ degrees and what keeps these bullets so small even when they are shocked and heated). We examine only the question whether the observed temporal structure can arise within these models.

### 5.1 External Shocks on a Clumpy Medium

Dermer and Mitman  claim that the simple efficiency argument of Sari and Piran  was flawed. They point out that if the direction of motion of a specific blob is almost exactly towards the observer the angular time would be of order $`d^2/cR`$. This is narrower by a factor $`d\mathrm{\Gamma }/R`$ than the angular time across the same blob that is located at a typical angle of $`\mathrm{\Gamma }^1`$. These special blobs would produce strong narrow peaks and will form a small region along a narrow cone with a larger covering factor. Dermer and Mitman  present a numerical simulation of light curves produced by external shocks on a clumpy inhomogenous medium with $`\delta t/T10^2`$ and efficiencies of up to $`10`$%.

A detailed analysis of the light curve poses, however, several problems for this model. While this result is marginal for bursts with $`\delta t/T10^2`$ it is insufficient for bursts with $`\delta t/T10^3`$. Variability on a time scale of milliseconds is observed  in many long GRBs (namely $`\delta t/T`$ can be as samll as $`10^4`$.). Moreover, in this case we would expect that earlier pulses (that arise from blobs along the direction of motion) would be narrower than latter pulses. This is in a direct contradiction with the observations . Finally there is no reason to expect the observed similarity between the pulse width and the pulse separation and the correlation between the pulse width and the preceding interval in this model.

### 5.2 The Shot-Gun and the Cannon Ball Models

Heinz and Begelman  suggested that the “inner engine” operates as a shot-gun emiting multiple narrow bullets with an angular size much smaller than $`\mathrm{\Gamma }^1`$ (see Fig 6). These bullets do not spread while propagating and they are slowed down rapidly by an external shock with a very dense circumburst matter. The pulses width is $`t_{ang}`$ or the slowing down time while the duration of the burst is determined by the time that the “inner engine” emits the bullets. While the cannon ball model of Dar and De Rujulla  is drastically different from a physical point of view it is rather similar in terms of its temporal features. Hence the following remarks apply to this model as well.

This model satisfies our NO GO theorem in the sense that also here the burst is not produced by a single explosion. Moreover the observed light curve represents also here the temporal activity of the source. However, it is based on external shocks.

The most serious problem is the fact that the width of the pulses here is determined by the angular time or the hydrodynamic time (which in turn depends on the external density profile of the circumburst matter) while the intervals between the pulses depend on the activity of the inner engine. There is no reason why the two distributions will be similar and why there should be a correlation between them.

## 6 An Internal Shocks Toy Model

The discovery  that the distribution of pulse widths and pulse separations are comparable and that there is a correlation between the pulse width and the preceding interval provides an independent evidence in favor of the internal shocks model. Furthermore it suggests that the different shells emitted by the internal engine are most likely “equal energy” rather than “equal mass” shells.

The similarity between the pulse width and the pulse separation distribution and the correlation is extremely unlikely to arise from a single shell passing thought a random distribution of clumps that surround the “central engine”. In this case the arrival time of individual pulses will depend on the position of the emitting clumps relative to the observers. Two following pulses could arise from two different clumps that are rather distant from each other. There is no reason why the pulses and intervals should be correlated in any way. A similar problem arises in the “shot-gun” or the “cannon-ball” models in which the pulse duration is determined by one parameter and the separation by another.

Both features arise naturally within the internal shocks model  in which both the pulse duration and the separation between the pulses are determined by the same parameter. We outline here the main arguments showing that. Consider two shells with a separation $`L`$. The slower outer shell Lorentz factor is $`\mathrm{\Gamma }_1=\mathrm{\Gamma }`$ and the inner faster shell Lorentz factor is $`\mathrm{\Gamma }_2=a\mathrm{\Gamma }`$ ($`a>2`$ for an efficient collision), both in the observer frame. The shells’ are ejected at $`t_1`$ and $`t_2t_1+L/c`$. The collision takes place at a radius $`R_s2\mathrm{\Gamma }^2L`$ (Note that $`R_s`$ does not depend on $`\mathrm{\Gamma }_2`$). The emitted photons from the collision will reach the observer at time (omitting the photons flight time, and assuming transparent shells):

$$t_ot_1+R_s/(2c\mathrm{\Gamma }^2)t_1+Lt_2.$$
(1)
The photons from this pulse are observed almost simultaneously with a (hypothetical) photon that was emitted from the “inner engine” together with the second shell (at $`t_2`$). This explains why Kobayashi et al  find numerically that for internal shocks the observed light curve replicates the temporal activity of the source.

Consider now four shells emitted at times $`t_i`$ ($`i=1,2,3,4`$) with a separation of the order of $`L`$ between them. Assume that there are two collisions - between the first and the second shells and between the third and the fourth shells. The first collision will be observed at $`t_2`$ while the second one will be observed at $`t_4`$. Therefore, $`\mathrm{\Delta }tt_4t_22L/c`$. Now assume a different collision scenario, the second and the first shells collide, and afterward the third shell takes over and collide with them (the forth shell does not play any roll in this case). The first collision will be observed at $`t_2`$ while the second one will be observed at $`t_3`$. Therefore, $`\mathrm{\Delta }tt_3t_2L/c.`$ Numerical simulations  show that more then 80% of the efficient collisions follows one of the two scenarios described above. Therefore we can conclude:

$$\mathrm{\Delta }tL/c.$$
(2)
Note that this result is independent of the shells’ masses.

The pulse width is determined by the angular time (ignoring the cooling time): $`\delta t=R_s/(2c\mathrm{\Gamma }_s^2)`$ where $`\mathrm{\Gamma }_s`$ is the Lorentz factor of the shocked emitting region. If the shells have an equal mass ($`m_1=m_2`$) then $`\mathrm{\Gamma }_s=\sqrt{a}\mathrm{\Gamma }`$ while if they have equal energy ($`m_1=am_2`$) then $`\mathrm{\Gamma }_s\mathrm{\Gamma }`$. Therefore:

$$\delta t\{\begin{array}{cc}\hfill R_s/2a\mathrm{\Gamma }^2cL/ac& \mathrm{equal}\mathrm{mass},\hfill \\ \hfill R_s/2\mathrm{\Gamma }^2cL/c& \mathrm{equal}\mathrm{energy}.\hfill \end{array}$$
(3)
The ratio of the Lorentz factors $`a`$, determines the collision’s efficiency. For efficient collision the variations in the shells Lorentz factor (and therefore $`a`$) must be large.

It follows from Eqs. 2 and 3 that for equal energy shells the $`\mathrm{\Delta }t`$-$`\delta t`$ similarity and correlation arises naturally from the reflection of the shells initial separation in both variables. However, for equal mass shells $`\delta t`$ is shorter by a factor of $`a`$ than $`\mathrm{\Delta }t`$. Since $`a`$ has a large variance this would wipes off the $`\mathrm{\Delta }t`$-$`\delta t`$ correlation. This suggests that equal energy shells are more likely to produce the observed light curves.

## 7 Conclusions

We cannot provide a recipe for a GRB “inner engine”. However we can list the specifications of this engine (for a long variable GRB). It must satisfy the following conditions:

* It should accelerate $`10^{51}`$ergs to a variable relativistic flow with $`\mathrm{\Gamma }>100`$.
* It should collimate this flow, with a varying degree of collimation (up to $`1^o`$).
* It should be active from several seconds up to several hundred seconds (according to the duration of the observed burst).
* It should vary on a time scale of seconds or less (corresponding to the duration of a typical pulse within the burst).
* Different shells of matter should have a comparable energy and their different Lorentz facors should arise due to a modulation of the accelerated mass.
* At times the engine should stop for several dozen seconds (resulting in a quiescent periods).

Before concluding we stress that these specification are for long variable burst (which compose the majority of the observed bursts). Many of these (but not all) apply also to short variable bursts (about two thirds of the short bursts). These specifications do not apply to smooth bursts (either short or long ones).

This research was supported by a grant from the US-Israel Binational Science Foundation.


text:
# Fabry-Perot observations at the 6m telescope

## 1.  SCORPIO design

The 2D spectroscopy observations with IFP at the 6m telescope using the CIGALE system were begun in early 80s by teams from Marseille Observatory (France) and from SAO (Boulesteix et al., 1982; Amram et al., 1992). In 1997 a CCD was attached on the old focal reducer instead of photon counter. A new multi-mode focal reducer SCORPIO (Spectral Camera with Optical Reducer for Photometrical and Interferometrical Observations) was developed in the SAO RAS in 2000 (Afanasiev et al., this Proceeding). Here we describe only Fabry-Perot mode of the SCORPIO. For more details see our WWW-site
`http://www.sao.ru/~moisav/scorpio/scorpio.html`.

The optics of SCORPIO compensates the aberrations of the main telescope mirror and provide a total focal ration $`F/2.9`$. The Queensgate IFPs ET-50 installs into the 35mm collimated beam. SCORPIO are under full remote control from a PC. A total DQE (telescope$`+`$filter$`+`$IFP$`+`$CCD) is about $`20\%`$ in the $`H_\alpha `$ spectral region. SCORPIO equipped with narrow filters ($`FWHM=1015`$ Å) for IFP observations in the $`H_\alpha `$ line for object’s velocities $`200\mathrm{}+10.000\text{km}\text{s}^1`$.

## 2.  Calibration and data reduction

Focal reducer mounted on the prime focus universal adapter which has two movable probes for offset guiding and continuum and spectral lamps for calibration. A He-Ne-Ar lamp and continuum spectrum lamp mounted in the integration sphere which produced a uniform lighting. The optics of calibration system forms the $`F/4`$ beam on the reducer entrance. There is no systematic shift in velocities because telecentrical beams used for wavelength calibration.

A CCD has a large advantage in quantum effeciency in comparison with a photon counter. However a readout noise and a reading time do not allow to do very short exposures of the object channels. Therefore exist a problem of the photometry correction of flux and seeing in the individual channels because only 1-2 scanning cycles per channel may be accepted on CCD observations, we should take in account a time dependence of a sky brightness and of an atmospheric extinction. Therefore wavelength object cube contains a set of narrow rings (Fig. 1c). This pattern limits the velocities measurements for the low brightness objects. We suggest a method to resolve this problem. Main idea is removing the sky counts before cube linerization. After standard CCD-reduction a mean radial profile of the sky lines is constructed by the azimuthal averaging of the interferometrical rings (Fig. 1a). The sky radial profile are subtracted from the object. Fig. 1b shows example of the sky subtracted images. Variations of an atmospheric extinction, seeing and channel offset shifts could be measured from the photometry of background stars and corrected in the cube. Unfortunately our method requires that about half of field of view will be free from the object emission, but it allows to obtain a good results in the object cube.

## References

Amram P., Marcelin M., Bonnarel F., Boulesteix J., Afanas’ev V.L., Dodonov S.N., 1992, A&A, 263, 69

Boulesteix J., Georgelin Y., Marcelin M., Fort J.A., 1982, in ”Instrumentation for astronomy with large optical telescopes”, 223