text:
# Jets in GRBs: Tests and Predictions for the Structured Jet Model

## 1 Introduction

The degree to which gamma-ray bursts (GRBs) and their afterglows are beamed is an important issue. A proper understanding of the geometry of the relativistic outflow affects the total energetics of GRB central engines and the GRB event rates, both of which are crucial parameters for constraining possible progenitor models. Evidence for non-isotropic outflows is believed to come from observations of achromatic breaks in afterglow light curves (e.g., Rhoads 1999, Sari et al. 1999) and the detection of polarized emission (Covino et al. 1999, Wijers et al. 1999). We should however note that, whereas this is a natural interpretation for the breaks, in some cases other interpretations have been proposed or different conclusions derived (e.g. Nicastro et al. 1999; Vrba et al. 2000). Here we adopt the point of view that the observed breaks are indeed manifestations of jets.

In early theoretical papers (e.g., Rhoads 1997, Sari et al. 1999) it was assumed for simplicity that the ejecta had to be distributed approximately uniform across the entire opening angle in the gamma-ray phase and that the majority of the explosive energy in the afterglow phase must have a single bulk Lorentz factor. In this “uniform” model, a break occurring in the light curves at $`t_{\mathrm{break}}`$ can be directly translated into a jet with an opening angle $`\theta _{\mathrm{break}}`$. Using this simple framework, Frail et al. (2001) carried out an analysis of all known afterglows and found that there was a distribution of jet opening angles leading to a reduction in the gamma-ray energy from its isotropic value with relatively small scatter. The observed distribution was shown to be heavily weighted towards small opening angles.

Rather than positing a uniform jet, it is equally reasonable to assume that GRB jets are structured in some fashion. In the collapsar progenitor model (e.g. Wang, Woosley & MacFadyen 2002) the Lorentz factor and energy are high near the rotation axis, but decrease off axis as the degree of entrainment increases. Salmonson (2000) has argued for such a jet structure to explain the empirical correlation between the GRB peak luminosity and pulse lag (Norris et al. 2000). In this case, the distribution of observed break times in afterglow light curves are not due to a distribution of opening angles but rather originate from variations in the viewing angle of a structured jet (Postnov, Prokhorov & Lipunov 2001). In two recent papers by Rossi, Lazzati & Rees (2002a) and Zhang & Meszaros (2002) it was shown that a jet with a universal beaming configuration could reproduce the near-constant energy result of Frail et al. (2001) provided that the energy per unit solid angle (and Lorentz factor) varied as the inverse square of angular distance from the jet symmetry axis.

Discriminating between the uniform and structured jet models is important as they yield different estimates for the true GRB event rate and the total energy, besides leading to clues on the physical mechanism producing the jet itself. Afterglow lighcurves have argued to be degenerate to the structure of the jet (Rossi et al. 2002a, Zhang & Meszaros 2002) and therefore they cannot be used as diagnostics. Rossi et al. (2002b) showed that detection of polarization can provide useful constraints. Here we concentrate on diagnostics based on geometrical effects. In particular, in the structured jet model, since the inferred opening angle is just a geometric effect of the viewing angle, it is possible to predict the distribution of angles and compare it to the observed distribution of $`\theta _{\mathrm{break}}`$ by Frail et al. (2001) and Bloom et al. (2003). Unfortunately, a specific prediction is not possible for the uniform model, since there is no framework for jet formation which yields the distribution of opening angles.

In this paper, we work out the expected distribution of burst opening angles, under the structured jet scenario and given various assumption of the star formation rate evolution in the universe. We find two competing effects. First, even though randomly oriented bursts would be rarely observed on axis, they are much brighter and therefore can be seen to larger distances rendering small opening angles common. Cosmological effects limit the volume at large redshifts leading to an effective cutoff at the small opening angle. As a result, we predict that the most common opening angle should be about 0.12 rad. Furthermore, under the structured jet model, we predict that more sensitive future missions like SWIFT, will find a much larger typical angle, about 0.3 rad.

The observed data so far is still a small sample. Selection effects are hard to quantify, especially for large and small opening angles. In addition, the variety of instruments with different sensitivity used to detect the bursts, makes a robust conclusion difficult at this stage. Yet, we show that the current set of data, as given in the most updated sample of Bloom et al. (2003), is at least in rough agreement with the distribution we predict for the structured jet model. If this conclusion would hold with much larger and less biased sample, it would be a strong support to the structured jet model.

## 2 Computation of the observed distribution of jet angles

### 2.1 Scalings for a Euclidian universe

If all bursts were observable, then we would expect that the number $`dn(\theta )/d\theta `$ of bursts with angle in the interval $`d\theta `$ around $`\theta `$ would be proportional to $`\theta `$. This implies that most of the observed bursts should have a large angle, which is in complete contradiction with observations. Zhang & Meszaros (2002) attributed this apparent discrepancy to the small sample size or afterglow selection effects. However, this argument does not account for the fact that bursts with small $`\theta `$ are brighter by a factor of $`\theta ^2`$, and therefore can be seen (in an Euclidian universe) up to a distance $`\theta ^1`$ farther, which contains a volume larger by a factor of $`\theta ^3`$. The expected distribution in Euclidean geometry is therefore expected to be $`dn/d\theta \theta ^2`$. Although this is closer to what observations suggest, we will show that this now exaggerates the number of small-angle bursts compared to a proper cosmological prediction: GRBs originate at redshifts of order unity, and therefore suffer considerable cosmological corrections. This is why their log $`N`$ \- log $`S`$ curve does not obey the Euclidian $`S^{3/2}`$, but is shallower at low $`S`$. It is for the same reason that the number of observed bursts of low $`\theta `$ will not be as high as predicted by the Euclidian $`\theta ^2`$. In the following, we work out these cosmological effects in detail.

### 2.2 Cosmological effects

Let $`R_{\mathrm{GRB}}(z)`$ be the GRB rate per unit comoving volume per unit time, then the total (i.e. over the all sky) rate of bursts with inferred jet angle between $`\theta `$ and $`\theta +d\theta `$ is given by

$$\frac{dn(\theta )}{d\theta }=\mathrm{sin}\theta _0^{z_{\mathrm{max}}(\theta )}𝑑z\frac{R_{\mathrm{GRB}}(z)}{(1+z)}\frac{dV(z)}{dz},$$
(1)
where $`z_{\mathrm{max}}(\theta )`$, is the maximal redshift up to which we can observe a burst with apparent angle $`\theta `$. This redshift is found by numerically inverting the equation

$$F_{ph,\mathrm{lim}}=\frac{L_{ph}(\theta )}{4\pi D^2(z_{\mathrm{max}})(1+z_{\mathrm{max}})^\alpha },$$
(2)
where $`F_{ph,\mathrm{lim}}`$ is the limiting photons flux (photons per unit area per unit time) that is detectable by the GRB detector with frequency range, $`\nu _l<\nu <\nu _u`$. $`L_{ph}(\theta )`$ is the photons luminosity (photons per unit time) in the same frequency rage but in the local frame of the burst, of a burst with an apparent angle $`\theta `$. A factor of $`(1+z)^{\alpha 1}`$ is a spectral correction, assuming that the GRB has a differential photon spectral index $`\alpha `$, and another factor of $`1+z`$ takes care of time dilation. For BATSE, $`\nu _l=50`$ keV and $`\nu _u=300`$ keV.

