text:
# 1 Introduction

## 1 Introduction

Many SUSY studies rely on computer codes that calculate the mass spectrum of the minimal supersymmetric standard model (MSSM), the couplings, branching ratios, etc., from given sets of model parameters. For the LHC, for instance, many simulations are done for particular benchmark scenarios or by mapping the $`(m_0,m_{1/2})`$ parameter plane. For such studies it is certainly important wheather a particular decay channel is open or not and what branching ratio it has. Also, theoretically or experimentally excluded regions depend on the details of the spectrum. Studies for an $`e^+e^{}`$ Linear Collider deal, in addition, with high precision measurements of (s)particle properties, with the determination of the underlying SUSY breaking parameters, their extrapolation to the GUT scale, model destinction, etc. Experimental accuracies of the per-cent or even per-mille level are expected. It is thus clear that we need theoretical predictions of a precision comparable to the experimental accuracy. However, it has been noticed  that different programs can give quite different results for the same set of input parameters.

In this article, we compare the mass spectrum calculations of four public codes: Isajet 7.63 , SuSpect 2.005 , SoftSusy 1.4 , and SPheno 1.0 , in the minimal supergravity (mSUGRA) framework. We discuss the renormalization group (RG) running and the implementation of radiative corrections, concentrating on the parameter regions where the largest differences are encountered: large $`\mathrm{tan}\beta `$ and large $`m_0`$. An overview of which corrections are implemented in each of the four programs is given in Table 1.

## 2 Large $`\mathrm{𝐭𝐚𝐧}𝜷`$

Large $`\mathrm{tan}\beta `$ has always been recognized as a difficult case since it requires a thorough treatement of the bottom Yukawa coupling $`h_b`$. It is well known  that $`h_b`$ gets large $`\mathrm{tan}\beta `$ enhanced corrections from SUSY loops, the dominant contributions coming from $`\stackrel{~}{b}\stackrel{~}{g}`$ and $`\stackrel{~}{t}\stackrel{~}{\chi }^+`$ exchanges. These generate a $`H_2^0b\overline{b}`$ coupling, which is forbidden at tree-level, $`h_bH_1^0b\overline{b}+\mathrm{\Delta }h_bH_2^0b\overline{b}`$. This modifies the tree-level relation between the bottom mass and Yukawa coupling, $`m_b=h_bv_1m_b=h_bv_1(1+\mathrm{\Delta }_b)`$ with $`\mathrm{\Delta }_b=(\mathrm{\Delta }h_b/h_b)\mathrm{tan}\beta `$. In the programs under discussion this is taken into account as

$$h_b(M_Z)=\widehat{m}_b^{\mathrm{MSSM}}(M_Z)/v_1(M_Z),\widehat{m}_b^{\mathrm{MSSM}}(M_Z)=\frac{\widehat{m}_b^{\mathrm{SM}}(M_Z)}{1+\mathrm{\Delta }m_b/m_b}.$$
(1)
Here $`\widehat{m}_b^{\mathrm{SM}}`$ is the $`\overline{\mathrm{DR}}`$ bottom mass in the Standard Model and $`\mathrm{\Delta }m_b=(\mathrm{\Delta }m_b)^{\stackrel{~}{b}\stackrel{~}{g}+\stackrel{~}{t}\stackrel{~}{\chi }^++\mathrm{}}`$ contains the SUSY-loop corrections. The complete 1-loop expression for $`\mathrm{\Delta }m_b`$ is given in .222Here note that  and  partly have different conventions, e.g., for the ordering of the squark mass eigenstates and the sign of $`\mu `$. Compared to the naive 1-loop expansion $`\widehat{m}_b^{\mathrm{MSSM}}=\widehat{m}_b^{\mathrm{SM}}(1\mathrm{\Delta }m_b/m_b)`$, eq. (1) makes a numerical difference of about $`10\%`$ in $`h_b`$ and about 10–30% in $`m_A`$ for large $`\mathrm{tan}\beta `$. The re-summation of SUSY threshold corrections  will be discussed elsewhere . Although all four programs now apply eq. (1), some numerical differences in $`h_b`$ remain. These are partly due to differences in $`\alpha _s`$: Suspect, SoftSusy and SPheno calculate $`\alpha _s`$ in the $`\overline{\mathrm{DR}}`$ scheme, Isajet uses the $`\overline{\mathrm{MS}}`$ value. Another reason is that Isajet uses $`m_b=m_b(M_{SUSY})`$ for the expression $`\mathrm{\Delta }m_b/m_b`$ in eq. (1), while the other programs use $`m_b(M_Z)`$ or the bottom pole mass; also the gluino masses differ by about 5%. Moreover, the vacuum expectation values $`v_{1,2}`$ are not running in Isajet.

The bottom Yukawa coupling has its largest effect in the Higgs sector. Figure 1 shows the running of $`m_{H_{1,2}}^2`$ for $`m_0=400`$ GeV, $`m_{1/2}=300`$ GeV, $`A_0=0`$, $`\mu >0`$, and the two cases $`\mathrm{tan}\beta =10`$ and $`\mathrm{tan}\beta =50`$. As one can see, there is good agreement for not too large $`\mathrm{tan}\beta `$. However, for $`\mathrm{tan}\beta =50`$, quite different results are obtained for $`m_{H_1}^2`$, whose evolution is driven by $`h_b`$:

$$\frac{dm_{H_1}^2}{dt}\frac{3}{8\pi ^2}h_bX_b+\mathrm{},X_b=(m_{\stackrel{~}{Q}}^2+m_{\stackrel{~}{D}}^2+m_{H_1}^2+A_b^2).$$
(2)
Note in particular the dotted line which shows the result obtained with Isajet 7.58. In this version, the SUSY corrections to $`h_b`$ were not yet resummed.

The differences in $`m_{H_1}^2`$ directly translate into $`m_A^2`$ and thus into the physical Higgs boson masses, since

$$m_A^2=\frac{1}{c_{2\beta }}\left(\overline{m}_{H_2}^2\overline{m}_{H_1}^2\right)+\frac{s_\beta ^2t_1}{v_1}+\frac{c_\beta ^2t_2}{v_2}M_Z^2.$$
(3)
Here $`\overline{m}_{H_i}^2=m_{H_i}^2t_i/v_i`$, $`i=1,2`$, and $`t_{1,2}`$ are the tadpole contributions. The self energies of $`Z`$ and $`A`$ have been neglected in eq. (3). We note that including only the tadpoles from the third generation is in general a good approximation. The remaining 1-loop contributions account for a $`𝒪(1\%)`$ correction.

Figure 2 shows the Higgs boson masses obtained by the four programs as a function of $`\mathrm{tan}\beta `$. The new Isajet version 7.63 has led to a majer improvement compared to the situation discussed in  (the results obtained by Isajet 7.58 are again shown as dotted lines in Fig. 2). For $`m_A`$ and $`m_{H^\pm }`$ there is now agreement within $`10\%`$ up to $`\mathrm{tan}\beta 45`$. Sources for the remaining differences are pointed out above. Moreover, it makes a difference whether one uses running couplings and/or masses for the tadpoles $`t_{1,2}`$. Here each program has a different approach. For the neutral scalars, however, the situation is not so good. Especially for $`m_{h^0}`$, a discrepancy of $`4`$ GeV is too large compared to the expected experimental accuracy. This discrepancy is mainly due to the different radiative corrections taken into account for the ($`h^0,H^0`$) system. They vary between 1- and 2-loop, effective potential and diagramatic calculations, see Table 1. Given the expected experimental accuracy for $`m_{h^0}`$ it is clear that the best available calculation should be used.

## 3 Large $`𝒎_\mathrm{𝟎}`$

For large $`m_0`$, the running of $`m_{H_2}^2`$ becomes very steep and very sensitive to the top Yukawa coupling $`h_t=\widehat{m}_t/v_2`$ :

$$\frac{dm_{H_2}^2}{dt}\frac{3}{8\pi ^2}h_tX_t+\mathrm{},X_t=(m_{\stackrel{~}{Q}}^2+m_{\stackrel{~}{U}}^2+m_{H_2}^2+A_t^2).$$
(4)
As a result, the $`\mu `$ parameter given by

$$\mu ^2=\frac{\overline{m}_{H_1}\overline{m}_{H_2}^2\mathrm{tan}^2\beta }{\mathrm{tan}^2\beta 1}\frac{1}{2}M_Z^2$$
(5)
becomes extremely sensitive to $`h_t`$. This is visualized in Fig. 3 where we show in (a) the running of $`m_{H_{1,2}}^2`$ for $`m_0=1450`$ GeV, and in (b) $`\mu `$ as a function of $`m_0`$. The other parameters are $`m_{1/2}=300`$ GeV, $`A_0=0`$, $`\mu >0`$, and $`\mathrm{tan}\beta =10`$. The large discrepancy in $`\mu `$ for $`m_0\stackrel{>}{}\mathrm{\hspace{0.33em}800}`$ GeV lead to completely different chargino/neutralino properties and likewise to very different excluded regions in Isajet compared to the other programs. For instance, radiative EWSB breaks down in Isajet for $`m_0\stackrel{>}{}\mathrm{\hspace{0.33em}1.5}`$ TeV. In SuSpect, SoftSusy, and SPheno, this happens only for $`m_0\stackrel{>}{}\mathrm{\hspace{0.33em}2.5}`$–2.8 TeV.333After the conference, a sign error was corrected in SPheno. As a consequence, its results for large $`m_0`$ now nicely agree with those of SoftSusy and SuSpect (contrary to what was presented in the talk).

In order to understand the behaviour in Fig. 3b it is useful to write eq. (5) in the form

$$\mu ^2c_1m_0^2+c_2m_{1/2}^20.5M_Z^2.$$
(6)
Approximate analytical expressions for $`c_1`$ and $`c_2`$ can be found e.g., in . For $`A_0=0`$ and $`\mathrm{tan}\beta =10`$ we get

$$c_1\left(\frac{\widehat{m}_t}{156.5\mathrm{GeV}}\right)^21,c_2\left(\frac{\widehat{m}_t}{102.5\mathrm{GeV}}\right)^20.52.$$
(7)
Since the Higgs potential is minimized at $`M_{SUSY}=\sqrt{m_{\stackrel{~}{t}_1}m_{\stackrel{~}{t}_2}}`$, we take $`\widehat{m}_t`$ in eq. (7) as $`\widehat{m}_t=\widehat{m}_t(M_{SUSY})`$. The $`m_0`$ dependence seen in Isajet is reproduced for $`\widehat{m}_t151`$ GeV. The one of SuSpect, SoftSusy and SPheno is reproduced for $`\widehat{m}_t155`$ GeV. Figure 4 shows a contour plot of $`\mu `$ in the ($`m_0,\widehat{m}_t`$) plane. Notice the fast increasing dependence on $`\widehat{m}_t`$ for increasing $`m_0`$. Notice also that for $`\widehat{m}_t156`$–157 GeV, $`\mu `$ becomes almost independent of $`m_0`$, which is the actual focus point condition.

There are some obvious differences in the calculations. For instance, $`M_{SUSY}`$, the scale where the SUSY parameters are frozen out and the Higgs potential is minimized, varies by about 100 GeV due to different radiative corrections to the stop masses, c.f. Table 1. In the loop corrections to $`m_t`$, analogous differences occur as discussed above for $`\mathrm{\Delta }m_b/m_b`$. Also the evolution of $`h_t`$ between $`M_Z`$ and $`M_t`$ and the inclusion of threshold effects are delicate points. However, this is not yet sufficient to explain the observed discrepancies. More work is needed to clarify the situation.

## 4 Conclusions

For the calculation of the SUSY mass spectrum from GUT scale boundary conditions, there are two particular difficult parameter regions where large numerical differences have been noticed: large $`\mathrm{tan}\beta `$ and large $`m_0`$. These regions are very sensitive to the bottom and top Yukawa couplings, respectively.

The inclusion of the SUSY 1-loop corrections to $`h_b`$ has led to a considerable improvement in the large $`\mathrm{tan}\beta `$ case. In particular, the four programs now agree on $`m_A`$ within $`\stackrel{<}{}\mathrm{\hspace{0.33em}10}\%`$ for $`\mathrm{tan}\beta \stackrel{<}{}\mathrm{\hspace{0.33em}45}`$. Further improvements are of course desirable.

For large $`m_0`$, there are still very large numerical discrepancies due to the corrections to $`h_t`$. As a matter of fact, $`h_t`$ is much smaller in Isajet than in the other programs. Some differences in the calculation of $`h_t`$ have been pointed out, but these do not satisfingly explain the observed discrepancies. Work is in progress to clarify the situation .

## Acknowledgements

We thank H. Baer, A. Djouadi, and J.-L. Kneur for useful discussions.


format: text

text:
# GRB 010921: Strong Limits on an Underlying Supernova from HST

## 1 Introduction

Since the discovery of gamma-ray burst (GRB) afterglows there has been growing evidence linking GRBs to massive stars: the host galaxies of GRBs are star-forming galaxies and the position of GRBs appear to trace the blue light of young stars (Bloom, Kulkarni & Djorgovski, 2002); some of the host galaxies appear to be dusty with star-formation rates comparable to ultra-luminous infrared galaxies (Berger, Kulkarni & Frail, 2001; Frail et al., 2002). On smaller spatial scales, there is growing evidence tying GRBs to regions of high ambient density (Galama & Wijers, 2001; Harrison et al., 2001) and the so-called dark GRBs arise in or behind regions of high extinction (Djorgovski et al., 2001; Piro et al., 2002).

However, the most direct evidence linking GRBs to massive stars comes from observations of underlying supernovae (SNe) and X-ray lines. The presence of X-ray lines would require a significant amount of matter on stellar scales (e.g. Piro et al. 2000), as may be expected in models involving the death of massive stars. However, to date, these detections (e.g. Piro et al. 2000; Reeves et al. 2002) have not been made with high significance.

If GRBs do arise from the death of massive stars, then it is reasonable to expect associated SNe. The GRB-SN link was observationally motivated by two discoveries: the association of GRB 980425 with the peculiar Type Ic SN 1998bw (Galama et al., 1998; Kulkarni et al., 1998) and an excess of red light superposed on the rapidly decaying afterglow of GRB 980326 (Bloom et al., 1999). However, these two discoveries were not conclusive. The SN association would require GRB 980425 to be extra-ordinarily under-energetic as compared to all other cosmologically located GRBs and the case for GRB 980326 is weakened by the lack of a redshift for the GRB or the host galaxy.

Nonetheless, the two discoveries motivated searches for similar underlying SN components. As summarized in section 5, suggestions of similar red “bumps” in the light curves of various other GRB afterglows have been made (to varying degrees of confidence).

However, there is little dispute that the well-studied red bump in the afterglow of GRB 011121 is most easily explained by an underlying supernova (Bloom et al., 2002; Garnavich et al., 2002). Furthermore, from radio and IR observations of the afterglow (Price et al., 2002a), there is excellent evidence that the circumburst medium was inhomogeneous with ambient density $`\rho r^2`$, as expected from a massive star progenitor (Chevalier & Li, 2000); here, $`r`$ is the distance from the progenitor.

These developments are in accordance with the expectation of the “collapsar” model (Woosley, 1993; MacFadyen & Woosley, 1999). In this model, the core of a rotating massive star collapses to a black hole which then accretes matter and drives a relativistic jet. Internal shocks within this jet first cause bursts of $`\gamma `$-rays and then subsequently result in afterglow emission as the jet shocks the ambient medium.

It is important to appreciate that the SN light is primarily powered by radioactive decay of the freshly synthesized 56Ni whereas the burst of $`\gamma `$-rays are powered by the activity of the central engine. In the current generation of collapsar models, there is sufficient flexibility to allow for a large dispersion of 56Ni and the energy of the engine. Thus, the next phase of understanding the GRB-SN connection111A class of models, known as “supranova” models, posit a supernova greatly in advance, many months, of the the GRB event (Vietri & Stella, 1999). The long delay was physically motivated to explain the X-ray lines as arising from a large spatial region. The current data (e.g. GRB 011121) do not allow for such long delays. will benefit from (and require) observational measures of these parameters.

Motivated thus, we have an ongoing program of searches for SNe in GRB afterglows with the Hubble Space Telescope (HST). Here, we present a systematic search for a SN underlying GRB 010921. In §2 we present our observations and the details of photometry in §3. We fit afterglow models and constrain the brightness of an underlying SN in §4. We then present an overview of previous such efforts and conclude in §5.

## 2 Observations and Reductions

GRB 010921 was detected by the High Energy Transient Explorer (HETE-2) satellite at 2001 September 21.219 UT (Ricker et al., 2002) and the position was refined by the InterPlanetary Network error-box (Hurley et al., 2001). Using the 5-m Hale telescope and the Very Large Array we discovered the afterglow of this event as well as the redshift of the host galaxy (Price et al., 2002b).

The low redshift of this event, $`z=0.451`$, made it a prime candidate for a search for an underlying SN. Accordingly, as a part of our large Hubble Space Telescope (HST) Cycle 9 program (GO-8867, P. I.: Kulkarni), we triggered a series of observations with the Wide Field Planetary Camera 2 (WFPC2) aboard HST. Owing to the lateness in identifying the afterglow candidate, the first observation was on day 35, slightly after the expected peak of the SN. At each of epochs 1–3 we obtained $`4\times 400`$ s exposures in each of five filters (F450W, F555W, F702W, F814W and F850LP) with a single diagonal dither by 2.5 pixels to recover the under-sampled point-spread function (PSF). The fourth epoch was optimized for photometry of the host galaxy and, accordingly, we increased the exposure time to $`4\times 1100`$ s.

We used “On-The-Fly” pre-processing to produce debiased, flattened images. The images were then drizzled (Fruchter & Hook, 2002) onto an image with pixels smaller than the original by a factor of 0.7 using a pixfrac of 0.8. After rotation to a common orientation the images were registered to the first epoch images using the centroids of common objects in the field. The typical r.m.s. registration errors were less than 0.15 drizzled pixels.

## 3 Host Subtraction and Photometry

The host galaxy of GRB 010921 has an integrated magnitude of $`R22`$ mag or about 5$`\mu `$Jy (Price et al., 2002b). Consequently great care has to be taken to properly photometer the fading afterglow. Below, we review various photometric techniques.

Total magnitudes: The simplest technique is to perform aperture photometry (e.g. Galama et al. 2000; Price et al. 2001). The afterglow flux is obtained by subtracting the host flux estimated from a very late time measurement. A major concern is that the host flux is dependent upon the choice of aperture (both center and size). Thus, if different images have different seeing then it is possible to obtain an artificial bump in the light curve.

Host subtraction: The above concern can be alleviated by subtracting a late-time image from the earlier images. The afterglow may then be easily photometered in the host-subtracted images. This method has been used with considerable success by those observing SNe Ia (e.g. Schmidt et al. 1998).

$`N\times (N1)/2`$ subtraction: In this technique, each image is subtracted from every other image and the afterglow residual photometered. The flux at each epoch can be fit through least-squares, assuming the flux at the final epoch is zero (Novicki and Tonry, personal communication). This method makes use of the fact that the host galaxy has not been observed only once at late times, but at each epoch and thus better S/N can be obtained from the over-constrained system.

We employed the $`N\times (N1)/2`$ subtraction technique to photometer the GRB 010921 afterglow in our HST images. The images were subtracted using a modified version of ISIS (Alard, 2000) and photometered using the analytic PSF-fitting routine within Vista (J. Tonry, personal communication). We used the synphot package within IRAF to calculate the response of the instrument and filter combination to a source with constant flux of 1 mJy; the resulting values are AB magnitudes (Fukugita et al., 1996), expressed as fluxes. Corrections were made for charge-transfer (in)efficiency (CTE) using the prescription of Dolphin (2000) and aperture-corrected to infinity.

We have also re-analyzed and photometered ground-based images (Price et al., 2002b; Park et al., 2002) of the afterglow, applying $`N\times (N1)/2`$ subtraction. Since this technique assumes that the flux of the afterglow in the final epoch is zero, which may not be correct for these images, we subtracted the appropriate fourth-epoch HST observation (which we have assumed contains no afterglow) from the final ground-based images, measured the flux of the afterglow and added this value to the fluxes derived from the $`N\times (N1)/2`$ subtraction.

The results of the photometry are host-subtracted fluxes for the afterglow in each of the images, under the assumption that the afterglow flux in the final HST image (2001 Dec 21) is zero (or negligible). These values are presented in Tables 1 and 2. The values in Table 2 supersede the corresponding measurements presented in Price et al. (2002b) and Park et al. (2002). We plot the afterglow light curves in Figure 1. The light-curves are monotonically decreasing (i.e. do not level off), and hence we deduce that our assumption of negligible flux in the final HST image is justified.

## 4 Discussion

Temporal breaks in optical light-curves have been seen in many afterglows and are usually attributed to a “jet” geometry (see Frail et al. 2001). We adopt a standard optical afterglow model, consisting of a broken power-law temporal decay with power-law indices $`\alpha _1`$ and $`\alpha _2`$, and a power-law spectral index, $`\beta `$ (Sari, Piran & Narayan, 1998; Sari, Piran & Halpern, 1999).

Each of these indices are functions of the electron energy distribution index, $`p`$, dependent upon the location of the cooling break relative to the optical bands, and so we consider two cases: the cooling break is redward of the optical (hereafter, case R); and the cooling break is blueward of the optical (case B). We consider, in addition to a constant circumburst medium, an inhomogeneous circumburst medium, $`\rho r^2`$ (see Chevalier & Li 2000 and Price et al. 2002a).

We apply the parametric extinction curves of Cardelli, Clayton & Mathis (1989) and Fitzpatrick & Massa (1988) using the interpolation calculated by Reichart (2001). These extinction curves are characterized by two values, the magnitude of the extinction in the rest-frame of the host galaxy, $`A_V^{\mathrm{host}}`$, and the slope of the UV extinction curve, $`c_2`$ (see Reichart 2001). Following Price et al. (2002b), we adopt $`c_2=4/3`$, corresponding to an LMC-like extinction curve. Adopting other extinction curves (e.g. MW, SMC) yields similar, but more-constraining results (i.e. any underlying SN must be even fainter than the upper limit we derive below); see Price et al. (2002b).