The normalization constant is determined by the condition

$$\frac{L_{\nu _1\nu _2}}{4\pi }\mathrm{\hspace{0.17em}2}\pi \theta ^2T=E$$
(3)
where $`E10^{51}\mathrm{ergs}`$ is the roughly constant energy of GRBs as inferred by (Bloom et al. 2003). That relation used the luminosity integrated over the frequency range $`\nu _1=20`$ keV to $`\nu _2=2000`$ keV as calculated by Bloom et al. (2001). For powerlaw photon spectrum, this is related to the photon luminosity in the triggering band $`\nu _l`$ $`\nu _h`$ by

$$L_{ph}(\theta )=\frac{2E}{\theta ^2Th\nu _l}\frac{\alpha 2}{\alpha 1}\left(\frac{\nu _l}{\nu _1}\right)^{\alpha +2}\frac{1(\nu _u/\nu _l)^{\alpha +1}}{1(\nu _2/\nu _1)^{\alpha +2}}$$
(4)
If we take a spectrum with $`\alpha 1`$, (which is appropriate for the frequency range $`50300`$ keV as reported by Mallozzi et al. 1996) then we obtain,

$$L_P(\theta )=1.1\times 10^{57}T^1\theta ^2\mathrm{ph}\mathrm{sec}^1.$$
(5)
It should be noted that $`T`$ is not the total duration of the burst, but an “effective” duration of convenience here, that is the duration that the burst would have if its energy output were constant at the peak value rather than highly variable. In the simplest version of our model (§2 & 3), we will assume a single value of $`T`$ for all the bursts; however we will explore (§4) how our results vary when a scatter in $`T`$ is introduced (which is a more realistic assumption).

The jet model with the energy profile $`\theta ^2`$ also makes detailed predictions for the observed GRB flux distribution. This, within our formalism, can be written as

$$\frac{dn(S)}{dS}=_0^{\mathrm{}}𝑑L_pf(L_p)\left[\frac{R_{\mathrm{GRB}}(z)}{(1+z)}\frac{dV(z)}{dz}\left|\frac{dz}{dS}(z,L_p)\right|\right]_{z=z(S,L_p)}$$
(6)
where, given Eq. (5), the luminosity function takes the form $`f(L_p)L_p^2`$.

In Eqs.(1) and (6), $`dV(z)/dz`$ is the comoving volume. In a flat cosmology with a cosmological constant it is given by

$$\frac{dV(z)}{dz}=4\pi D^2(z)\frac{dD(z)}{dz}$$
(7)
where $`D(z)`$ is the comoving distance,

$$D(z)=\frac{c}{H_0}_0^z\frac{dz^{}}{\sqrt{(1+\mathrm{\Omega }_mz)(1+z^{})^2\mathrm{\Omega }_\mathrm{\Lambda }(2z^{}+z^2)}},$$
(8)
We assume a cosmological model with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`H_0=71`$ km s-1 Mpc-1.

We assume that GRBs trace the star formation history and we adopt, as our “standard” model for $`z10`$, the Rowan-Robinson star formation rate which can be fitted with the expression