The main purpose of this analysis is to determine whether the light curves contain an SN component. To this end, we use the observations of SN 1998bw for an SN template since it is one of the well observed bright Ib/c SNe which may be related to a GRB (Galama et al., 1998). Specifically, we used the $`UBVRI`$ photometry of Galama et al. (1998) and derived the flux distribution of SN 1998bw, using the zero-points and filter curves of Bessell (1990). The resulting low resolution spectrum (consisting of 5 points at the effective wavelength of each broadband filter), is redshifted to $`z=0.451`$ (Price et al., 2002b), assuming a flat lambda cosmology with $`\mathrm{\Omega }_\mathrm{M}=0.3`$ and $`H_0=65`$ km s-1 Mpc-1. The redshifted spectrum, which represents what SN 1998bw would look like at cosmological distances, is integrated with the appropriate filter curve to derive the apparent brightness at this redshift.

SN 1998bw at $`z=0.451`$ would peak in the rest-frame $`I`$-band at approximately 4 $`\mu `$Jy. It is evident from Figure 1 that the afterglow is much fainter than this, and, further, that there is no clear bump in the afterglow light curve. We therefore allow the SN component to be scaled by $`\delta `$ magnitudes in our model. The SN is placed behind the same foreground (i.e. Milky Way) and host galaxy extinction as the afterglow (which can be inferred by demanding that the temporal and intrinsic spectral slopes, which both depend on the electron distribution index, $`p`$, be consistent; see e.g. Price et al. 2002a).

To calculate the SN detection limit of our observations, we fit the model by minimizing $`\chi ^2`$. The afterglow was not detected in any of the F450W images, and so we exclude them from our analysis. Subtracting the host F450W image from our ground-based $`g^{}`$ image left a large residual at the position of the host galaxy (not of the OT). This poor subtraction is likely due to the filter mis-match, and so we do not include this point in our analysis.

Our analyses are summarized in Table 3. In short, we find no evidence for an underlying SN. In order to calculate the formal limits, we re-fit the data for a range of values of the SN brightness and computed the probability distribution from the resultant $`\mathrm{\Delta }\chi ^2`$. As can be seen from Table 4, the least constraining limit comes from the case where the afterglow evolves in a wind-stratified medium with the cooling break redward of the optical band, and even in this case, a SN brighter than $`\delta =1.40`$ mag is excluded at 99.7% confidence, and a SN as bright as SN 1998bw ($`\delta =0`$ mag) is ruled out at greater than 99.999% confidence.

The peak brightness and the time scales for SNe Ib/c are generally correlated such that fainter SNe may peak earlier (Iwamoto et al., 1998). It may be important to take this into account for our analysis, since the observations most sensitive to the presence of an underlying SN are all after the peak. To do this, we shifted the $`UBVRI`$ photometry of the (intrinsically-)fainter Type Ic SN 1994I (Richmond et al., 1996) to $`z=0.451`$, and derived the transformation between the redshifted SN 1998bw and 1994I light curves using a similar method as Bloom et al. (2002). This method is analogous to the “stretch” method for SNe Ia (Perlmutter et al., 1999). If we use this transformation in our model to transform the redshifted SN 1998bw light curve to the light curve of a SN fainter than SN 1998bw by $`\delta `$ magnitudes, then our least-constraining limit on an underlying SN becomes $`\delta =1.34`$ mag fainter than SN 1998bw (at 99.7% confidence). The agreement with the above limit indicates that the uncertainty in our knowledge of the the light-curve shape and luminosity scaling light-curve is not important for this analysis.

Leaving aside the SN issue, our fits provide a jet-break time of approximately 35 days. From the FREGATE 8 – 400 keV fluence of $`1.5\times 10^5`$ erg cm-2, we calculate the $`k`$-corrected isotropic-equivalent energy release (Bloom, Frail & Sari, 2001) in the $`\gamma `$-ray band, $`E_\gamma (1.3\pm 0.3)\times 10^{52}`$ erg. Applying the geometric correction from our measurement of the jet break (using the formulation and normalization of Frail et al. 2001), we obtain a jet opening angle of $`18^{}`$. Thus the true energy release is $`(6.5\pm 1.6)\times 10^{50}`$ erg — consistent with the clustering of energy releases around $`5\times 10^{50}`$ erg (Frail et al., 2001).

## 5 Conclusions

Here we report the search for an underlying SN in the afterglow of GRB 010921. Thanks to the superb photometric stability of HST and the $`N\times (N1)/2`$ subtraction technique, we have been able to trace the light curve of the afterglow of GRB 010921 over two months. The resulting photometry is unbiased by aperture effects that are so prevalent in simple aperture and PSF-fitting photometry. We report two results.

First, we find a jet break time of 35 days, using only optical data. Second, we find no evidence for an SN. A SN, if present, must be fainter than SN 1998bw by $`>1.34`$ mag at 99.7% confidence. To our knowledge, to date, this is the most stringent limit for an underlying SN associated with a cosmologically located GRB.

As noted in §1, the collapsar model as currently understood has little power in predicting the dispersion in the amount of 56Ni synthesized as compared to the energy in relativistic ejecta. Underlying SNe are directly powered by the former whereas the GRB is powered by the latter. Observations are needed to start mapping the distribution in these critical explosion parameters. Progress can be expected with such observational inputs accompanied by further refinements in the model. Motivated thus, we summarize in Table 5 the status of SN searches for all Table 5 all known GRBs with redshift222Beyond a redshift of $``$ 1.2, the distinctive and strong absorption blueward of 4000Å is redshifted out of the optical bands. The higher sensitivity of the optical bands thus favor searches for SNe below this redshift. less than 1.2.

The most secure case for an SN is that for GRB 011121 (Bloom et al., 2002; Garnavich et al., 2002). GRB 980326 shows a strong red excess at about a month but unfortunately a redshift is lacking. GRB 970228 shows a less clear excess but benefits from a known redshift. Stated conservatively, a SN as bright as that of SN 1998bw can be ruled out in GRB 000911. In all cases, save that of GRBs 980326 and 011121, the presence of a host with a magnitude comparable to the brightness of the peak of the SN, makes it difficult to identify an SN component. As noted in §3, “bumps” can arise from host contamination. Combining HST and ground based measurements (as is the case for GRB 970228) is prone to considerable errors (§3).

In summary, there is good evidence for an SN comparable in brightness to SN 1998bw in GRB 011121 (Bloom et al., 2002). For GRB 010921, using the HST observations reported here, we constrain any putative underlying SN to be 1.34 mag fainter than SN 1998bw. In the collapsar framework, this absence could be most readily attributed to the well known dispersion of the peak luminosity of Type Ib/c SNe.

An alternative possibility is that there may be more than one type of progenitor for long duration GRBs. Along these lines we note that Chevalier & Li (2000) claim that some afterglows (e.g. GRB 990123) are incompatible with a $`\rho r^2`$ inhomogeneous circumburst distribution whereas other afterglows (e.g. GRBs 970228 and 970508) are better explained by invoking an inhomogeneous circumburst medium. Progress requires both searches for underlying SNe as well as characterizing the circumburst medium via modeling of the early-time afterglow (e.g. GRB 011121, see Price et al. 2002a).

Finally, we note that the afterglow of GRB 010921 (and any coincident SN) was extincted by $`A_V^{\mathrm{MW}}0.5`$ mag of dust in the foreground, and $`A_V^{\mathrm{host}}1`$ mag of dust in the host galaxy (Table 3). Thus, in the future, using ACS aboard HST it should be possible to extend SN searches to at least 3 mag fainter than SN 1998bw, at which point it will be possible to detect more typical SNe Ib/c coincident with GRBs.

We thank Pete Challis for helpful discussions about WFPC2 reduction, and Megan Novicki and John Tonry for an advance copy of their $`N\times (N1)/2`$ subtraction paper. SRK and SGD thank NSF for supporting our ground-based GRB observing program. BPS and PAP thank the ARC for supporting Australian GRB research. Support for Proposal number HST-GO-08867.01-A was provided by NASA through a grant from Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA Contract NAS5-26555. KH is grateful for support under grant HST-GO-09180.07-A.


text:
# Study of Inclusive Production of Charmonium Mesons in 𝐵 Decays

## I Introduction

Studies of the inclusive production of charmonium mesons in $`B`$ decays provide insight into the physics of the underlying production mechanisms. Non-relativistic QCD (NRQCD) ref:bodwin , which may provide an explanation ref:psi2sbrat  for the unexpectedly large production of $`J/\psi `$ and $`\psi (2S)`$ mesons observed in $`p\overline{p}`$ collisions ref:abe , is an example of such a mechanism. NRQCD calculations use phenomenological matrix elements that should be applicable to a variety of production processes ref:leibovich , including such kinematically different regimes as hadron collisions at the Tevatron collider and $`e^+e^{}`$ collisions at the $`\mathrm{{\rm Y}}(4S)`$.

This paper presents an analysis of $`J/\psi `$, $`\chi _{c1}`$, $`\chi _{c2}`$, and $`\psi (2S)`$ mesons produced in $`B`$ decays at the $`\mathrm{{\rm Y}}(4S)`$ resonance. $`J/\psi `$ and $`\psi (2S)`$ mesons are reconstructed in the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ decay modes, and the $`\chi _{c1}`$ and $`\chi _{c2}`$ in the $`J/\psi \gamma `$ final states. We also reconstruct $`\psi (2S)J/\psi \pi ^+\pi ^{}`$.

The results include new measurements of previously observed decays ref:cleo ; ref:bellechic  as well as momentum and helicity distributions not previously measured.

## II The BaBar Detector and the PEP-II Collider

The BABAR detector is located at the PEP-II $`e^+e^{}`$ storage rings operating at the Stanford Linear Accelerator Center. At PEP-II, 9.0$`\mathrm{GeV}`$ electrons collide with 3.1$`\mathrm{GeV}`$ positrons to produce a center-of-mass energy of 10.58$`\mathrm{GeV}`$, the mass of the $`\mathrm{{\rm Y}}(4S)`$ resonance.

The BABAR detector is described elsewhere ref:babar ; here we give only a brief overview. Surrounding the interaction point is a 5-layer double-sided silicon vertex tracker (SVT), which provides precision spatial information for all charged particles, and also measures their energy loss ($`\mathrm{d}E/\mathrm{d}x`$). The SVT is the primary detection device for low momentum charged particles. Outside the SVT, a 40-layer drift chamber (DCH) provides measurements of the transverse momenta $`p_T`$ of charged particles with respect to the beam direction. The resolution of the $`p_T`$ measurement for tracks with momenta above 1 $`\mathrm{GeV}/c`$ is parameterized as

$$\frac{\sigma (p_T)}{p_T}=0.13p_T\%+0.45\%,$$
(1)
where $`p_T`$ is measured in $`\mathrm{GeV}/c`$. The drift chamber also measures $`\mathrm{d}E/\mathrm{d}x`$ with a resolution of 7.5%. Beyond the outer radius of the DCH is a detector of internally reflected Cherenkov radiation (DIRC), which is used primarily for charged hadron identification. The detector consists of quartz bars in which Cherenkov light is produced as relativistic charged particles traverse the material. The light is internally reflected along the length of the bar into a water-filled stand-off box mounted on the rear of the detector. The Cherenkov rings expand in the stand-off box and are measured with an array of photomultiplier tubes mounted on its outer surface. A CsI(Tl) crystal electromagnetic calorimeter (EMC) is used to detect photons and neutral hadrons, as well as to identify electrons. The energy resolution of the calorimeter is parameterized as:

$$\frac{\sigma (E)}{E}=\frac{2.3\%}{E^{\frac{1}{4}}}1.9\%,$$
(2)
where the energy $`E`$ is measured in $`\mathrm{GeV}`$. The EMC is surrounded by a superconducting solenoid that produces a 1.5-T magnetic field. The instrumented flux return (IFR) consists of multiple layers of resistive plate chambers (RPC) interleaved with the flux return iron. The IFR is used in the identification of muons and neutral hadrons.

Data acquisition is triggered with a two-level system. The first level uses fast algorithms implemented in hardware that examine tracks in the DCH and energetic clusters in the EMC. The second level retains events in which the track candidates point back to the beam interaction region, or in which the EMC cluster candidates are correlated in time with the rest of the event and exceed the energy of a minimum ionizing particle. Over 99.9% of $`B\overline{B}`$ events pass the second level trigger. A fraction of all events that pass the first level are passed through the second to allow monitoring of its performance.

## III Coordinate System and Reference Frames

We use a right-handed coordinate system with the $`z`$ axis along the electron beam direction and the $`y`$ axis upwards, with origin at the nominal beam interaction point. The polar angle $`\theta `$ is measured from the $`z`$ axis and the azimuthal angle $`\varphi `$ from the $`x`$ axis. Unless otherwise stated, kinematic quantities are calculated in the rest frame of the detector. The other reference frame we commonly use is the center of mass of the colliding electrons and positrons, which we call the center-of-mass frame.

## IV Event Selection

The data used in these analyses were collected between October 1999 and October 2000 and correspond to an integrated luminosity of 20.3$`\text{ fb}^1`$ taken on the $`\mathrm{{\rm Y}}(4S)`$ and 2.6$`\text{ fb}^1`$ taken off-resonance at an energy 0.04$`\mathrm{GeV}`$ lower than the peak, which is below the threshold for $`B\overline{B}`$ production and therefore includes only continuum processes.

We use an equivalent luminosity of simulated data ref:geant , including both $`B\overline{B}`$ and continuum, to study the efficiency of the analysis.

We require events to satisfy criteria that are intended to have high efficiency for $`\mathrm{{\rm Y}}(4S)`$ events while rejecting a significant fraction of continuum events and strongly suppressing beam gas events. The event must satisfy either the DCH components of both trigger levels or the EMC components of both, and have three or more high-quality tracks in the angular region with full tracking acceptance, $`0.41<\theta <2.54`$. Reconstructed charged particles are considered high-quality tracks if they have at least 12 hits in the DCH, $`p_T>100`$$`\mathrm{MeV}/c`$, and a point of closest approach to the beam spot of $`<3`$$`\mathrm{cm}`$ in $`z`$ and $`<1.5`$$`\mathrm{cm}`$ in $`xy`$.

The ratio of the second to the zeroth Fox-Wolfram moment $`R_2`$ ref:fox  measures how uniformly the energy in the event is distributed, distinguishing the more spherical $`B\overline{B}`$ events with small $`R_2`$ from the more jet-like continuum events at large $`R_2`$ (Fig. 1). We require $`R_2<0.5`$.

The primary event vertex is obtained from all charged particles with $`0.41<\theta <2.54`$. Particles contributing a large $`\chi ^2`$ are removed from the vertex and the process is iterated until stable. To reject events due to a beam particle striking the beam pipe or a residual gas molecule, we require the primary event vertex to be within 0.5$`\mathrm{cm}`$ of the beam spot in $`xy`$ and 6$`\mathrm{cm}`$ in $`z`$.

Finally, the event must include total visible energy (charged particles plus unassociated EMC clusters above 30$`\mathrm{MeV}`$) greater than 4.5$`\mathrm{GeV}`$.

The efficiency for simulated $`\mathrm{{\rm Y}}(4S)B\overline{B}`$ events is 95.4%. More relevant is the ratio $`C_E`$ of this efficiency to that for events containing the charmonium decay of interest. $`C_E`$ is calculated for each of the final states using simulated data. Inaccuracies in the simulation of the tracking systems produce an uncertainty of 1.1%, common to all modes.

### IV.1 Determination of Number of $`\mathrm{{\rm Y}}(4S)`$ Mesons

The number of $`\mathrm{{\rm Y}}(4S)`$ events satisfying the above selection criteria ($`N_\mathrm{{\rm Y}}`$) is obtained from the total number of events satisfying the criteria by subtracting the component due to the continuum.

$$N_\mathrm{{\rm Y}}=N_{on}\kappa R_{off}M_{on},$$
(3)
where

is the number of events satisfying the criteria in the on-resonance data set;

is the number of muon pairs in the on-resonance data set;

is the ratio of the number of events satisfying the hadronic selection criteria to the number of muon pairs in the off resonance data; and

$`=1.0000\pm 0.0025`$ allows for differences in the ratios of continuum and muon-pair cross sections and efficiencies between on-resonance and off-resonance data.

The two highest-momentum tracks in muon pair events must both deposit less than 1$`\mathrm{GeV}`$ in the EMC and satisfy $`|\mathrm{cos}\theta |<0.7485`$ in the center-of-mass frame (to be within the region of full tracking efficiency). The tracks must be within $`10^{}`$ of being back-to-back, and must have a combined mass greater than 7.5$`\mathrm{GeV}/c^2`$. Since the number of muon pairs in the on and off resonance data sets appears only as a ratio in Equation 3, $`N_\mathrm{{\rm Y}}`$ depends only weakly on the details of the selection criteria. Varying these over reasonable ranges changes $`N_\mathrm{{\rm Y}}`$ by 0.5%, which we take as a systematic error. The on-resonance data contains 7.8 times as many muon pairs as the off-resonance data.

Changes to the trigger configuration caused $`R_{off}`$ to vary from 4.89 to 4.94 during the period in which the data was collected. For the purposes of this calculation, the data is grouped into periods of compatible $`R_{off}`$. Combining the data in a single group changes $`N_\mathrm{{\rm Y}}`$ by 0.3%, which is taken as a systematic error. In principle, $`R_{off}`$ could also vary due to beam gas backgrounds. However, the $`z`$ distribution of the primary event vertex indicates that less than 0.1% of events are due to beam gas.

We find $`N_\mathrm{{\rm Y}}=(21.26\pm 0.17)\times 10^6`$. The 0.8% uncertainty is systematic and includes the 0.5%, 0.3% and 0.1% contributions from muon pairs, $`R_{off}`$, and beam gas described above. The largest component is from the 0.25% uncertainty on $`\kappa `$, which corresponds to a 0.6% uncertainty on $`N_\mathrm{{\rm Y}}`$.

## V $`J/\psi `$ Production

### V.1 $`J/\psi `$ Reconstruction

We reconstruct $`J/\psi `$ candidates in selected events using the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ final states. The leptons are required to satisfy the track quality criteria listed earlier. Electron candidates are further required to have a distance of closest approach to the beam line of less than 0.25$`\mathrm{cm}`$ to reject electrons produced by photon conversions.

Both leptons are required to fall in the angular range $`0.410<\theta <2.409`$$`\mathrm{rad}`$ (the overlap of the SVT and EMC coverage). Simulation indicates that $`(75.3\pm 0.9)`$% of $`J/\psi `$ decays give both leptons in this region. As described in Sec. IX.1, the $`J/\psi `$ momentum distribution in the simulation ref:jetset  is slightly different from the observed distribution. This difference, and a small variation of efficiency with momentum, produce the uncertainty in the efficiency.

We obtain the efficiency for the leptons to satisfy the quality criteria by comparing the performance of the independent SVT and DCH tracking systems in hadronic events. The corresponding uncertainty in the efficiency is 2.4% per $`J/\psi `$.

One particle in a $`J/\psi e^+e^{}`$ candidate must satisfy the “very tight” electron identification criteria described below. The other must satisfy the “tight” criteria displayed in square brackets.

* Difference between measured and expected energy loss in the DCH between $`2\sigma `$ and $`+4\sigma `$, where $`\sigma `$ is the measurement error \[$`3\sigma `$ and $`+7\sigma `$\];
* Ratio of energy measured in EMC to measured momentum $`E/p`$ in range 0.89 to 1.2 \[0.75 to 1.3\];
* Associated EMC cluster must include at least four crystals \[same\];
* Lateral energy distribution LAT ref:lat  of EMC cluster in range 0.1–0.6 \[0.0–0.6\];
* The $`A_{42}`$ Zernike moment ref:zern  of the EMC cluster $`<0.42`$ \[no requirement\]; and
* DIRC Cherenkov angle within $`3\sigma `$ of expected value \[no requirement\].

LAT is a measure of the radial energy profile of the cluster, and is used to suppress clusters from electronic noise (very low LAT) or hadronic interactions (high LAT). $`A_{42}`$ measures the azimuthal asymmetry of the cluster about its peak, distinguishing electromagnetic from hadronic showers.

We reduce the impact of bremsstrahlung by combining photons radiated by electron candidates with the track measured in the tracking system (“bremsstrahlung-recovery”). Such photons must have EMC energy greater than 30$`\mathrm{MeV}`$ and a polar angle within 35$`\mathrm{mrad}`$ of the electron direction. The azimuthal angle of the photon must be within 50$`\mathrm{mrad}`$ of the electron direction at the beamspot or be between this direction and azimuthal location of the electron shower in the EMC.

For $`J/\psi \mu ^+\mu ^{}`$ candidates, one muon candidate must satisfy “tight” criteria, while the other satisfies “loose” criteria (shown in square brackets):

* Energy in calorimeter between 0.05 and 0.40$`\mathrm{GeV}`$ \[$`<0.50`$$`\mathrm{GeV}`$ \];
* Number of IFR layers $`N_{IFR}2`$ \[same\];
* Particle penetrates at least 2.2 interaction lengths $`\lambda `$ of detector material \[$`2\lambda `$\];
* Measured penetration within $`\pm 0.8\lambda `$ of the value expected for that momentum \[$`\pm 1.0\lambda `$\];
* An average of less than 8 hit strips per IFR layer \[same\];
* The RMS of the hits per IFR layer less than 4 \[same\];
* For candidates in the forward endcap, the number of hit IFR layers divided by the total number of layers between the first and last hit layers must be $`>0.34`$ to reject beam background in the outermost layer \[$`>0.30`$\];
* $`\chi ^2`$ of the track fit in the IFR $`<3\times N_{IFR}`$ \[same\]; and
* $`\chi ^2`$ of match between track from SVT and DCH and that found in the IFR $`<5\times N_{IFR}`$ \[same\].

The particle identification efficiencies are obtained by comparing the yield of $`J/\psi `$ mesons applying various criteria to one or both tracks. The efficiency for $`J/\psi e^+e^{}`$ satisfying the angular acceptance and track quality criteria is 90.5% with a systematic error of 1.8%. For $`J/\psi \mu ^+\mu ^{}`$, it is 71.7% with a systematic error of 1.4%. The systematic errors are somewhat conservative, in that they include a component due to the $`J/\psi `$ statistics. Misidentification of hadrons as muons is higher than misidentification as electrons, producing $`J/\psi `$ background levels that are approximately a factor of two higher.