$$R_{\mathrm{GRB}}(z)=\{\begin{array}{cc}R_0\mathrm{\hspace{0.33em}10}^{0.75z},\hfill & z<z_{\mathrm{peak}}\hfill \\ R_0\mathrm{\hspace{0.33em}10}^{0.75z_{\mathrm{peak}}},\hfill & \text{ }zz_{\mathrm{peak}}\hfill \end{array},$$
(9)
where $`z_{\mathrm{peak}}2`$. For $`z10`$, in our standard model we use an interpolation that follows the star formation history derived in numerical simulations by Gnedin & Ostriker (1997). While using this SFR as our working model, we also explore the effects on the predicted distribution $`dn/d\theta `$ of different star formation histories, and in particular we consider two opposite extremes, one in which the SFR does not rapidly decline for $`z10`$ as implied in the numerical simulations of Gnedin & Ostriker, and another, the Madau curve (Madau 1996), in which the SFR rapidly declines at redshifts $`z3`$.

For each model for the SFR, the normalization constant $`R_0`$ is determined by the condition

$$_0^{\pi /2}𝑑\theta \frac{dn(\theta )}{d\theta }=R_{\mathrm{GRB}}^{\mathrm{obs}},$$
(10)
where $`R_{\mathrm{GRB}}^{\mathrm{obs}}=667`$ yr-1 is the observed BATSE rate, and $`F_{\mathrm{lim}}`$ in Eq. (2) is the BATSE threshold flux for which this rate has been measured. We adopted the 90% efficiency peak flux threshold for BATSE, that is $`F_{\mathrm{lim}}=0.424`$ ph/sec (e.g. Mallozzi, Pendleton & Paciesas 1996).

## 3 Comparison with data and predictions for more sensitive surveys

In order to compare the theoretical distributions derived in §2.2, we require a sample of gamma-ray bursts whose values of $`\theta _{\mathrm{break}}`$ have been measured. The largest published sample of $`\theta _{\mathrm{break}}`$ values at the time of our work comes from the analysis of 28 bursts with redshifts and well-studied afterglow light curves by Bloom et al. (2003). This list of 28 includes all bursts with measured redshifts at this time. For all but four, some limit on $`\theta _{\mathrm{break}}`$ was derived. The exceptions are GRB 970228, in which sparse data together with likely contribution from a supernovae makes it difficult to interpret the lightcurve, GRB 990506 where no sufficient data exists, GRBs 980425 which had no optical counterpart, and GRB 021211 which is being analyzed at the time of writing. For eight of the remaining 24, upper or lower limits where put on the opening angle. We have not used these limits in our comparison. One could think that this would tend to have the effect of narrowing our sample distribution. However, the upper and lower limits do not tend to be at the edge of the distribution of the measured opening angles (see Table 2 of Bloom et al. 2003). They therefore do not necessarily reflect extreme cases, but cases with lower quality data. However, it should be noted that in three cases (GRBs 970828, 991216, 990705) the confidence that indeed a jet has been identified is weak, since the break was observed in a single frequency only. Finally, we should remark that the inferred values of $`\theta _{\mathrm{break}}`$ have some uncertainties. These values are indeed computed using the expression given in Frail et al. (2001): $`\theta _{\mathrm{break}}t_{\mathrm{break}}^{3/8}(1+z)^{3/8}E_{\mathrm{iso}}^{1/8}\eta _\gamma ^{1/8}n^{1/8}`$. The measured values of $`t_{\mathrm{break}}`$ have relative errors within 30%, while for the densities a value of $`n=10`$ cm-3 is assumed for the 5 (out of 16) cases for which the data quality did not allow a self-consistent determination of the density through broad-band afterglow modelling. Finally, Bloom et al. assumed a constant value ($`\eta _\gamma =0.2`$) for the efficiency of the bursts. Whereas the dependence on $`\eta _\gamma `$ is weak, a large spread in this not-well constrained quantity would introduce a further source of error in the determination of $`\theta _{\mathrm{break}}`$. The combination of these caveats prevents a solid comparison of the data with our predictions, and the following comparison should be taken as a general guide for this type of analysis, while showing that the observed distribution so far does not seem to be in contradiction with the prediction of the structured jet hypothesis.

We performed a Kolmogorov-Smirnov (KS) test to assess the compatibility of the theoretical distributions with the data, and we found (see Figure 1), for the Rowan-Robinson SFR, a probability of $`90\%`$ that the data are drawn from the theoretical distribution. In order to have a reasonable agreement with the data, a value $`T8`$ sec is needed in the theoretical model.

As shown in Figure 1, there is little difference in the results (and in the required value of $`T`$ to fit the data) between the case where the SFR drops rapidly after $`z10`$ (model 1), and where it keeps constant also at higher redshifts (model 2). The similarity between the distributions in model 1 and model 2 is a consequence of the combined effect of the decrease in volume at those high redshifts and the increased time dilution of the observed rate. On the other hand, for the same value of $`T`$ (or equivalently peak luminosity), the distribution that uses the Madau SFR (model 3) predicts significantly more events with large angles. This is because this star formation rate drops abruptly at a redshift $`z34`$, and therefore there is no much gain in the number of small-$`\theta `$ (i.e. brighter) events which can be seen at higher redshifts. A KS test showed that the model 3 distribution is consistent with the data (at the 40% level) if a value $`T25`$ is used. A lower value for the peak luminosity is needed to shift the $`dn/d\theta `$ distribution to lower $`\theta `$’s.

For a jet model to be self-consistent it is necessary that, if the distribution $`dn/d\theta `$ has a good agreement with the data for certain model parameters, the corresponding $`dn/dS`$ has to have a good agreement with the corresponding data for the same set of parameters. We found that with $`T8`$ (as required in model 1) the peak fluxes are within a factor of a few for those bursts with measured $`\theta `$. More generally, when comparing the theoretical $`dn/dS`$ distribution with the all BATSE catalogue we find a very good agreement for peak fluxes in the range $`1P15`$ ph/cm2/s, and a departure (as an overprediction) at higher fluxes. However, it should be noted that, given the large number of bursts in the BATSE catalogue, a comparison with the $`dn/dS`$ data distribution is much more sensitive to the model parameters than the $`dn/d\theta `$ comparison. In a recent analysis, Lloyd et al. (2001) found that a good fit can be obtained with a luminosity function $`L^{2.2}`$ (which would reduce the number of high $`P`$ bursts with respect to our model) and a redshift evolution. Such a detailed analysis is beyond the scope of this paper given the lack of comparable wealth of data for the $`\theta `$ distribution 111It should also be noted that the value of the peak luminosity that best fits the $`dn/d\theta `$ distribution depends on the adopted value of the detection efficiency, which may vary from burst to burst.. On the other hand, we should remark that, if we adopt the value $`T25`$ required to fit the $`dn/d\theta `$ distribution with the Madau SFR, then the luminosities that we predict are generally smaller than those measured for the bursts with known redshift. Therefore, in the following we will only use our model 1 with $`T=8`$ for further calculations. The corresponding cumulative distribution is shown with the solid line in the two panels in Figure 2 where it is compared to the binned, cumulative data. We should note that the size of the bins in the figure has been chosen so that there is an equal increment in the distribution for every new data point. This gives a better visual idea of the data distribution, but does not reflect the actual magnitude of the errors as described above.

Figure 3 shows the predicted evolution of the distribution of observed break angles with increasing sensitivity of the survey. In particular, we considered the 100% efficiency sensitivities corresponding to BATSE (solid line), HETE-2 and Swift (dashed line), and an intermediate sensitivity (dotted line). An interesting prediction of the luminosity distribution in Eq.(5) is that the average observed jet angle is an increasing function of the survey sensitivity. There are two counteracting effects that determine the average observed jet angle $`<\theta >`$ of the sample as a function of the survey sensitivity. As $`F_{\mathrm{lim}}`$ decreases, $`z_{\mathrm{max}}`$ increases, bringing into the sample a fraction of bursts with larger redshifts and correspondingly smaller $`\theta `$ (which are the most luminous). On the other hand, the higher sensitivity also brings into the sample a fraction of bursts with larger $`\theta `$ at lower redshifts, and this latter effect dominates over the former, partly as a result of the volume-reduced and time-dilated rate of the high-$`z`$ bursts.

## 4 Extensions

All the above results have been produced under the assumption of a strict correlation between the total energy of the burst and its peak flux. However, as discussed in §2, this relation has a scatter; therefore we have also investigated the extent to which our results change when a dispersion in the distribution of the values of $`T`$ (i.e. in the relation 5) is introduced. If we parameterize the scatter in $`T`$ with a probability distribution, $`P(T)`$, then the distribution of jet angles (1) is generalized to

$$\frac{dn(\theta )}{d\theta }=2\pi \mathrm{sin}\theta _0^{\mathrm{}}𝑑TP(T)_0^{z_{\mathrm{max}}(\theta ,T)}𝑑z\frac{R_{\mathrm{GRB}}(z)}{(1+z)}\frac{dV(z)}{dz}.$$
(11)
We took the probability distribution for the scatter, $`P(T)`$, to be a log-gaussian distribution with mean equal to $`T=8`$, and studied the dependence of the break-angle distribution $`dn/d\theta `$ on the width of the distribution $`\sigma _T`$. A K-S test shows that, with a scatter $`\sigma _T=0.3`$, the theoretical distribution is still compatible with the data at the 80%, while with $`\sigma _T=0.5`$ the agreement is at the 40% level.

All the results so far have been derived under the assumption of an energy distribution from the jet axis $`\theta ^2`$ in the interval $`0\theta \pi /2`$. However, close to the axis this divergence must naturally have a cutoff which we represent by a core of size $`\theta _c`$. We now explore how our results vary by allowing for the presence of a core in the inner part of the jet and an outer cutoff at some large angle $`\theta _j`$, where the luminosity drops rapidly to zero rather than following the profile in Eq. (5). In this case the peak photon luminosity is given by $`L_P(\theta )`$ as in Eq. (5) for $`\theta _c<\theta <\theta _j`$, and by $`L_P(\theta =\theta _c)`$ for all the angles $`\theta \theta _c`$. Hence the cumulative $`\theta `$ distribution is given by

$$N(<\theta )=\{\begin{array}{cc}0\hfill & \theta <\theta _c\hfill \\ _0^{\theta _c}𝑑\theta ^{}\frac{dn[\theta ^{},L(\theta _c)]}{d\theta ^{}}N_c\hfill & \text{ }\theta =\theta _c\hfill \\ N_c+_{\theta _c}^\theta 𝑑\theta ^{}\frac{dn[\theta ^{},L(\theta ^{})]}{d\theta ^{}}N_j\hfill & \text{ }\theta _c<\theta <\theta _j\hfill \\ N_j\hfill & \theta >\theta _j\hfill \end{array}.$$
(12)
For the case where the increase in luminosity saturates at an angle $`\theta _c`$ from the jet axis (left panel of Fig. 2), the number of observed bursts with break angle $`\theta <\theta _c`$ is smaller than the corresponding number at the same observed $`\theta `$ for the distribution with $`\theta _c=0`$ (solid line in the figure). This is because $`z_{\mathrm{max}}`$ saturates to the value given by the solution of Eqn. (2) with $`L=L_c`$ for all $`\theta \theta _c`$, whereas when there is no core222Note that when we say here “no core” or $`\theta _c=0`$, we mean an infinitesimally small core, as there would be a formal divergence in the energy if $`\theta _c`$ were precisely equal to zero., $`z_{\mathrm{max}}`$ is larger for the jet angles $`\theta <\theta _c`$. The situation is reversed in the case of a jet with $`\theta _c=0`$ but a total aperture $`\theta _j<\pi /2`$ (right panel of Fig.2). The number of bursts with observed break angle $`\theta >\theta _j`$ is zero. Therefore, in order for the normalization (i.e. total number of observed bursts) to be the same for any $`\theta _j`$, the cumulative number $`N(<\theta )`$ for the case $`\theta _j<\pi /2`$ must be larger than the corresponding number for the case $`\theta _j=\pi /2`$ (solid line in the figure).

The probability remains of the same order by varying the core angle $`\theta _c`$ in the range $`0\theta _c0.055`$, and drops as $`\theta _c`$ is increased above 0.055 (which is the smallest angle in the data set). Similarly, no significant variation in the probability is found as the outer boundary of the jet angle, $`\theta _j`$ is decreased from $`\pi /2`$ to 0.55 rad, which is the largest observed break angle. In short, the current data is too poor to constrain either the size of the core, $`\theta _c`$, or the outer size of the jet, $`\theta _j`$, beyond the trivial statement that this range must include the range of observed opening angles. It should however be remarked that, whereas the current data on the observed $`\theta _{\mathrm{break}}`$ do not allow to pin down the values of the model parameters $`\theta _c`$ and $`\theta _j`$, the type of analysis that we are proposing has the potentiality to further constrain details of the model once a larger sample of jet opening angles is gathered.

## 5 Conclusions

A natural framework for the interpretation of achromatic breaks in the afterglow light curves is the presence of jets in the GRB ejecta. There is however a certain degree of degeneracy in the resulting light curves between the case where the jet has uniform energy profile and a sharp edge, and where it has instead a variable energy profile. Distinguishing between the two scenarios is especially important in order to have a proper estimate of the GRB event rates, and a better understanding of the physics of the GRB explosion. A particularly useful discriminant of the two jet scenarios is the distribution of the observed break times, which can be theoretically predicted in the structured jet model, and compared to the data.

In this paper, we have derived such distribution and compared it to the observed sample of data on $`\theta _{\mathrm{break}}`$. We found that the observed data set, altough a small sample, is consistent with the predictions of the structured-jet model. We should however remark that the alternative “uniform” model, in which there is a direct correlation between the observed opening angle and the physical jet angle, cannot be ruled out by this approach. However, if, as more data becomes available, the agreement with the predictions of the structured jet model remains intact, it would be very contrived to justify it within the framework of the uniform jet. This model, in fact, does not make any prediction for the distribution of opening angles. One would then need to invoke an explanation for why there is the same number of bursts for each logarithmic interval of openening angle.

Besides performing a first attempt to test the structured-jet model, we have shown that this model offers a number of predictions which may be testable with future, larger datasets. The predicted opening angle distribution shown in Figure 1 (which assumes the BATSE detection threshold), has a distinct peak at $`\theta 0.12`$ rad. We predict that future, more sensitive missions will predict more bursts with large opening angles. Specifically, we estimate that the opening angle distribution for the SWIFT mission will peak near $`\theta 0.3`$ rad. Surprisingly, in this model the average redshift is only a weak function of the sensitivity and consequently we expect no increase in the average redshift detected by the next generation of gamma-ray instruments.

Another prediction that the structured-jet model makes, regards the number of GRBs in the local universe. Our “standard” SFR model yields $`R_{\mathrm{GRB}}(z=0)0.5`$ Gpc-3 yr-1. Combining this with the local galactic density $`0.0048`$ Mpc-3 (Loveday et al. 1992) one obtains $`0.1`$ GEM (galactic events per Myr). This is an interesting number for detection of local GRB remnants (Loeb & Perna 1998; Efremov et al. 1998; Perna et al. 2000). For the uniform jet model, precise rates are more difficult to estimate, as they depend on the assumed luminosity function, or equivalently the intrinsic jet angle distribution which is not known apriori for this model. If one were to assume that $`dn/d\theta \theta `$ as in the structured jet model, this would result in a logarithmic correction factor, $`(1+\mathrm{log}[(\pi /2)/\theta _c])`$, to the rate that we estimated here for the structured jet model333The correction is expected to be less than an order of magnitude as $`\theta _c1/\mathrm{\Gamma }`$.. The total energy output of GRBs would however remain the same, as in the structured jet model the total energy of each burst is also corrected by a logarithmic factor with respect to the total energy of each burst in the uniform case.

###### Acknowledgements.

We thank the referee for a very careful and thoughtful review of our manuscript. RP and DAF thank the California Institute of Technology for its kind hospitality during the time that part of this work was carried out. RS holds a senior Sherman Fairchild fellowship. This research was partially supported by a NASA grant to RS. The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.


format: text

text:
# A Planetary Companion to 𝛾 Cephei A

## 1 Introduction

The first high-precision radial velocity survey for planetary companions to nearby stars was conducted with the use of a HF gas absorption cell on the Canada-France-Hawaii Telescope (Campbell & Walker 1979; Walker et al. 1995). The $`\gamma `$ Cephei system (HR 8974 = HD 222404 = HIP 116727) was one of 16 stars on the observing list of this program. Campbell, Walker, & Yang (1988) reported that two of the stars on the observing list ($`\gamma `$ Cephei and $`\chi ^1`$ Orionis) were previously unknown single-lined spectroscopic binaries. In the case of $`\gamma `$ Cephei, Campbell, Walker, & Yang (1988) found evidence for radial velocity “bumps” superimposed on the large amplitude binary motion. The residuals after removing the binary orbit yielded a period of about 2.7 years and a velocity semiamplitude of about 25 m s-1. The authors examined several possible causes of the observed radial velocity variability, and finally concluded that the system had a “probable third body.” If substantiated, this would have been the first detection of an extrasolar planetary system. However, Bohlender, Irwin, & Yang (1992) classified $`\gamma `$ Cep as a K0 III star, raising the possibility that the observed radial velocity variations were due to the recently discovered (Walker et al. 1989; Hatzes & Cochran 1993) long-period radial velocity variability of most K giants. A detailed analysis of the CFHT radial velocity data on $`\gamma `$ Cep by Walker et al. (1992) showed a clear low-amplitude (27 m s-1) signal with a 2.52 year period superimposed on the binary orbital motion of indeterminate period. They also found evidence of a possible variation of the Ca II $`\lambda `$8662 Å emission line index with the same 2.5 year period, leading them to conclude that the observed low-amplitude radial velocity variation was most likely due to K-giant variability at the period of the stellar rotation.

Here we present new high-precision radial velocity data on the $`\gamma `$ Cephei system obtained from McDonald Observatory. When combined with the CFHT data, we show that the 2.5-year low-amplitude radial velocity variability of $`\gamma `$ Cep has remained constant for over 20 years. A simultaneous orbital solution for both the binary star and the 2.5-year low amplitude variability is computed. We demonstrate that there is no correlation between the low-amplitude radial velocity variations and the Mt. Wilson Ca II S-index, nor is there any indication of photospheric absorption line profile variability. These results, combined with the current best classification of $`\gamma `$ Cep as a K1 IV star (Fuhrmann 2003), all indicate that the preferred interpretation of the data is that $`\gamma `$ Cep A has a planetary-mass companion in a 2.5 year period orbit.

## 2 The Radial Velocity Data Sets

We consider four independent sets of high precision radial velocity data for $`\gamma `$ Cephei, covering the interval from 1981 through 2002. The first of these is from the CFHT survey Campbell, Walker, & Yang (1988), Walker et al. (1995), with the velocities taken from Table 1 of Walker et al. (1992). All of the other data were from the McDonald Observatory Planetary Search (MOPS) Program (Cochran & Hatzes 2000). Phase I of the MOPS program used the telluric O2 lines near 6300Å as the velocity metric, a technique suggested by Griffin & Griffin (1973). A single order of the 2.7m coudé spectrometer 6-foot camera with the echelle grating was isolated onto a Texas Instruments (TI) $`800\times 800`$ CCD at $`R=210,000`$. This system gave 15–20 m s-1 precision on stars down to about $`V=6`$, but suffered from systematic velocity errors, most likely due to prevailing atmospheric winds. In 1992 the program switched to a temperature stabilized I2 cell (Koch & Wöhl 1984; Libbrecht 1988) as the velocity metric for Phase II of the MOPS ($`R=210,000`$ for this data as well). This eliminated the systematic errors, and gave a routine radial velocity precision of 15 m s-1. This precision was limited by the 9.6 Å bandpass of the spectrum, and by the poor charge-transfer and readout properties of the TI CCD. To solve these problems, and to achieve substantially improved precision, we began Phase III of the radial velocity program in July 1998, using the I2 cell with the newly installed 2dcoudé cross-dispersed echelle spectrograph (Tull et al. 1994). This instrument when used with a Tektronix 2048$`\times `$2048 detector provides a nominal wavelength coverage of 3600 Å – 1$`\mu `$m at a resolving power of $`R=60,000`$. The Tektronix CCD also had significantly better charge transfer and readout properties than the TI device. The complete spectral coverage of 2dcoudé gives us two advantages: first, we can utilize the full reference spectrum of the I2-cell for the RV determination and second, we can simultaneously determine the stellar chromospheric emission in the cores of the Ca II H&K lines to use as stellar chromospheric activity indicators. We will discuss in detail the question whether there is a correlation between the RV-results and the activity indices for $`\gamma `$ Cep in section §4.

To extract the RV-information from the I2 self-calibrated spectra taken during Phase III we employed our $`Austral`$ RV-code, which uses a Maximum Entropy Method deconvolution in order to obtain a higher resolved stellar template spectrum and several reconstruction algorithms which model the shape and symmetry of the instrumental profile (IP) at the time of observation. A detailed description of the $`Austral`$ code can be found in Endl, Kürster, & Els (2000). The algorithm follows in general the modeling idea first outlined by Butler et. al. (1996), and the IP reconstruction techniques by Valenti et al. (1995). Using the 2dcoudé spectrometer in I2-cell self-calibration mode and the $`Austral`$ code for the analysis, we obtain a long term RV precision of $`515\mathrm{m}\mathrm{s}^1`$ on a routine basis for stars down to a magnitude of $`V=9.0`$.

All of the MOPS data from Phases I, II and III are given in Tables 1–3. The uncertainties quoted there for Phase II and III data are the “internal” errors, as represented by the rms of the individual spectral chunks about the mean value. We regard these as a lower limit on the actual uncertainties, since these values do not include the effects of any residual systematic errors that may be present. It is difficult to estimate internal errors for the Phase I data because these are not analyzed in “chunks” like the Phase II and III measurements. Phase I observations of a constant star ($`\tau `$ Cet) show an rms scatter of 23 m s-1. The rms scatter of the Phase I measurements about the final orbital solution is 17.4 m s-1. We list a slightly higher value of 19 m s-1 as the “error” of the Phase I measurements in Table 1. This error was the most reasonable error to assign based upon $`\chi ^2`$ tests of the data and was the value assumed for the Phase I measurements in the final orbital solutions. The true error for an individual Phase I measurement is almost surely higher due to different signal-to-noise ratios of individual spectra and systematic errors due to wind and temperature and pressure changes in the earth’s atmosphere.

The RV measurements for all data sets are shown in Figure 1. A different velocity offset (see below) had to be applied to each data set so that they would all have the same zero point. One can clearly see that the “wiggles” superimposed on the binary variations, that were first reported by Walker et al. (1992) are still present in the final data set (Phase III) taken 20 years later.

Figure 2 shows a sequence of Lomb-Scargle periodograms (Lomb 1976, Scargle 1982) of the RV measurements after subtraction of the velocity variations due to the binary companion (see below). The top panel is only for the CFHT data. The central panel is for the CFHT + McDonald data, excluding the Phase III data. The lower panel is the periodogram for all of the RV measurements. The increase in power at a period of $``$ 2.5 yrs with the addition of each data set is the first indication that these residual RV variations are long-lived and coherent. (The false alarm probability of the peak for the full dataset is $``$ 10-20).

## 3 Orbital Solutions

The program GaussFit (Jefferys et al. 1988; McArthur et al. 1994) was used to fit simultaneously all orbital parameters, both for the stellar and the presumed sub-stellar companion, using nonlinear least squares with robust estimation. Because each instrument produces relative radial velocities with their own arbitrary zero point, the individual velocity offsets were a free parameter in the least squares solution. Table 4 gives the velocity offsets of the different data sets. Subtracting these velocity offsets from an individual data sets will place them all on the same radial velocity scale. The line in Figure 1 shows the combined orbital solution to $`\gamma `$ Cep. Figure 3 shows the “planet only” orbital solution after subtracting the contribution of the binary orbit to the RV measurements.

Table 5 lists the orbital parameters for the the stellar companion. The errors listed are the correlated errors that are produced by GaussFit. The largest error occurs in the period since our observations span about 20 years, or about one-third of the binary orbit. Note that the reduced $`\chi ^2`$ of the binary “only” orbital solution is rather large (4.36) indicating the presence of additional RV variations.

Griffin et al. (2002) presented an orbital solution for the stellar companion to $`\gamma `$ Cep. They combined radial velocity measurements spanning over 100 years taken at five observatories. These measurements included the CFHT data as well as some of our McDonald Observatory measurements (read from an enlarged copy of a published graph!). The parameters of their orbital solution is also listed in Table 5. In spite of the more limited time span of our observations our orbital solution agrees quite well with the Griffin et al. orbit.

Table 6 lists the GaussFit orbital parameters for the sub-stellar companion. Also listed is the rms scatter of the individual datasets about the combined orbital solution. The period of 2.48 years and semi-amplitude of 27.5 m s-1 are consistent with the values found by Walker et al. (1992) using only the CFHT data set. Note that the reduced $`\chi ^2`$ is significantly lower ( $`\chi ^2`$ = 1.47) when including the planet in the orbital solution. Fuhrmann (2003) estimates the primary mass of $`\gamma `$ Cep as $`M`$ $`=`$ 1.59 $`\pm `$ 0.12 $`M_{}`$. This results in a planetary mass of $`M_p`$ $`sin`$$`i`$ $`=`$ 1.7 $`\pm `$ 0.4 $`M_{Jupiter}`$ and an orbital semimajor axis of 2.13 AU.

Figure 4 shows the phase diagram of the individual data sets. The phase, amplitude, and overall shape of the RV curves hold up remarkably well over a 20 year time span, again evidence for long-lived and coherent variations that are consistent with the presence of a planetary companion.

## 4 The Nature of the RV Variations

The question naturally arises as to whether some phenomena instrinsic to the star (spots, convective shifts, pulsations) are responsible for the 2.48 yr period. After all, $`\gamma `$ Cep is most likely a subgiant (see below), a class of stars for which magnetic activity, surface structure, convection, and pulsations are poorly known. Moreover, Walker et al. (1992) did find evidence for weak periodic variations in the equivalent width of the Ca II $`\lambda `$8662 Å line with the same period as the planet. In spite of the intrinsic variability of giant stars planetary companions have been found around other giant stars (Frink et al 2002; Setiawan et al. 2003) and one subgiant (Butler et al. 2001). So planets around giant evolved stars is not an unreasonable hypothesis. Before concluding that the planet hypothesis is the most likely explanation for the residual RV variations we must demonstrate that $`\gamma `$ Cep does not exhibit any significant variations with the planet period in other quantities. Here we examine whether $`\gamma `$ Cep also has spectral and photometric variations.

### 4.1 Bisector Analysis

The spectral line shapes can also provide evidence in support of the planetary hypothesis. Surface features (Hatzes 2002) or nonradial pulsations (Hatzes 1996, Brown et al. 1998) should produce changes in the spectral line shapes with the same period as the RV variations. A convenient means of measuring the asymmetry of a spectral line is the line bisector, or the line segments connecting the midpoints of the spectral line from the core to the continuum. Line bisector studies have been used to establish the planetary nature of 51 Peg (Hatzes, Cochran, & Johns-Krull 1997; Hatzes, Cochran, & Bakker 1998a,b) and the starspot nature of the planet-like RV variations in HD 166435 (Queloz et al. 2001).

The McDonald Observatory Phase I observations taken in the 6300 Å region provide an excellent data set for studying possible line profile variations in $`\gamma `$ Cep. The data has very high spectral resolution ($`R`$ = 210,000) and there are several strong spectral lines that were velocity shifted clear of the telluric features because of the Earth’s barycentric velocity. The Fe I $`\lambda `$6301.5 Å feature was free of telluric lines over the entire data set, Fe I $`\lambda `$6297.8 for 30 observations, and Fe I $`\lambda `$6302.5 Å for 8 observations. Both the bisector velocity span and curvature were measured for these lines. The velocity span is defined as the velocity difference between two arbitrary points on the line bisector while the curvature is the difference between the velocity span of the top half of the bisector minus the velocity span of the lower half of the bisector.

A 5th order polynomial was fit to each measured line bisector and the velocity span points taken at 0.4 and 0.8 of the continuum. (In measuring the line bisector one should avoid both the core and the continuum where the errors in the bisector measurement become large.) The additional point required for the curvature measurements was taken at 0.6 of the continuum value. The mean value of the bisector span and curvature for each line was then subtracted from the individual measurements and the residual span and velocity measurements were then averaged, weighted by the rms scatter of the bisector measurements for each spectral line.

Figure 5 shows the phase-binned averages ($`\mathrm{\Delta }\varphi `$ $``$ 0.05) of the bisector span and curvature variations for $`\gamma `$ Cep phased to the planet orbital period. The short-dashed horizontal line gives the zero-point reference and the long-dashed lines show the velocity extrema of the residual RV variations due to the planetary companion. There are no obvious phase variations of the bisector quantities with the planet period. The least squares sine fit to the bisector quantities assuming a period of 906 days yields amplitudes of 4.8 $`\pm `$ 4.4 m s-1 and 2.4 $`\pm `$ 4.5 m s-1 for the bisector and span variations, respectively. These amplitudes are consistent and are significantly less than the observed 27 m s-1 radial velocity variations of the planet.

The McDonald Phase I data had the largest errors of all the data sets and the phase variations of these show the least evidence for the 906 day period (Fig. 4). One could argue that there are no residual RV variations with which to correlate to the bisector measurements. However, the Phase I data is still consistent with the presence of the planet signal. Not only is this signal present in the Phase I and II data alone, but the power increases by almost a factor of two over the individual periodograms when the two data sets are combined. Although not obvious to the eye, the 906-day period is still present in the Phase I data. On the other hand, the averaged, phase-binned bisector measurements have variations significantly less than the 906-day RV amplitude and we believe this adds additional evidence (along with the photometric and Ca II analysis presented below) in support of the planet hypothesis. If the bisector variations had scatter comparable to the RV amplitude, then one could at least make a plausible argument in favor of possible bisector variability. Figure 5 excludes that.

### 4.2 Photometric Variations

The Hipparcos mission (Perryman 1997) took high precision photometric measurements of $`\gamma `$ Cep contemporaneously with the measurements used for our RV study. If cool spots on the stellar surface, or some form of stellar pulsations were causing the residual RV variations, then this should be evident in the Hipparcos photometry.

Figure 6 shows the Hipparcos photometry taken from 1989.9 to 1993.2 phased to the planet orbital period. Crosses represent individual measurements while solid points represent phase-binned averages ($`\mathrm{\Delta }\varphi `$ $`<`$ 0.05). Error bars represent the rms scatter of the measurements used for each bin. There are no obvious photometric variations with the planet period as confirmed by a periodogram analysis of the daily averages of the Hipparcos photometry (Figure 7). A least squares sine fit to the phased photometric data yields an amplitude of $`\mathrm{\Delta }V`$ = 0.001 $`\pm `$ 0.0009 mag for any possible photometric variations. (A fit to the phase-binned measurements yields $`\mathrm{\Delta }V`$ = 0.0003.) The lack of photometric variations in $`\gamma `$ Cep is also consistent with the planet hypothesis for the residual RV variations for this star.

### 4.3 Ca II Variations

#### 4.3.1 S-index Measurements

Stellar activity in $`\gamma `$ Cep could induce significant periodic centroid shifts in photospheric absorption lines which could be confused for perturbations made by planetary companions (Saar & Donahue 1997; Saar & Fischer 2000, Queloz et al. 2001). For the Phase III data an instrumental setup was chosen so as to include the Ca ii H and K lines on the detector.

To measure stellar chromospheric activity, the Mt. Wilson $`S`$ index was adopted. This index is defined (e.g. Baliunas et al. 1992) as a quantity proportional to the sum of the flux in 1 Å FWHM triangular bandpasses centered on the Ca II H and K lines divided by the sum of the flux in 20Å bandpasses in the continuum at 3901 and 4001Å (Soderblom, Duncan, & Johnson 1991).

The four quantities to be measured (the two calcium line core fluxes plus the two continuum bandpass fluxes) are spread across 3 echelle spectral orders which overlap by several Å. In order to be consistent with activity monitoring with our other programs (Paulson et al. 2002), we did not use measurements in the bluemost order, i.e. the continuum region centered on 3901Å. In addition, we did not measure the Ca II H line (at 3968.47Å) because the wings of strong Balmer H$`ϵ`$ feature (at 3970.07Å) are within the measured Ca ii bandpass. This can adversely affect the measurement of the Ca ii H line flux. Therefore, we have defined an index $`S_{\mathrm{McD}}`$ which is the ratio of the flux in a 1 Å triangular bandpass centered on the Ca ii K line to the flux in a 20 Å bandpass centered on the redward continuum at 4001 Å. Thirty of our program stars (for the McDonald Observatory Planet Search) have previously been measured as part of the Mt. Wilson survey (Baliunas et al. 1995; Duncan et al. 1991). We find a linear relationship between our measurements and those of the Mt. Wilson survey of the form:

$$S_{\mathrm{Mt}.\mathrm{Wilson}}=0.038(\pm 0.006)+1.069(\pm 0.040)\times S_{\mathrm{McD}}$$
(1)
Both the slope and the intercept of the formal fit are within 1 $`\sigma `$ of $`S_{\mathrm{McD}}=S_{\mathrm{Mt}.\mathrm{Wilson}}`$. We are thus able to transform our data into a standard Mt. Wilson $`S`$ index scale using Equation (1).

The periodogram of the S-index measurements (Figure 8) shows no significant power at the orbital frequency of the planet. The false alarm probability of the highest peak, assessed using a bootstrap randomization process (Murdoch et al. 1993; Kürster et al. 1997), is 50%. The lack of variability at the planet period is substantiated by phasing these variations to the 906 day period (Figure 9). Clearly no significant sinusoidal variations in this chromospheric index are present. The transformed Ca II K S-index for $`\gamma `$ Cep is $``$5.3 which implies a level of magnetic activity less than that of the Sun.

#### 4.3.2 Analysis of the Ca II $`\lambda `$8662 Å data from Walker et al. (1992)

A conclusion regarding the planetary nature of the residual RV variations of $`\gamma `$ Cep would not be complete without a discussion of the Ca II variations found by Walker et al (1992). Although the McDonald S-index measurements do not show evidence for rotational modulation, Ca II $`\lambda `$8662 Å equivalent width measurements ($`W_\lambda `$) of Walker et al. (1992) did show a hint of sinusoidal variations when phased to the planetary orbital period. Although this on its own does not completely refute the existence of a planet, having Ca II variations with the same period would cast more doubt on this hypothesis. For these reasons we made a careful examination the significance of the Ca II variations reported by Walker et al. (1992).

Figure 10 shows the correlation between the RV and the changes in the Ca II $`\lambda `$8662 Å equivalent width. The two quantities show no obvious correlations. The correlation coefficient, $`r`$, is only 0.08 and the probability that the two quantities are not correlated is 0.52. However, there are 3 obvious outliers in the figure ($`|\mathrm{\Delta }W_\lambda |>`$ 3%). Eliminating these 3 points increases the correlation coefficient to $`r`$ = 0.23 with a probability of 0.08 that the quantities are uncorrelated. However, this is still not a strong correlation.

The McDonald S-index data also do not seem to be correlated to the RV measurements. Figure 11 shows the correlation between the S-index and radial velocity measurements from the McDonald data. The correlation coefficient in this case is $`r`$ = $``$0.3 and the probability that the two quantities are uncorrelated is 0.13.

Figure 12 shows the Lomb-Scargle periodogram of the variations in the Ca II $`\lambda `$8662 Å equivalent width measurements of Walker et al. (1992). (There was a clear outlier in the Walker et al. Ca II data with $`\mathrm{\Delta }EW`$ = $``$4.3. This was also evident in Fig. 3 of Walker et al. (1992). This data point was eliminated before performing the periodogram analysis.) Although the most power occurs at a periods of $``$ 15 days (top panel), there appears to be significant power at the planet period in the expanded scale in the lower panel (the vertical line indicates the location of the planet RV frequency). This was the basis for the Walker et al. favoring rotational modulation as the cause of the RV variations.

A more detailed examination of the CFHT Ca II $`\lambda `$8662 Å equivalent measurements show that the variation in these are not long-lived and may have no relation to the observed 906-day RV period. We divided the CFHT Ca II data into two data sets, one spanning 1981–1986.5 and the second 1986.5–1992. The division of the data set was taken so as to maximize the power in the long period variations found in the second set. (Our results do not change substantially if we divide the data set into two equal numbered points.) This resulted in 21 data points in the first data set and 29 in the second set.

Figure 13 shows the periodograms of these two data sets. The 1981-1986.5 set shows no power in the frequency interval $`0<`$ $`\nu `$ $`<`$ 0.01 c d-1. The 1986.5–1992 data set shows significant power, but at a frequency corresponding to a period of 781 $`\pm `$ 116 days, significantly less than the planet orbital period. Extending the periodogram to higher frequencies shows that this is the highest peak out to the Nyquist frequency, in contrast to the case for the full data set where the low frequency feature was the second highest peak. A periodogram analysis of the full Ca II data set lowers the Lomb-Scargle power near $`\nu `$ = 0.011 c d-1 from about 9.0 to 6.7 and increases the period slightly. Unlike the case for the RV data, increasing the number of measurements does not increase the power in the periodogram of the feature of interest. Clearly this signal is not long-lived.

The statistical significance of this signal was examined using the bootstrap randomization technique. The Ca II $`\lambda `$8662 Å equivalent width measurements over the time span 1986.5 – 1992 were randomly shuffled keeping the observed times fixed. A periodogram was then computed for each “random” data set. The fraction of the periodograms having power higher than the data periodogram in the range 0.0005 $`<`$ $`\nu `$ $`<`$ 0.01 c d-1 is the false alarm probability that noise would create the detected signal. After 105 shuffles there was no instance of a random periodogram having power higher than that found in the data periodogram. The false alarm probability is thus $`<`$ 10-5. This signal, in spite of only being present in the data for the last half of the data set, is highly significant.

The period of the Ca II $`\lambda `$8662 Å measurements is significantly less, by 1$`\sigma `$, from the residual RV period. Figure 14 shows the phase diagrams of the Ca II $`\lambda `$8662 Å measurements. The top panel are the measurements during 1986.5-1992 phased to the 781 day period found in the Ca II data. The central panel is the same data phased to the 906 day residual RV period. Although there are slight sinusoidal variations when phased to the 906 day period, phasing to the 781 day period produces significantly less scatter. By comparison, the lower panel shows the Ca II measurements from 1981 – 1986.5 phased to the 906 day period. Crosses are individual measurements whereas solid points represent phase-binned averages (the error bar represents the scatter of the data used for each phase bin). There are no significant sinusoidal variations in the 1981 – 1986.5 data set.

A least squares sine fit to the the phased Ca II yields an amplitude of 1.67 $`\pm `$ 0.25 mÅ for the dates spanning 1986.5 – 1992 (1.03 $`\pm `$ 0.25 assuming a 906 day period). For the dates covering 1981 – 1986.5 the amplitude of the $`EW_{8662}`$ variations is 0.48 $`\pm `$ 0.33 mÅ. Clearly there are amplitude variations in the Ca II $`\lambda `$8662.

We conclude that the periodic variations in Ca II found by Walker et al. (1992) are not long-lived and were only present during 1986.5–1992. This period was clearly not present in our Ca II S-index measurements spanning 1998–2002. Since the RV variations are unchanged during the entire time span 1981–2002 we do not believe that the 906 day RV period is related to the cause of the Ca II $`\lambda `$8662 Å variability.

### 4.4 The Spectral Classification of $`\gamma `$ Cep

One reason that Walker et al. (1992) favored rotational modulation for the residual RV variations was the reclassification of $`\gamma `$ Cep as a K0III. This was based on a visual comparison of the spectrum of $`\gamma `$ Cep to other K giants (Bohlender et al. 1992). Recently, Fuhrmann (2003) presented an extensive spectral analysis of nearby stars in the galactic disk and halo. Included in this study was $`\gamma `$ Cep. Table 7 list the stellar properties derived by Fuhrmann. Also listed is the Hipparcos distance ($`d_{Hip}`$) and distance determined spectroscopically ($`d_{Sp}`$) using the derived stellar parameters. These distances agree to within 3% indicating that we can have some confidence in the results of the spectral analysis. Fuhrmann’s classification of K1IV for $`\gamma `$ Cep is consistent with the sub-giant status for $`\eta `$ Cep (K0IV) which has a comparable effective temperature (4990 K), gravity (log $`g`$ = 3.4), radius (4.14 $`R_{}`$), and bolometric magnitude ($`M_{bol}`$ = 2.3). The most current and best analysis of $`\gamma `$ Cep supports this being a subgiant star.

## 5 Searching for additional companions

The long time-baseline of RV-monitoring of $`\gamma `$ Cep and the emergence of several multiple extrasolar planetary systems from Doppler surveys (e.g. Butler et al. 1999), encouraged us to search for additional companions in this system. For this purpose we performed a period search within the RV-residuals, after subtracting both the binary orbital motion and the first planetary signal.

Fig. 15 shows the Lomb-Scargle periodogram of the RV-residuals in the period range of $`2`$ to $`7900`$ days (frequency = 0.00013 – 0.5 c d-1). The strongest peak is found at $`P11`$ days but its significance level is very low (false-alarm-probability (FAP) is $`6.5\%`$). The FAP-levels shown in the figure were determined by $`10,000`$ runs of a bootstrap randomization scheme.

It is clear from this analysis that no additional periodic signal above the noise level is present in the RV-residuals of $`\gamma `$ Cep and we thus conclude that the 906-day period planet is the only giant planet in this system evident in our data..

## 6 Discussion

Precise stellar radial velocity measurements of $`\gamma `$ Cep from 4 independent data sets spanning over 20 years show long-lived, low-amplitude RV variations ($`K`$ = 27.5 m s-1) superimposed on the larger radial velocity variations due to the reflex motion caused by a stellar companion. We interpret the low amplitude, shorter period variations as due to the presence of a planetary companion with $`M`$ $`sin`$$`i`$ = 1.7 $`M_{Jupiter}`$ and an orbital semi-major axis, $`a`$ = 2.13 AU. Our conclusion that these short period variations are not due to rotation, pulsations, or changes in the convection pattern of the star is based on several facts:

1. The 2.48 yr period has been present for over 20 years with no changes in phase or amplitude during this time.
2. No variations with this period are seen in the McDonald Ca II S-index measurements.
3. Spectral line bisector span and curvature measurements for $`\gamma `$ Cep that are constant to less than 5 m s-1 over an orbital cycle of the planet.
4. Contemporaneous Hipparcos photometry that is constant to less than 0.001 mag over an orbital cycle.
5. The periodic variations in the Ca II $`\lambda `$8662 Å found by Walker et al. (1992) were only present during 1986.5–1992 and are thus most likely not associated with the residual RV variability observed for this star.

The CFHT Ca II data do show evidence for long period variability. Although this signal is weakly detected in the full data set, it is much stronger in only the last half of the data set spanning 1986.5–1992, and completely absent in the data from 1981 – 1986.5. Furthermore, our S-index measurements made during 1997 – 2003 show no modulation. If rotational modulation was responsible for the residual RV variations then we would have seen RV amplitude changes. However, the amplitude and phase for the residual RV variations remained constant. It therefore seems unlikely that rotational modulation is the cause of the 906 day period. Furthermore, the best fit period to the long term Ca II variations during 1986.5 – 1992 is 781 $`\pm `$ 116 days, which is significantly less than the 906 day RV period.

Gamma Cep shows no Ca II modulation during epochs centered on 1994 and 2000, and modulation (781 day period) during a time centered on 1989. This indicates a possible “activity cycle” period of 10 – 15 years, much longer than the presumed planet period.

Our spectral line bisector measurements, which are constant, also support our conclusion that we are not seeing rotational modulation. One could argue that the residual RV variations are due to changes in the convection pattern in the star and not due to magnetic structure (plage, spots). For instance, if the ratio of areas of convective, hot rising cells and intergranule, sinking lanes changes with an activity cycle (in this case with a period of $``$ 900 days), then the amount of convective blue (red) shift would change periodically resulting in a measured RV signal. If the stellar magnetic fields are strong enough to alter the convection pattern of the star, but too weak to cause chromospheric structure, then we would see RV variations without strong variations in the Ca II emission. However, in this case we should also see changes in the spectral line bisectors. The lack of significant variability in the spectral line bisectors does not support the hypothesis of a changing convection pattern on the star, or a least changes that can influence the RV measurements.

The inclination angle of the rotation axis can be estimated by comparing the expected equatorial rotational velocity to the projected rotational velocity (v sin $`i`$), assuming that 781 days (period of the Ca II $`\lambda `$8662 Å variations) is the true rotation period of $`\gamma `$ Cep. Fuhrmann (2003) estimates a radius $`R`$ = 4.66 $`R_{}`$ which gives an equatorial rotational velocity of 4.9 km s-1. We have measured a v sin $`i`$ $`=`$ 1.5 $`\pm `$ 1.0 km s-1 which is consistent to the value determined by Fuhrmann (2003). This yields a sin $`i`$ = 0.5 – 0.1. Assuming that the orbital and rotation axis of the star is aligned, the true mass for the planetary companion to $`\gamma `$ Cep is $``$ 3–16 $`M_{Jupiter}`$.

Given that the 906-day RV variations are due to a sub-stellar companion, the question arises whether this system is stable. After all, both the planet and the stellar companion have modestly eccentric orbits. Preliminary indications are that this system is indeed stable (Dvorak et al. 2003). (The study by Dvorak et al. was made prior to the publication of our paper and it used preliminary orbital parameters published only in conference abstracts. The stability analysis should be re-done using the final parameters and error estimates we have presented here.)

The $`\gamma `$ Cep system presents a very interesting system for the study of planet formation. Planets have been found around host stars that are in a binary system, but these are widely separated pairs so the presence of the stellar companion may not have an influence on the planet formation around one of the stars. The $`\gamma `$ Cep binary is the shortest period binary for which a planetary companion has been found in orbit around one of the components. This implies that even in such a relatively close binary the presence of the stellar companion does not hinder the process of planet formation.

One final comment. RV surveys have had a stunning success at finding planets in orbit around other stars. Yet this is still an indirect detection method with the disadvantage that various stellar phenomenon can mimic a planetary signal. RV searches for extrasolar planets are, and should continue to be a careful process with many candidates undergoing the painstaking process of confirmation. Several extrasolar planet candidates have a proved to be due to stellar phenomena (Queloz et al. 2001; Butler et al. 2002) and more false detections will undoubtedly be uncovered in the future. Likewise, radial velocity variations that were initially attributed to stellar phenomena in the past may ultimately prove to be due to a planetary companion, as is apparently the case for $`\gamma `$ Cep. RV searches for extrasolar planets must be accompanied by careful studies of the host star and the measurement of a various photometric and spectroscopic (Ca II emission, line shapes) quantities. These measurements are important for they could have confirmed the planet around $`\gamma `$ Cep seven years prior to the discovery of 51 Peg b.

This material is based upon work supported by the National Aeronautics and Space Administration under Grant NAG5-9227 issued through the Office of Space Science, and by National Science Foundation Grant AST9808980. We thank Sebastian Els for his attempt to image the stellar secondary using adaptive optics observations with the 4.2m-William Herschel telescope.