Finally, the $`J/\psi `$ candidate must have momentum in the center of mass $`p^{}<2.0`$$`\mathrm{GeV}/c`$. This requirement is fully efficient for a $`J/\psi `$ meson from $`B`$ decay but rejects approximately 74% of those produced in the continuum ref:cont .

More than one $`J/\psi `$ candidate may be found in an event. A second candidate is observed in 0.8% of events that include at least one.

### V.2 Extraction of Number of $`J/\psi `$ Mesons

The mass of the $`J/\psi `$ candidate is obtained after constraining the two leptons to a common vertex. In less than 1% of events, the vertex fit does not converge, and we instead use four-vector addition to obtain the candidate mass. The mass resolution is poorer by approximately 1% for these events. Figure 2 shows the mass distribution of the selected candidates in the two lepton modes.

The number of $`J/\psi `$ mesons in the mass window used in the fit (2.6–3.3$`\mathrm{GeV}/c^2`$ for the $`e^+e^{}`$ mode, 2.8–3.3$`\mathrm{GeV}/c^2`$ for $`\mu ^+\mu ^{}`$) is determined by a binned likelihood fit to the distribution. The background is represented by a third-order Chebychev polynomial. The probability distribution function (pdf) for the $`\mu ^+\mu ^{}`$ signal is the distribution from simulation convolved with a Gaussian distribution, with mean (allowing for a systematic shift between simulation and data) and width (allowing for poorer resolution) free to float in the fit. The additional smearing is required because the simulation underestimates the actual amount of material in the detector. The fit returns an offset of 3$`\mathrm{MeV}/c^2`$ and an additional resolution of 7.8$`\mathrm{MeV}/c^2`$. The total mass resolution in data is approximately 12$`\mathrm{MeV}/c^2`$. The simulation includes final state radiation ref:photos  and predicts that $`(2.7\pm 0.1)`$% of reconstructed $`J/\psi \mu ^+\mu ^{}`$ candidates fall outside the mass range used in the fit.

To include the impact of bremsstrahlung, the $`J/\psi e^+e^{}`$ signal pdf includes four components, corresponding to mesons where neither electron has undergone bremsstrahlung or at least one, and mesons for which the bremsstrahlung-recovery process has located a photon or not. The shapes are derived from simulated data, but the relative weights of three of the four components are allowed to float in the fit. The fraction of events that did not undergo bremsstrahlung but had a photon assigned by the bremsstrahlung-recovery process is nominally fixed to the value predicted by simulation, although it is varied to obtain a systematic error.

An estimated $`(7.3\pm 0.8)`$% of reconstructed $`J/\psi e^+e^{}`$ candidates fall outside the mass range, which can be compared to the value of 6.1% if the relative pdf weights are fixed to the values predicted by simulation.

We perform similar fits to the off-resonance data with all signal fit parameters fixed except for the number of mesons. The result is scaled by the ratio of on- to off-peak luminosity and subtracted to obtain the number of mesons attributable to $`B`$ decay, which appears in Table 1 as the net meson yield.

The fitting procedure is validated and systematic errors on its results are obtained by comparing the generated number of events with the fit number of events for many simulated mass distributions convolved with a Gaussian distribution. For the purposes of this test, we increase the statistics of the simulated data by relaxing the particle identification and track-quality requirements. We also test second and fourth-order Chebychev polynomials for the background pdf. We perform these tests for the $`\chi _{c1}`$, $`\chi _{c2}`$, and $`\psi (2S)`$ mass distributions as well.

We vary fit parameters that are fixed during the fit to obtain an additional systematic contribution. In the case of the $`J/\psi `$, we vary the bremsstrahlung-recovery error rate from one-half to twice its nominal value. The systematic errors on the fit yields are 0.7% for $`J/\psi \mu ^+\mu ^{}`$ and 4.1% for $`J/\psi e^+e^{}`$.

### V.3 Determination of the $`BJ/\psi X`$ Branching Fraction

We calculate values for the $`BJ/\psi X`$ branching fraction using the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ final states separately, then combine the two. The equation for the branching fraction is the same for both cases (and for the other mesons studied):

$$=\frac{N_\psi }{2N_\mathrm{{\rm Y}}ϵ_C_c}C_E,$$
(4)
where

is the net number of mesons in the mass fit range after continuum subtraction;

is the efficiency for a meson to satisfy the selection criteria, including the requirement that the mass fall in the mass fit range;

is world average ref:pdg2000  for the relevant secondary charmonium branching fraction. For the $`J/\psi `$, this is for $`J/\psi e^+e^{}`$ or $`J/\psi \mu ^+\mu ^{}`$;

corrects for the difference in event selection efficiency (Sec. IV) between generic $`B\overline{B}`$ events and charmonium events. It is equal to the efficiency for generic $`B\overline{B}`$ events divided by the efficiency for the relevant charmonium final state.

Table 1 summarizes the meson yields and efficiencies. It also presents the branching fraction product $`_c`$, an experimental quantity that does not depend on the secondary charmonium branching fractions. There is a 3.1% systematic error common to both modes—and, in fact, to all final states we study—due to acceptance (1.2%), track quality selection (2.4%), uncertainty on $`C_E`$ (1.1%), and number of $`\mathrm{{\rm Y}}(4S)`$ (0.8%).

The separate $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ branching fraction measurements are averaged to obtain the final result (Table 2). Each measurement is weighted in the average by the inverse of the square of the statistical error plus the square of the systematic errors unique to that mode. The common systematic error is the largest component of the $`(BJ/\psi X)`$ uncertainty.

## VI $`\chi _c`$ Production

We reconstruct $`\chi _{c1}`$ and $`\chi _{c2}`$ mesons in the $`J/\psi \gamma `$ decay mode. The $`J/\psi `$ mesons must satisfy the criteria listed in Sec. V.1, with the additional requirement that the $`J/\psi `$ candidate mass $`m`$ satisfy $`3.05<m<3.12`$$`\mathrm{GeV}/c^2`$ for $`e^+e^{}`$ decays and $`3.07<m<3.12`$$`\mathrm{GeV}/c^2`$ for $`\mu ^+\mu ^{}`$. An estimated $`(74.0\pm 0.4)`$% of $`J/\psi e^+e^{}`$ and $`(91.4\pm 0.3)`$% of $`J/\psi \mu ^+\mu ^{}`$ mesons that satisfy all other criteria satisfy this additional mass selection.

Photon candidates are EMC clusters in the angular range $`0.41<\theta <2.409`$$`\mathrm{rad}`$ with energy between 0.12$`\mathrm{GeV}`$ and 1.0$`\mathrm{GeV}`$. Hadronic showers are suppressed by requiring LAT less than 0.8 and $`A_{42}<0.15`$, while clusters from nearby hadronic showers are suppressed by requiring that candidates be at least 9 from all charged tracks.

Most photons satisfying these requirements are produced in $`\pi ^0`$ decay. We reject a candidate that, when combined with any other photon, produces a mass between 0.117$`\mathrm{GeV}/c^2`$ and 0.147$`\mathrm{GeV}/c^2`$. The second photon must have energy greater than 30$`\mathrm{MeV}`$ and LAT$`<0.8`$, with no requirement on $`A_{42}`$ or distance from charged tracks.

A systematic error due to the photon selection criteria is obtained by comparing the branching ratio $`\tau ^+h^+\pi ^0\pi ^0/\tau ^+h^+\pi ^0`$ in data to that in simulation, where $`h^+`$ is any charged track. We also vary the minimum energy requirement from 0.10 to 0.14$`\mathrm{GeV}`$ and test alternative $`\pi ^0`$ veto regions. We obtain a systematic error of 3.1% for the $`\chi _{c1}`$ and 4.4% for the $`\chi _{c2}`$ common to both the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ final states. An additional component to the systematic error arises from changes in the shape of the background, shown in Fig. 3, affecting the fit results. It is specific to each final state and mode, and amounts to 2.2% ($`e^+e^{}`$) and 1.3% ($`\mu ^+\mu ^{}`$) for the $`\chi _{c1}`$, and 9.2% ($`e^+e^{}`$) and 11.9% ($`\mu ^+\mu ^{}`$) for the $`\chi _{c2}`$.

The photon is constrained to originate at the $`J/\psi `$ vertex in the calculation of the $`\chi _c`$ four-momentum. We require $`p^{}<1.7`$$`\mathrm{GeV}/c`$, a requirement that is satisfied by $`\chi _{c1}`$ or $`\chi _{c2}`$ mesons from $`B`$ decays.

We determine the number of mesons from a fit to the plot of the mass difference between the candidate and the daughter $`J/\psi `$ masses (Fig. 3). We use different signal pdfs for the $`\chi _{c1}`$ and $`\chi _{c2}`$. These are formed by convolving the pdf calculated by simulation with a Gaussian distribution, where the offset and sigma are constrained to be the same for the $`\chi _{c1}`$ and $`\chi _{c2}`$. The background is described by a third-order Chebychev polynomial. Systematic errors on the fit are obtained as for the $`J/\psi `$. The correlation coefficient between the number of $`\chi _{c1}`$ and $`\chi _{c2}`$ mesons obtained from the fit is 0.19.

Equation 4 is used to determine inclusive $`B\chi _{c1}X`$ and $`B\chi _{c2}X`$ branching fractions separately for $`J/\psi e^+e^{}`$ and $`J/\psi \mu ^+\mu ^{}`$. In this case, $`_c=(\chi _{cJ}J/\psi \gamma )(J/\psi \mathrm{}^+\mathrm{}^{})`$, where $`\chi _{cJ}`$ is $`\chi _{c1}`$ or $`\chi _{c2}`$ and $`\mathrm{}`$ is $`e`$ or $`\mu `$. Table 1 summarizes the yields, efficiencies, uncertainties and branching fraction products $`(B\chi _{cJ}X)(\chi _{cJ}J/\psi \gamma )(J/\psi \mathrm{}^+\mathrm{}^{})`$. The 4.5% ($`\chi _{c1}`$) and 5.3% ($`\chi _{c2}`$) systematic errors common to both $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ include photon reconstruction in addition to the $`J/\psi `$ reconstruction items.

As with the $`J/\psi `$, the inclusive $`B`$ branching fractions to the $`\chi _{c1}`$ and $`\chi _{c2}`$ are calculated separately using the $`J/\psi e^+e^{}`$ and $`J/\psi \mu ^+\mu ^{}`$ decays. The two values are then combined, distinguishing uncertainties common to both from those unique to a single final state.

The branching fraction obtained for the $`\chi _{c2}`$ is comparable to that for the $`\chi _{c1}`$ and is summarized in Table 2. This result is consistent with a prediction from a color octet calculation ref:bodwinc2 , and is in contrast to the expectation of a null result in a factorization calculation ref:kuhn .

## VII $`\psi (2S)`$ Production

The reconstruction of the $`\psi (2S)`$ in the $`\mathrm{}^+\mathrm{}^{}`$ final state is very similar to the $`J/\psi `$ reconstruction outlined in Sec. V.1, with the $`p^{}`$ requirement tightened to $`p^{}<1.6`$$`\mathrm{GeV}/c`$. Figure 4 shows the resulting candidate mass distribution. A fit to extract the number of mesons in each plot is performed as for the $`J/\psi `$, but with the resolution and bremsstrahlung parameters fixed to the values found in the higher-statistics $`J/\psi `$ channels. These parameters are varied according to their uncertainties as one contribution to the systematic error on the fit; the remaining contributions are determined as for the $`J/\psi `$.

These data are used to calculate the branching fraction product $`(B\psi (2S)X)(\psi (2S)\mathrm{}^+\mathrm{}^{})`$, and are later used in the determination of the $`p^{}`$ distribution of $`\psi (2S)`$ mesons produced in $`B`$ decay. However, the extraction of the $`B\psi (2S)X`$ branching fraction requires the use of $`\psi (2S)e^+e^{}`$ and $`\psi (2S)\mu ^+\mu ^{}`$ branching fractions. Since this same data set has previously been used to measure these branching fractions ref:psi2s , we do not use $`\psi (2S)\mathrm{}^+\mathrm{}^{}`$ events to find $`(B\psi (2S)X)`$.

Instead, we use $`\psi (2S)J/\psi \pi ^+\pi ^{}`$ for this purpose. The reconstruction of a $`\psi (2S)`$ candidate in this final state starts with a $`J/\psi `$ candidate satisfying the tighter mass constraints used in $`\chi _{cJ}`$ reconstruction. All charged particles, including those failing the “good-track” criteria, are assumed to be pion candidates. The pion pair is required to be oppositely charged and to have a mass, calculated by four-vector addition, in the range 0.45 to 0.60$`\mathrm{GeV}/c^2`$. The mass distribution from simulation is compared to the measured ref:bes  distribution to obtain a systematic error of 0.5% on reconstruction efficiency. Finally, the $`p^{}`$ of the $`\psi (2S)`$ candidate is required to be less than 1.6$`\mathrm{GeV}/c`$.

Figure 5 displays the mass difference between the $`\psi (2S)`$ and the $`J/\psi `$ candidates separately for $`J/\psi e^+e^{}`$ and $`J/\psi \mu ^+\mu ^{}`$. As for the other final states, the distributions are fit to obtain the number of mesons. The resolution smearing parameters are not required to be the same for the two plots, but are consistent: $`1.5\pm 0.8`$$`\mathrm{MeV}/c^2`$ ($`e^+e^{}`$) and $`1.8\pm 0.5`$$`\mathrm{MeV}/c^2`$ ($`\mu ^+\mu ^{}`$). The secondary branching fractions in Eq. 4 are in this case $`_c=(\psi (2S)J/\psi \pi ^+\pi ^{})(J/\psi \mathrm{}^+\mathrm{}^{})`$.

## VIII Direct Branching Fractions

To obtain the branching fraction for $`J/\psi `$ mesons produced directly in the decay of $`B`$ mesons, we subtract the feeddown contributions to the inclusive branching fraction due to the decay of $`\chi _{c1}`$, $`\chi _{c2}`$, and $`\psi (2S)`$ mesons. For the $`\chi _{c1}`$ and $`\chi _{c2}`$, the feeddown branching fraction is $`(B\chi _{cJ}X)(\chi _{cJ}J/\psi \gamma )`$, while for the $`\psi (2S)`$, it is $`(B\psi (2S)X)(\psi (2S)J/\psi X)`$.

Similarly, the feeddown from the $`\psi (2S)`$ to the $`\chi _{c1}`$ and $`\chi _{c2}`$ is $`(B\psi (2S)X)(\psi (2S)\chi _{cJ}\gamma )`$.

Note that a number of uncertainties are common to both the inclusive and feeddown components, including track quality and particle identification criteria, $`(\chi _{c1}J/\psi \gamma )`$, and $`(\chi _{c2}J/\psi \gamma )`$. We use world average values ref:pdg2000  for the $`\psi (2S)`$ branching fractions. The resulting direct branching fractions are summarized in Table 2.

## IX $`p^{}`$ Distributions

The momentum distributions of charmonium mesons provide an insight into their production mechanisms. Since we do not fully reconstruct the $`B`$ meson, we cannot determine the meson momentum in the $`B`$ rest frame and instead use $`p^{}`$, the value in the $`\mathrm{{\rm Y}}(4S)`$ center-of-mass frame. The difference, due to the motion of the $`B`$ in the center-of-mass frame, has an rms spread of 0.12$`\mathrm{GeV}/c`$.

### IX.1 Inclusive $`p^{}`$ Distributions

To measure the $`p^{}`$ distributions of $`J/\psi `$, $`\chi _{c1}`$, $`\chi _{c2}`$, and $`\psi (2S)`$ mesons produced in $`B`$ decays, we create mass or mass-difference histograms of on-resonance candidates with $`p^{}`$ in the desired range. The $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ final states are again treated separately. The distributions are then fit, with all signal pdf parameters (other than the number of mesons) fixed to the values obtained from the earlier fits. The fits are performed for 100$`\mathrm{MeV}/c`$ wide $`p^{}`$ ranges, and in each case the sum of the yields differs from the original fit by fewer than ten events.

In the case of the $`J/\psi `$, we perform similar fits on the off-resonance data and perform a continuum subtraction for each $`p^{}`$ bin. Since there are no statistically significant off-resonance $`\chi _{c1}`$, $`\chi _{c2}`$, or $`\psi (2S)`$ signals, we do not perform a continuum subtraction in these cases.

The yield in each bin is corrected by the reconstruction efficiency obtained from simulated data, which decreases by approximately 10% between 0 and 2$`\mathrm{GeV}/c`$. The yield is then multiplied by an overall normalization factor for that particular final state and mode, which adjusts the sum of all bins to the earlier branching fraction measurement. We then perform a weighted average of the two distributions for the $`J/\psi `$, $`\chi _{c1}`$, or $`\chi _{c2}`$, or the four distributions for the $`\psi (2S)`$, to obtain the distributions shown in Fig. 68. For this purpose, we use the $`\psi (2S)e^+e^{}`$ and $`\psi (2S)\mu ^+\mu ^{}`$ branching fractions from Ref. ref:psi2s . In all cases, the distributions that are combined are consistent within statistical errors.

### IX.2 Direct $`p^{}`$ Distributions

The $`J/\psi `$ $`p^{}`$ distribution (Fig. 6) includes components due to mesons from the decays $`\chi _{c1}J/\psi \gamma `$, $`\chi _{c2}J/\psi \gamma `$, and $`\psi (2S)J/\psi X`$. To measure these distributions, we repeat the analysis with the data binned by the $`p^{}`$ of the $`J/\psi `$ daughter. The resulting $`J/\psi `$ feeddown distributions are presented in Fig. 9.

Note that we are using only the $`J/\psi \pi ^+\pi ^{}`$ decay mode to obtain the $`J/\psi `$ distribution from $`\psi (2S)`$ decay. In fact, 10.5% of $`\psi (2S)J/\psi X`$ decays are modes other than $`J/\psi \pi \pi `$. If we instead use the simulated $`J/\psi `$ distribution for this 10.5%, Fig. 9c changes by no more than a small fraction of the statistical error bar in any bin.

Subtracting these three components from the inclusive $`J/\psi `$ distribution in Fig. 6 leaves the contribution due to the $`J/\psi `$ mesons produced directly in $`B`$ decay (Fig. 10).

The superimposed histogram is a calculation of the expected distribution, which includes color octet and color singlet components. We use a recent NRQCD calculation ref:benekep  for the color octet component. The authors attribute the singlet component to $`J/\psi K^{()}`$ production, which we obtain from simulation. The two are normalized to obtain the best fit to our data. Possible sources of the apparent excess at low momentum are an intrinsic charm component of the $`B`$ ref:hou , the production, together with the $`J/\psi `$, of baryons ref:brodsky , or an $`s\overline{d}g`$ hybrid ref:eilam .

The small feeddown contribution to $`\chi _{c1}`$ and $`\chi _{c2}`$ from $`\psi (2S)`$ decay is calculated by simulation and is shown in Fig. 7.

## X $`J/\psi `$ Helicity

The helicity $`\theta _H`$ of a $`J/\psi \mathrm{}^+\mathrm{}^{}`$ candidate is the angle, measured in the $`J/\psi `$ rest frame, between the positively charged lepton and the flight direction of the $`J/\psi `$ in the $`\mathrm{{\rm Y}}(4S)`$ center-of-mass frame. A more natural definition would use the $`B`$ rest frame, but it cannot be determined in this analysis. Simulation indicates that the rms spread of the difference between the two definitions is 0.085 in $`\mathrm{cos}\theta _H`$.

### X.1 Inclusive Helicity Distribution

We proceed as for the $`J/\psi `$ $`p^{}`$ distribution, with data categorized into ranges of width 0.1 in $`\mathrm{cos}\theta _H`$ for two different momentum ranges, which we choose as $`p^{}<1.1`$$`\mathrm{GeV}/c`$ and $`1.1<p^{}<2.0`$$`\mathrm{GeV}/c`$. We fit the on- and off-resonance mass distributions to obtain yields in each bin and perform a continuum subtraction. We correct using the reconstruction efficiency obtained from simulation for that range, although we observe little dependence of efficiency on helicity. We then apply separate normalization factors to the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ data such that the total branching fraction (summed over the two $`p^{}`$ ranges) agrees with the value obtained earlier for that mode. The distributions from $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ are consistent and are averaged to obtain the helicity distributions for each of the two $`p^{}`$ ranges (Fig. 11).

We fit each distribution with a function $`1+\alpha \mathrm{cos}^2\theta _H`$ to obtain the polarization $`\alpha `$, where $`\alpha =0`$ indicates the sample is unpolarized, $`\alpha =1`$ transversely polarized, and $`\alpha =1`$ longitudinally polarized. The high $`p^{}`$ region, which includes the two-body $`B`$ decays, is more highly polarized, $`\alpha =0.592\pm 0.032`$, than the lower $`p^{}`$ region, $`\alpha =0.196\pm 0.044`$.

We assign a systematic error of 0.008 to these polarizations by instead considering the reconstruction efficiency to be independent of helicity.

### X.2 Direct $`J/\psi `$ Helicity

We determine the helicity distributions of $`J/\psi `$ mesons produced in the decay of $`\chi _{c1}`$, $`\chi _{c2}`$, and $`\psi (2S)`$ in the same way we calculate the $`p^{}`$ feeddown. Because of the limited statistics of these samples, we combine the two momentum regions used in the inclusive analysis. The resulting feeddown helicity distributions are shown together with the polarization fits in Fig. 12. We subtract these from the sum of the two distributions in Fig. 11 to obtain the helicity distribution for the $`J/\psi `$ produced directly in $`B`$ decay (Fig. 13). The polarization, $`\alpha =0.46\pm 0.06`$, is slightly out of the range $`0.33`$ to 0.05 predicted by an NRQCD calculation ref:fleming , but other authors have argued ref:ma  that relativistic corrections reduce the reliability of the calculation. The systematic uncertainty of 0.008 obtained above is small compared to the statistical error. This result is difficult to compare directly with that from CDF ref:cdfpol , due to the different mixture of $`b`$ mesons and baryons, and the distinction between the effective helicity calculated there and the true helicity.

## XI Summary

We have reported new measurements of $`B`$ meson decays to final states including charmonium mesons, which are summarized in Table 2. We have presented a number of momentum distributions. The distributions of the feeddown $`J/\psi `$ daughters of the $`\chi _{c1}`$, $`\chi _{c2}`$, and $`\psi (2S)`$ have not previously been measured, and allow us to more accurately determine the distribution for $`J/\psi `$ mesons produced directly in $`B`$ decay. The direct $`J/\psi `$ distribution is compared to a recent NRQCD calculation and appears to indicate an excess at low momentum.

The $`J/\psi `$ helicity distribution, which has also has not previously been published, indicates that the polarization of direct $`J/\psi `$ mesons is slightly out of the range predicted by an NRQCD calculation.

We are grateful for the excellent luminosity and machine conditions provided by our PEP-II colleagues, and for the substantial dedicated effort from the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and kind hospitality. This work is supported by DOE and NSF (USA), NSERC (Canada), IHEP (China), CEA and CNRS-IN2P3 (France), BMBF and DFG (Germany), INFN (Italy), NFR (Norway), MIST (Russia), and PPARC (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, Research Corporation, and Alexander von Humboldt Foundation.


text:
# Edge overload breakdown in evolving networks

## I Introduction

Large sparse networks are the underlying structure for transportation or communication systems, both man made (like computer networks AJB ; otherCN  or power grids PG ) or natural (like neural networks NN  or biochemical networks BioC ). These networks displays both randomness and some self-induced structure influencing the flow of transport and robustness against congestion or breakdown in the network. One of the most conspicuous structures among real-world communication networks is a highly skewed distribution of the degree (the number of neighbors of a vertex) AJB ; otherCN ; otherSF .

Avalanching breakdown in networks where the edges or vertices are sensitive to overload is a serious threat for real-world networks. A recent example being the black-out of 11 US states and two Canadian provinces the 10th August 1996 blackout . Recently the overload breakdown problem for vertices in growing networks with an emerging power-law degree distribution has been studied EGO2 . In the present paper we investigate the overload breakdown problem when edges (rather than vertices) are sensitive to overloading. We use the standard model for such networks—the Barabási-Albert (BA) model BA1 ; BA2 , but with a maximum load capacity assigned to each edge. The load is defined by means of the betweenness centrality—a centrality measure for communication and transport flow in a network BC . The procedure enables us to study overload breakdown triggered by the redistribution (and increase) of load in a growing network. This is in contrast to earlier models of cascading breakdown phenomena, all dealing with vertex breakdown, that has taken a fixed network as their starting point WATTS ; MGP .

## II Definitions

We represent networks as undirected and unweighted graphs $`G=(V,E)`$ where $`V`$ is the set of vertices, and $`E`$ is the set of unweighted edges (unordered pairs of vertices). Multiple edges between the same pair of vertices are not allowed.

### II.1 The Barabási-Albert model of scale-free networks

The standard model for evolving networks with an emerging power-law degree distribution is the Barabási-Albert model. In this model, starting from $`m_0`$ vertices and no edges, one vertex with $`m`$ edges is attached iteratively. The crucial ingredient is a biased selection of what vertex to attach to, the so called “preferential attachment:” In the process of adding edges, the probability $`P_u`$ for a new vertex $`v`$ to be attached to $`u`$ is given by BA3

$$P_u=\frac{k_u+1}{_{wV}(k_w+1)},$$
(1)
where $`k_u`$ is the degree of the vertex $`u`$. To understand the effect of preferential attachment, we will also investigate networks grown with a unbiased random attachment of vertices. Without the preferential attachment the networks are known to have an exponential tail of the degree distribution BA2 . The time $`t`$ is measured as the total number of added edges, which is different by factor $`m`$ from Refs. BA1 ; BA2  where $`t`$ is defined as the number of added vertices.

It should be noted that in very large communication networks, such as the Internet, the users can process information about only a subset of the whole network. How this affects the dynamics of network formation is investigated in Ref. mossa . In the present work we neglect such effects and assume linear preferential attachment.

### II.2 Load and capacity

To assess the load on the vertices of a communication network, or any network where contact between two vertices are established through a path in the network, a common choice is the betweenness centrality BC  which often is seen as a vertex quantity but has a natural extension to edges $`eE`$ GN :

$$C_B(e)=\underset{vV}{}\underset{wV\{v\}}{}\frac{\sigma _{vw}(e)}{\sigma _{vw}},$$
(2)
where $`\sigma _{vw}(e)`$ is the number of geodesics between $`v`$ and $`w`$ that contains $`e`$, and $`\sigma _{vw}`$ is the total number of geodesics between $`v`$ and $`w`$. $`C_B(e)`$ is thus the number of geodesics between pairs of vertices passing $`e`$; if more than one geodesics exists between $`v`$ and $`w`$ the fraction of vertices containing $`e`$ contributes to $`e`$’s betweenness.

In Ref. EGO2  (see also Ref. KAHNG ) the use of betweenness centrality as a load measure is given thorough motivations. These arguments are readily generalized to the case of edges sensitive to overloading: Suppose that $`\mathrm{\Lambda }`$ is the set of pairs of vertices with established communications through shortest paths at a given instant foot . Then let $`\lambda (e)`$ denote the load of $`eE`$ defined as the number of geodesics that contains $`e`$. Then we assume the effective load to be the average

$$\lambda (e)_\mathrm{\Omega }=\frac{1}{|\mathrm{\Omega }|}\underset{\mathrm{\Lambda }\mathrm{\Omega }}{}\lambda (e),$$
(3)
where $`\mathrm{\Omega }`$ is an ensemble of $`\mathrm{\Lambda }`$. To proceed, we restrict $`\mathrm{\Omega }`$ according to:

$$\mathrm{\Omega }=\{\mathrm{\Lambda }:|\mathrm{\Lambda }|=AN(N1)\},$$
(4)
where $`A`$ is constant with respect to $`N`$. This is to be interpreted that an element of $`\mathrm{\Omega }`$ is a set of $`AN(N1)`$ pairs of distinct vertices chosen uniformly at random, and thus corresponds to the case where the number of established communication routes ending at a specific vertex in average increases with $`N`$. This case can for example be expected in the early days of the Internet where the launches of new sites made the users browse a larger average number of sites. The case where the users at average connects to an $`N`$-independent number of others is discussed in Appendix A. The largest approximation, when using the betweenness as a load measure, is probably that routing protocols of e.g. the Internet has implicitly implemented load balancing foot ; OSPF ; huitema .

To introduce overloading to the dynamics we assign a capacity, or maximum value, $`\lambda ^{\mathrm{max}}(e)`$ to the load, the same for each edge, and say that the edge $`e`$ is overloaded if $`\lambda ^{\mathrm{max}}(e)<\lambda (e)_\mathrm{\Omega }`$. From the definition of $`\mathrm{\Omega }`$ we can see that our situation corresponds to having a maximum capacity on the betweenness centrality of the edges so that an edge is overloaded if $`C_B(e)>C_B^{\mathrm{max}}`$ (where $`C_B^{\mathrm{max}}`$ is constant). If an edge is overloaded it is simply removed from the graph, and the betweenness recalculated, if then another edge becomes overloaded it is removed, and so on. If more than one edge is overloaded at a time we choose the one to remove randomly. Multiple breakdowns during one time step defines a “breakdown avalanche.”

### II.3 Quantities for measuring network functionality

To measure the network functionality we consider three quantities—the number of edges $`L`$, inverse geodesic length $`l^1`$, and the size of the largest connected subgraph $`S`$: For the original BA model the number of edges increases linearly as $`L(t)=t`$ (i.e. one edge is added in unit time). But if an overload breakdown occurs in the system $`L`$ decreases, making it a suitable simplest-possible-measure of the network functionality. In a functional network a large portion of the vertices should have the possibility to connect to each other. In percolation and attack vulnerability studies of random networks one often uses $`S`$ to define the system as ‘percolated’ (or functioning), when the size of the largest connected subgraph $`S`$ scales as $`N`$ AJB ; ATTPER . One of the characteristic features of the BA model networks, as well as many real-world communication networks, is a less than algebraically increasing average geodesic length $`l`$. As the average geodesic length is infinite when the network is disconnected (as could be the case when an overload breakdown has occurred) we study the average inverse geodesic length NEWMAN :

$$l^1\frac{1}{d(v,w)}\frac{1}{N(N1)}\underset{vV}{}\underset{wV\{v\}}{}\frac{1}{d(v,w)},$$
(5)
which has a finite value even for the disconnected graph if one defines $`1/d(v,w)0`$ in the case that no path connects $`v`$ and $`w`$. To monitor the fragmentation of the network we will also measure the number of connected subgraphs $`n`$.

## III Simulation results

For relative small $`m`$, typical runs are exemplified in Fig. 1. For both random and preferential attachment $`S`$ reaches a critical time where after the network starts to break down, eventually $`S`$ reaches a steady state value. The breakdown develops differently in the two cases: For the random attachment the breakdown is relatively slow and the steady state value is high compared to the preferential attachment case where large successive avalanches fragments the network. The other two quantities reflects the same behavior: While the initial vertices gets joined into the network $`l^1`$ increases to an early maximum. After the decrease corresponding to the increase of $`l`$, $`l^1`$ decreases rapidly when the network becomes fragmented. $`L`$ shows the jagged shape, as expected correlated with that of $`S`$. As seen in Fig. 1(a) and (b), the discontinuity in $`L`$ (in the preferential attachment case), is less pronounced than that in $`S`$, so a small number of overloaded edges can be enough to cause large decrease in $`S`$. The reason for this behavior is that bridges (single edges interconnecting connected subgraphs) have a high betweenness and thus are prone to overloading. The number of connected subgraphs behaves qualitatively the same for random and preferential attachment. For other runs of the algorithm the breakdown can qualitatively be described as above. The averaged quantities varies relatively little, for example the peak-time for $`S`$ has a standard deviation of $`3`$%.

The corresponding overload case for vertices studied in Ref. EGO2  shows a similar time development with an period of incipient scale-freeness, an intermediate regime of breakdown and recovery (although the period of recovery is not as large for edges as for vertices), and a final breakdown to a large-$`t`$ state of disconnected clusters. One major difference between overload breakdown for vertices and edges is that the difference between random and preferential attachment is larger for edge overloading—edge robustness benefits more than vertex robustness from the geometry arising from random attachment.

Next we investigate the $`m`$-dependence. As seen in Fig. 2 the system becomes more and more robust when $`m`$ increases. This is of course expected since with a higher average degree more edges shares the load, so the maximal load can be expected to decrease. For high enough $`m`$ there are no avalanches, the largest connected component remains of the same size $`S=C_B^{\mathrm{max}}+1`$. When $`S=C_B^{\mathrm{max}}+1`$ the next edge attaching a new vertex will have $`C_B(e)=C_B^{\mathrm{max}}+1`$, and thus be overloaded. In most cases this will lead to removal of the newly added edge—otherwise another edge has to be overloaded at the same time, which is decreasingly likely with increasing $`m`$. In Fig. 2(b) we can see one exception to this interpretation at $`m=6`$: Here $`S`$ reaches $`C_B^{\mathrm{max}}+1`$ but starts to decay slowly at around $`t=8000`$. As mentioned, the largest connected subgraph is expected to become more stable as $`m`$ increases. If there is an $`m`$ above which $`S=C_B^{\mathrm{max}}+1`$ for arbitrary large $`t`$ above some $`t_0`$ is an open question. Comparing Fig. 2(a) and Fig. 2(b) shows that random attachment and preferential attachment have similar $`m`$-dependence behavior—the major difference being that preferential attachment has a much sharper increase of $`S`$; to be more precise the $`m`$ values that does not reach $`S=C_B^{\mathrm{max}}+1`$ for any $`t`$, have a lower value in the large-$`t`$ limit.

To get another angle of the mechanisms of the breakdowns for small $`m`$, we consider histograms of degree $`k_v`$ and betweenness $`C_B(e)`$. Figs. 3 and 4 shows these histograms both before and after the large drop in $`S`$ for $`m=m_0/2=4`$ and $`C_B^{\mathrm{max}}=500`$. (In the random attachment case this drop occurs at $`t_{\mathrm{drop}}1600`$, the corresponding value for preferential attachment is $`t_{\mathrm{drop}}2000`$.) For random attachment the difference between the histograms before and after the $`S`$-drop is distinctively smaller than for preferential attachment, just as expected from Fig. 2. The random attachment curves in Fig. 3(a) has a degree distribution of truncated exponential form both at the earlier and later times. In Fig. 3(a) it is exponential over two decades of $`P(k_v)`$, but falls off faster than exponentially for higher $`k_v`$. For preferential attachment the degree distributions (Fig. 3(b)) have a distinct difference—at $`t<t_{\mathrm{drop}}`$ there is an emergent power-law shape of the $`P(k_v)`$-curve, whereas at $`t>t_{\mathrm{drop}}`$ the shape is exponential, $`\mathrm{exp}(0.62k_v)`$, over five decades. To summarize, the degree distributions before and after the $`S`$-peak illustrates the same behavior as the time evolution of $`S`$—the breakdown in the preferential attachment case is both faster and more restructuring than in the random attachment case.

The betweenness distributions of Fig. 4 shows a peak that moves to higher $`C_B`$, as $`t`$ grows, until it reaches its maximal value at the time of the drop in $`S`$ and starts to decrease. For random attachment (Fig. 4 (a)) the shape of the distribution looks qualitatively the same before and after the drop, but for preferential attachment (Fig. 4 (b)) $`P(C_B)0`$ for betweenness smaller than the peak. The vertex betweenness distribution of the BA model is known to be strictly decreasing KAHNG , which would imply that the low-$`C_B`$ tails in Fig. 4 (b) (and most likely in Fig. 4 (a) as well) comes from a spread of the size of the largest cluster, rather than from a tail in the largest cluster’s betweenness distribution. Another feature of the betweenness histograms of Fig. 4 is the smaller peaks at low $`C_B`$ for $`t<t_{\mathrm{drop}}`$. These peaks corresponds to a sharp peak of the cluster size distribution just after the $`S`$-peak (see Fig. 5). Such smaller clusters have small average degree with many $`k_v=1`$ vertices, which all contributes to a peak at $`s`$ of the betweenness histograms. This explains the peak at $`C_B45`$ in the $`t=5000`$ curve of Fig. 4 (b).

The distribution of cluster sizes displayed in Fig. 5 gives some further insights: For $`t>t_{\mathrm{drop}}`$ of the random attachment curves shows a bimodal distribution as $`P(s)`$ is zero in the interval $`60s290`$. The preferential attachment curves, in contrast, has a long tail. Both the large $`s`$ peak for random attachment and the tail of preferential attachment corresponds to one single cluster. This is in striking contrast to the vertex overload case EGO2  where the network looses the unique largest component after the breakdown avalanches. As $`t`$ evolves well beyond $`t_{\mathrm{drop}}`$ the largest component peak decreases, and does thus not represent a giant component (a largest cluster proportional to $`N`$). The picture for both random and preferential attachment is thus that the system does not loose its unique largest cluster in a single breakdown avalanche—an avalanche rather results in a few isolated vertices or smaller clusters getting disconnected from the largest connected component.

The overall picture of the time evolution of $`S`$, $`L`$ and $`l^1`$ (Fig. 1), the $`m`$-dependence (Fig. 2), as well as the histograms of Figs. 3, 4 and 5 is that for small $`m`$, avalanching breakdowns fragments the network to a state from which it never recovers. For preferential attachment the newly fragmented network contains a single largest cluster with a well defined size, and the emergent scale-free degree distribution before $`t_{\mathrm{drop}}`$ is replaced by an exponential distribution. The breakdown for the random attachment case turns out to be less violent, and does not cause any major structural change. Furthermore, the difference between the random and preferential attachment cases is larger for edge breakdown than for the corresponding vertex breakdown model studied in Ref. EGO2 .

## IV Summary and conclusions

We have studied networks grown by the Barabási-Albert model for networks with emergent scale-freeness and edges sensitive to overloading. Except the preferential attachment defining the BA model, we also study an unbiased random attachment. We focus on the case where the number of established connections to random other vertices of the network scales linearly with the number of vertices in the network.

We find that for intermediate values of $`m`$ (the number of edges added per vertex) the network grows like the BA model up to a point where is starts to break down. After a number of avalanching breakdowns the network reaches a state characterized by many disconnected clusters from which a giant component never re-emerges (although, in the preferential attachment case, there will always be one single largest cluster much larger than any other). If the growth is by random attachment, the breakdown is less violent with smaller avalanches and no pronounced structural change. For large $`m`$ the steady state at large times is characterized by a constant largest cluster size.

In context of real world communication networks one can conclude that these would benefit from being grown by random rather than preferential attachment (and this difference being larger for edge overload than for vertex overload studied in Ref. EGO2 ). In the vertex overload case avalanches proceeds until the network is fragmented into small clusters; in the edge overload problem there is still one large component after the breakdowns, thus we infer that for real-world communication networks, vertex overloading is a greater threat than edge overloading, and congestion control in telecommunication networks kihl  and Internet routing protocols huitema  should focus on balancing the vertex rather than edge load. Only if the capacity of vertices (servers etc.) grows significantly faster than the capacity for edges (cables etc.), edge overload breakdown is a potential threat for avalanching breakdowns that is triggered by the change of load in a growing network.

## Acknowledgements

The author thanks Beom Jun Kim for discussions. This work was partially supported by the Swedish Natural Research Council through Contract No. F 5102-659/2001.

## Appendix A Intrinsic communication activity

This paper deals mainly with the case where the average user of a growing communication network communicates with a number of others that increases linearly with $`N`$. One can also imagine a case where, even though the network grows, the users in average communicates with a network size independent number of others; which is the topic of the present Appendix. (In Ref. EGO2  this scenario was termed “intrinsic communication activity.”) The behavior of real communication networks lies, presumably, between these two extremes.

### A.1 Definitions

To implement the situation of intrinsic communication activity, we modify Eq. 4 to:

$$\mathrm{\Omega }^{}=\{\mathrm{\Lambda }:|\mathrm{\Lambda }|=A^{}N\}$$
(6)
where $`A^{}`$ is constant with respect to $`N`$. I.e. the users has the $`N`$-independent average number $`A^{}`$ of established contacts through shortest routes. Averaging the load over $`\mathrm{\Omega }^{}`$ according to Eq. 3 gives:

$`\lambda (e)_\mathrm{\Omega }^{}`$ $`=`$ $`{\displaystyle \frac{1}{|\mathrm{\Omega }^{}|}}{\displaystyle \underset{\mathrm{\Lambda }\mathrm{\Omega }^{}}{}}{\displaystyle \underset{(w,w^{})\mathrm{\Lambda }}{}}{\displaystyle \frac{\sigma _{ww^{}}(e)}{\sigma _{ww^{}}}}`$ (7)
$`=`$ $`{\displaystyle \frac{A^{}}{N}}C_B(e).`$
From this we see that having a constant capacity for the load $`\lambda (e)`$ corresponds to having a limit on $`C_B(e)`$ that increases with $`N`$. Thus we view $`e`$ as overloaded if $`C_B(e)`$ exceeds $`C_B^{\mathrm{max}}=Nc_B^{\mathrm{max}}`$ (where $`c_B^{\mathrm{max}}`$ is constant).

### A.2 Results

In the vertex overload breakdown problem, the case of intrinsic communication activity has a more complex dynamics than the extrinsic communication activity case (studied in the main part of the text), with giant components forming only occasionally for some sets of parameter values EGO2 . For edge overload breakdown, on the other hand, the dynamics of a system with intrinsic communication activity seems very simple with no avalanching breakdowns and no qualitative difference between preferential and random attachment, see Fig. 6. We can also notice that the measured quantities has a power-law dependence of $`t`$. (Fig. 6 is constructed from one run with random and preferential attachment respectively.) For large times ($`1000t5000`$) the exponent $`\alpha `$ for the time development of the respective quantity is (in the large $`t`$ limit): $`\alpha _{l^1}0.6`$, $`\alpha _L=1.0`$, and $`\alpha _S=\alpha _n=0.50`$ for both (a) and (b). Initially $`\alpha _{l^1}`$ is closer to zero, for $`100t1000`$ we have $`\alpha _{l^1}0.5`$. To illustrate the consistency of the exponents we note that

$$\underset{eE}{}C_B(e)=\underset{vV}{}\underset{wV\{v\}}{}d(v,w)=n\frac{N}{n}\left(\frac{N}{n}1\right)l_{\mathrm{CS}},$$
(8)
where $`l_{\mathrm{CS}}`$ is the average geodesic length for a connected subgraph, and $`d(v,w)=0`$ if $`v`$ and $`w`$ are disconnected. This yields

$$C_B(e)m\left(\frac{N}{n}1\right)l_{\mathrm{CS}}\underset{vV}{\mathrm{max}}C_B(e).$$
(9)
If one assumes that $`C_B(e)\mathrm{max}_{vV}C_B(e)`$ and $`Nn`$ we have $`l_{\mathrm{CS}}n`$. Making the crude approximation $`l^1l_{\mathrm{CS}}^1`$ gives $`\alpha _{l^1}\alpha _n`$, which holds well for small $`t`$. As $`t`$ increases the spread in shape of the connected subgraphs becomes larger so the $`l^1l_{\mathrm{CS}}^1`$ approximation becomes worse which is seen as a slight increase of the slope $`\alpha _{l^1}`$. That the approximation $`\alpha _{l^1}\alpha _n`$ is rather good throughout the range of $`t`$ is also reflected in that the average size of connected components $`N/n`$ is never very far from $`S`$: At $`t=5000`$ we have (see Fig. 6(a)) $`N/n50`$ and $`S=57`$ for random attachment, and $`N/n53`$ and $`S=56`$ for preferential attachment. In this approximation we see that $`lnN^{1/2}`$ so the small average geodesic length is lost within the connected subgraphs. If $`c_B^{\mathrm{max}}`$ is chosen larger, so the network initially grows without edges breaking, there are no large avalanches but a crossover to the behavior seen in Fig. 6.