text:
# THE PROPAGATION AND ERUPTION OF RELATIVISTIC JETS FROM THE STELLAR PROGENITORS OF GRBs

## 1 INTRODUCTION

It is now generally acknowledged that “long-soft” gamma-ray bursts (henceforth, just “GRBs”) are a phenomenon associated with the deaths of massive stars. In addition to the observed association with star-forming regions in galaxies (Vreeswijk et al., 2001; Bloom, Kulkarni, & Djorgovski, 2002; Gorosabel et al., 2003); “bumps” observed in the afterglows of many GRBs (Reichert, 1999; Galama et al., 2000; Bloom et al., 2002; Garnavich et al., 2003; Bloom, 2003); spectral features like a WR-star in the afterglow of GRB 021004 (Mirabal et al., 2002); and the association of GRB 980425 with SN 1998bw (Galama et al., 1998), there is now incontrovertible evidence that at least one GRB (030329) was accompanied by a bright, energetic supernova of Type Ic, SN 2003dh (Price et al., 2003; Hjorth et al., 2003; Stanek et al., 2003). Thus some, if not all GRBs, are produced when the iron core of a massive star collapses either to a black hole (Woosley, 1993; MacFadyen & Woosley, 1999) or a very rapidly rotating highly magnetic neutron star (Wheeler et al., 2000), producing a relativistic jet. Other alternatives in which the GRB occurs after the supernova (Vietri & Stella, 1998) cannot simultaneously explain GRB 030329 and SN 2003dh.

At the same time, the general class of high energy transients once generically called “gamma-ray bursts” has been diversifying. For some time, soft gamma-ray repeaters (SGRs) and “short hard gamma-ray bursts” have enjoyed a distinct status and, presumably, a separate origin. In addition, we now have cosmological “X-ray flashes” (“XRFs”; Heise et al, 2001; Kippen et al., 2003), “long, faint gamma-ray bursts” (in’t Zand et al., 2003), and lower energy events like GRB980425 (Galama et al., 1998). Is a different model required for each new phenomenon, or is some unified model at work, as in active galactic nuclei, a model whose observable properties vary with its environment, the angle at which it is viewed, and perhaps its redshift?

The answer is probably “both”. Even within the confines of the collapsar model, there are Types I, II and III (Heger et al., 2003). Type I happens only in massive stars that make their black holes promptly (MacFadyen & Woosley, 1999). This is what people generally have in mind when they use the word “collapsar”. In Type II though, a similar mass black hole forms by fall back and the burst lasts much longer (MacFadyen, Woosley, & Heger, 2001). Type III, happens only at very low metallicity and requires Pop III stars as its progenitors (Fryer, Woosley, & Heger, 2001). The event is very energetic, but highly redshifted. These other models may be particularly appropriate for unusually long bursts. A jet that wavered in its orientation in any of these events would emerge with a greater load of baryons. So would a jet whose power declined in a time less than that required to traverse the star. A jet in a blue supergiant would become choked because the power of the central engine has decreased substantially when the jet is still deep inside the star. Thus any sort of collapsar happening in a blue supergiant instead of a Wolf-Rayet star would produce a different sort of transient. Finally, it may just be that jets are just born with different baryon loadings, perhaps relating to the mass of the presupernova star and its distribution of angular momentum and magnetic field. Our present understanding of nature allows for all these solutions.

However, it also seems to reason that not all jets are the same, even for Type I collapsars, and that even if they were, different phenomena would be seen at different angles. We are thus motivated to consider the observational consequences of highly relativistic jets as they propagate through, and emerge from massive stars. What would they look like if seen from different angles? In particular, what is the distribution with polar angle of the energy and Lorentz factor? Might softer, less energetic phenomena be observed off axis and with what frequency?

Jets inside massive stars have been studied numerically in both Newtonian (MacFadyen & Woosley, 1999; MacFadyen, Woosley, & Heger, 2001; Khokhlov et al., 1999) and relativistic simulations (Aloy et al., 2000; Zhang, Woosley, & MacFadyen, 2003, henceforth Paper 1) and it has been shown that the collapsar model is able to explain many of the observed characteristics of GRBs. These previous studies have also raised issues which require further examination, especially with higher resolution. For instance, the emergence of the jet and its interaction with the material at the stellar surface and the stellar wind could lead to some sort of “precursor” activity. The cocoon of the jet will also have different properties from the jet itself and shocks within the cocoon or with external matter could lead to $`\gamma `$-ray and hard X-ray transients (Ramirez-Ruiz, Celotti, & Rees, 2002). There is also the overarching question of whether jets, calculated in two dimensions with assumed axial symmetry, are stable when studied in three dimensions.

Here we examine in two- and three-dimensional numerical studies, the interaction of relativistic jets with the outer layers of the Wolf-Rayet stars thought responsible for GRBs. We follow the emergence of the jets for a sufficient length of time to ascertain the properties of the cocoon explosion that surrounds the GRB and we study the dependence of the results on the dimensionality of the calculation. Such explosions are visible to a much greater angle and to a much larger number of observers. The observed phenomenon could be a hard X-ray flash or a weak GRB. We also find that the properties of the emergent jet are quite sensitive to whether the jet maintains its orientation (to within a critical angle) for durations of $`10\mathrm{s}`$ or so, the time it takes the jet to traverse the star. Jets that waiver may make hard transients of some sort, but not GRBs.

## 2 Stellar Model and Numerical Methods

### 2.1 Progenitor Star

A collapsar is formed when the iron core of a rotating massive star collapses to a black hole and an accretion disk. The interaction of this disk with the hole, through processes that are still poorly understood, produces jets with a high energy to mass ratio. Here we are concerned primarily with the propagation of these relativistic jets, their interactions with the star and the stellar wind, and the observational implications, and not so much with how they are born. In fact, our results here would be the same if the jet were produced by some other means (Wheeler et al., 2000) inside the same star.

Unlike in our previous study (Paper 1), the initial stellar model is taken directly from a presupernova star with a very finely zoned surface (Table 1; Fig. 1). Because we are only interested in events that happen outside about $`10^{10}\mathrm{cm}`$ on a time scale of, at most, tens of seconds, modification of the initial star by core collapse and rotation is negligible. Specifically, our initial model is a $`15\mathrm{M}_{}`$ helium star calculated by Heger & Woosley (2003). This helium star has been evolved to iron-core collapse while following the transport of angular momentum and including the effects of rotational mixing as in Heger, Langer, & Woosley (2000). The initial star was a nearly pure helium star rotating rigidly with a surface velocity equivalent to 10% Keplerian. This is a typical value in the study of Heger, Langer, & Woosley (2000). Presumably the star lost its envelope to a companion early in helium burning. Further mass loss and the transport of angular momentum by magnetic fields were ignored so as to form a big iron core that would very likely collapse to a black hole and give ample angular momentum for a disk. The radius of the helium star at onset of collapse is $`8.8\times 10^{10}\mathrm{cm}`$ and the surface of the star is very finely zoned (surface zoning $`10^{21}\mathrm{g}`$). Despite the fact that angular momentum transport was followed in the initial model to show its suitability as a collapsar (Fig. 1), rotation plays no role in the present study.

Outside the star, the background density, which comes from the stellar wind, is assumed to scale as $`R^2`$. The actual value of this density, unless it is very high, plays no role in determining the final answer at the small scales considered here ($`<2\times 10^{12}`$ cm), but finite value is needed in order to stabilize the code. We assumed a value at $`R=10^{11}\mathrm{cm}`$ of $`5\times 10^{11}\mathrm{g}\mathrm{cm}^3`$, where $`R`$ is the distance from the center of the star. This corresponds to a mass loss rate of $`1\times 10^5\mathrm{M}_{}\mathrm{yr}^1`$ for a wind velocity of $`1000\mathrm{km}\mathrm{s}^1`$ at $`10^{11}\mathrm{cm}`$. This mass loss rate is rather low for typical Wolf-Rayet stars of this mass in our galaxy (Nugis & Lamers, 2000), but might be more appropriate for the low metallicity in the GRB neighborhood (Crowther et al., 2002). A small value was also taken in order to be consistent with the assumed progenitor structure that was calculated with zero mass loss.

### 2.2 Computer Code

We employed a multi-dimensional relativistic hydrodynamics code that has been used previously to study relativistic jets in the collapsar environment (Paper 1). Briefly, our code employs an explicit Eulerian Godunov-type shock-capturing method (Aloy et al., 1999). The governing equations of relativistic hydrodynamics with a causal equation of state can be written as a system of conservation laws for rest mass, momentum, and energy. To solve the equations numerically, each spatial dimension is discretized into cells. Using the method of lines, the multi-dimensional time-dependent equations can then be evolved by solving the fluxes at each cell interface. In order to achieve high order accuracy in time, the time integration is done using a high order Runge-Kutta scheme (Shu & Osher, 1988), which involves prediction and correction. An approximate Riemann solver (Aloy et al., 1999) using the Marquina’s algorithm is used to compute the numerical fluxes from physical variables: pressure, rest mass density, and velocity at the cell interface. The values of the physical fluid variables at the cell interface are interpolated using reconstruction schemes (e.g., piecewise parabolic method, Colella & Woodward, 1984). This reconstruction procedure ensures high order accuracy in spatial dimensions. The conserved variables: rest mass, momentum and energy are evolved directly by the scheme. In each time step, physical variables, such as pressure, rest mass density and velocity, which are necessary in calculating the numerical fluxes, can be recovered from conserved variables by Newton-Raphson iteration. The code can operate in Cartesian, cylindrical, or spherical coordinates. Approximate Newtonian gravity is implemented by including source terms in the equations. Total energy, which includes kinetic and internal energy, in the laboratory frame are not exactly conserved due to gravitational potential energy and the way of implementing gravity as source terms. However, the code conserves energy to machine precision when gravity is turned off. Furthermore, the gravitational potential energy is negligible for relativistic material in our calculations. For instance, assuming a $`10\mathrm{M}_{}`$ central point mass, a fluid element with a velocity of $`0.5c`$ at the inner boundary of our computational grids, $`10^{10}\mathrm{cm}`$ (§ 3.1, § 4.1), has potential energy of less than 0.1% of its kinetic energy. Thus we expect that our calculations conserve energy to an accuracy of better than 0.1%. In simulations present in this paper, A gamma-law equation of state with $`\gamma =4/3`$ is used for the simulations presented in this paper.

## 3 Two-Dimensional Models

### 3.1 Model Set Up and Definitions

The mass interior to $`1.0\times 10^{10}\mathrm{cm}`$ is removed from the presupernova star and replaced by a point mass. No self gravity is included. This should suffice since we are studying phenomena that happen on a relativistic time scale and the speed of sound is very sub-luminal. While jets presumably go out both axes, we follow here only one, assuming symmetry in the other hemisphere. Jets are injected along the rotation axis (the center of the cylindrical axis of the grid) through the inner boundary. Each jet is defined by its power (excluding rest mass energy), $`\dot{E}`$, its initial Lorentz factor, $`\mathrm{\Gamma }_0`$, and the ratio of its total energy (excluding rest mass energy) to its kinetic energy, $`f_0`$. Our previous studies in Paper 1 have shown that a jet that starting with a half-opening angle of $`20\mathrm{°}`$ and a high Lorentz factor, $`\mathrm{\Gamma }50`$ at $`2000\mathrm{km}`$, will be shocked deep inside of the star. By the time it reaches $`10^{10}\mathrm{cm}`$, a jet should have a large ratio of internal energy to rest mass, a half-opening angle of about $`5\mathrm{°}`$, and a Lorentz factor, $`\mathrm{\Gamma }510`$.

Though GRBs observationally have highly variable properties, the jet power was taken here to be constant for the first $`20\mathrm{s}`$, then turned down linearly during the next 10 seconds. That is, during the interval 20 to 30 s, the power scaled as $`(30\mathrm{s}t)/10\mathrm{s}`$ decreasing to zero at $`30\mathrm{s}`$. During the declining phase, the pressure and density remained constant, and the Lorentz factor was calculated from the internal energy, density, and power. After 30 seconds, a pure outflow boundary condition was used for the inner boundary. The axis of the jet is defined as the $`z`$-axis in all cases and the jets were initiated parallel to that axis in a region that subtended a half-angle of 5 degrees as viewed from the origin.

Four models were calculated which span a range of energies and Lorentz factor (Table 2). In Model 2A, a total energy of $`5\times 10^{51}\mathrm{erg}`$ is injected, comparable to the results of Frail et al. (2001), but larger than the results of Panaitescu & Kumar (2001) by an order of magnitude (Table 3). In Models 2B and 2C, higher and lower energy deposition rates were employed. For the parameters chosen, a jet with an initial Lorentz factor of 5 or 10 at $`10^{10}\mathrm{cm}`$ should have a final Lorentz factor of $`180200`$ if all internal energy is converted into kinetic energy. The zoning employed in the models is given in Table 4. Typically over three million zones were used.

Model 2C used an extended grid in the $`r`$-direction, which allowed the greater lateral expansion of the jet to be followed at late times. This was necessary because of its lower energy and greater expansion during the time it took to traverse the extent of the $`z`$-grid. Model 2T used conditions like those of Model 2B, but a grid that was both smaller and coarser, chosen to be equivalent to that used in the three-dimensional calculations described later in the paper (§ 4).

### 3.2 Results in Two Dimensions

The relativistic jet begins to propagate along the $`z`$-axis shortly after its initiation. In agreement with previous studies (Aloy et al., 2000, and Paper 1), the jet consists of a supersonic beam, a shocked cocoon, a bow shock, and is narrowly collimated. Some snapshots of Models 2A, 2B, and 2C are given in Figs. 2, 3, and 4. The resulting “equivalent isotropic energies” for matter exceeding a certain Lorentz factor are also given, at $`70\mathrm{s}`$, for the same three models in Figs. 5, 6, and 7. The latter figures also give the estimated terminal Lorentz factor assuming that internal energy along the radial line of sight is converted into kinetic energy. The equivalent isotropic energy is defined as that energy an isotropic explosion would need in order to give the calculated flux of energy along the line of sight. It is a way of mapping the highly asymmetric two-dimensional results into equivalent one-dimensional models. The change in Lorentz factor between $`70\mathrm{s}`$ and infinity is not very great, except in situations where matter that was highly relativistic to begin with receives an additional boost in its frame by expansion. The duration of the calculation was set by how long it took the jet to reach the end of the simulation grid and would be costly to increase, but as the figures show, especially in Model 2C, there was still an interesting amount of internal energy at $`70\mathrm{s}`$.

The fractions of energies inside a certain angle to the total energies on the grid are given in Fig. 8. In all cases studied the high Lorentz factor characteristic of common GRBs is confined to a narrow angle of about 3 to 5 degrees with a maximum equivalent isotropic energy in highly relativistic matter along the axis of $`3\times 10^{53}`$ to $`3\times 10^{54}\mathrm{erg}`$. At larger angles there is significant energy, though, and Lorentz factor, $`\mathrm{\Gamma }`$ 10 to 20. At an angle of 10 degrees for example, the equivalent isotropic energy in matter with $`\mathrm{\Gamma }>20`$ is $`10^{52}\mathrm{erg}`$ in Model 2A and even $`10^{51}\mathrm{erg}`$ in Model 2C. At larger angles there is less energy, but still the possibility of low power transients of hard radiation. Model 2B did not eject much material with $`\mathrm{\Gamma }10`$ at angles $`10`$ degrees. The equivalent isotropic energy at larger angles ($`>2\mathrm{°}`$) for all three models can be fit well by a simple power-law, $`E_0\times (\theta /2\mathrm{°})^3\mathrm{erg}`$, where $`E_0`$ is $`1.5\times 10^{54}`$, $`4.5\times 10^{54}`$ and $`6.8\times 10^{53}`$, for models 2A, 2B and 2C, respectively. The ratios of the values of $`E_0`$, $`1.5:4.5:0.68`$, are very close to those of the energy deposition rates, $`1.0:3.0:0.5`$. Inside $`2\mathrm{°}`$, the distributions of energy and Lorentz factor are roughly flat, $`E_0`$. The values of $`E_0`$ can be roughly estimated from the total injected energy and the percentage of the energy in the jet core. For example, about $`40\%`$ of the total energy on the grid is contained inside 2 degrees for Model 2A (Fig. 8). About $`5\times 10^{51}\mathrm{erg}`$ (Table 2) was injected for both jets in Model 2A. One can estimate that the equivalent isotropic energy inside 2 degrees is $`1.6\times 10^{54}\mathrm{erg}`$ for Model 2A. More models need to be calculated to find how the distributions of energy and Lorentz factor at breakout depend on initial conditions.

Table 3 gives the energies for various components of the ejecta at a time $`40\mathrm{s}`$ after the jet was initiated. There still remains considerable internal energy available for conversion to kinetic at this time (which was chosen as a time when most of the relativistic matter was still on the grid), but this will affect mostly the matter with high Lorentz factor, $`\mathrm{\Gamma }>10`$. The results show that most of the energy injected as jets at $`10^{10}\mathrm{cm}`$ emerges in matter that is still relativistic. Only a small fraction goes into sub-relativistic expansion, nominally $`v<0.5c`$ and consequently, only a little of the jet energy is used to blow up the star. If there were no other energy sources the kinetic energy of the supernova would be $`2\times 10^{51}`$ erg.

The total energy of relativistic ejecta is also useful for comparison to radio observations of the afterglows of GRBs. For example, Li & Chevalier (1999) have placed limits on the total energy of ejecta with $`v>0.5c`$ in SN 1998bw. The limit, $`3\times 10^{50}\mathrm{erg}`$, may be an approximate estimate of the actual energy. The corresponding energy for Model 2C here is four times larger (Table 3). It seems likely that lower energy models than Model 2C could be constructed that would still give high Lorentz factors and equivalent isotropic energies in a narrow range of angles around the polar axis. That is, GRB 980425 may have been a harder, more energetic GRB seen off axis. However, the low energy of that burst, some five orders of magnitude less that that expected for a centrally observed GRB (Fig. 4) remains surprising. Either GRB 980425 was observed at a polar angle larger than 15 degrees, or the burst was weaker at all angles because of baryon loading (§ 1; Woosley, Eastman, & Schmidt 1999).

The possibility that the off-axis emission from material with $`\mathrm{\Gamma }`$ 10 to 20 corresponds to XRFs is discussed in § 5.

In 2D simulations with cylindrical coordinates, there is an imposed symmetry axis of the coordinate system. It is very important to repeat these calculations in 3D to ensure that our 2D results are valid and to examine three-dimensional jet instabilities. Three-dimensional calculations are, however, computationally expensive. So we have to use lower resolution for 3D calculations. In order to compare 2D and 3D results with the same resolution and identical initial parameters, we did a lower resolution 2D run, Model 2T (Tables 2 and 4). A comparison of two- and three-dimensional results will be discussed in § 4. It is also interesting to compare Model 2B and 2T, which had identical parameters but in different resolution. Qualitatively the results are similar. In both models, the jet emerges from the star with a cocoon surrounding the jet beam and a dense “plug” at the head of the jet (Fig. 9). There are, however, noticeable differences, as we expected. In particular, it takes less time for the jet in Model 2T to emerge from the star than that in Model 2B, because they have different numerical viscosity. The larger numerical viscosity in Model 2T is due to its lower resolution. Despite morphologic differences between Model 2B and 2T (Fig. 9), the energetics of these two models are quite similar (Table 5, Fig. 10). The distributions of equivalent isotropic energy versus angle for the jet core ($`3`$ degrees) in Models 2B and 2T are very similar. However, Model 2B has more mildly relativistic material ($`2<\mathrm{\Gamma }<10`$) in its cocoon than Model 2T.

## 4 Three-Dimensional Models

For the three-dimensional models, the same helium star was remapped into a three-dimensional Cartesian grid. Because of the greater cost of computations, these studies covered only an interval of $`10\mathrm{s}`$, adequate to watch the jet propagate, develop a cocoon and break out, but not long enough to study the the expansion after break out to any great extent.

### 4.1 Model Definitions

The parameters of the jet, its power, Lorentz factor and energy loading, were all identical to those of Model 2B, that is $`\mathrm{\Gamma }=5`$, $`\dot{E}=3\times 10^{50}\mathrm{erg}\mathrm{s}^1`$ and $`f_o=40`$ (Table 6). The grid employed in all cases was Cartesian with 256 zones each along the $`x`$\- and $`y`$-axes and 512 along the $`z`$-axis (jet axis) (Table 7).

The initial model and the initiation of the jet in Model 3A were perfectly symmetric with regard to the axis of the jet. The perfect symmetry was maintained in the calculation because our numerical scheme did not break any symmetry. In order to break the perfect symmetry of the cylindrical initial conditions, an asymmetric jet was injected in Models 3BS and 3BL. At the base of the jet in Model 3BS, pressure and density were 1% more than those of Model 3A if $`y\mathrm{tan}\alpha x`$, and 1% less otherwise, where $`\alpha =40`$ degrees. Hence the jet in Model 3BS has a $`\pm 1\%`$ imbalance in power. In Model 3BL, a $`\pm 10\%`$ imbalance in power was employed.

Several other models were calculated to explore the collimation properties of non-radial jets that precessed. These jets were initiated at the same $`10^{10}\mathrm{cm}`$ inner boundary with a half-angle of 5 degrees as measured from the center. However, they were given a non-$`z`$ (symmetry axis) component of momentum that would have resulted, in a vacuum, in a propagation vector inclined to the $`z`$-axis by 3 degrees (Model 3P3), 5 degrees (Model 3P5), and 10 degrees (Model 3P10). The jet precessed around the $`z`$-axis with a period, in all cases of $`2\mathrm{s}`$. This period was chosen to be short compared with the time it took the jet to traverse the remaining star between $`10^{10}\mathrm{cm}`$ and the surface, yet long enough that the jet would still be distinct from a cone.

### 4.2 Breakout in Three Dimensions

As expected, Model 3A closely resembles Model 2T (Fig. 11). The initial stellar model, the zoning (Tables 4 and 7), and the jet parameters (Tables 2 and 6) were all the same. In particular, the initial conditions for the jet have cylindrical symmetry in both cases. It is gratifying that the answer is insensitive to the dimensionality of the grid, which in Model 2T, was two-dimensional cylindrical and, in Model 3A, was three-dimensional Cartesian. This implies that the bulk properties of jets calculated in Paper 1 - including collimation and modulation - would be the same in 3D.

However, Model 3A, by assuming a jet that initially has perfect cylindrical symmetry does not fully exercise the 3D code. Aside from numerical noise, the equations conserve the 2D symmetry of the initial conditions, even on a 3D grid. Fig. 11 also compares, at breakout, the properties of jets that were, initially, nearly identical. In particular, Model 3BS differed from 3A only in a 1% asymmetry in input energy from one side of the jet to the other, yet the structure of the emergent jet and cocoon is strikingly different. In Model 3BL with a 10% asymmetry the difference is even more striking. The plots show a cross section of the Cartesian grid along the initial jet axis in the $`x`$ = 0 plane. Because of the 40 degree offset (§ 4.1), there is somewhat more energy in the top half of the figure than the bottom. Though both jets retain high Lorentz factors in their cores, similar to Models 2T and 2B, the cocoon explosion is a little earlier and larger on the top. More dramatic is the difference in the high density “plug” among Models 3A, 3BS, and 3BL. In the latter two where the 2D symmetry was mildly broken, the plug has a much lower density and is not prominent. In 2D models and symmetric 3D Model 3A, the plug is held by a concave surface of the highly relativistic jet core. Because of the imposed axisymmetry, the plug cannot easily escape and is pushed forward by the jet beam. Whereas the story of the plug is different in asymmetric 3D Models 3BS and 3BL. In these models, where the 2D symmetry was initially mildly broken, instabilities will develop. The forward bow shock is no longer symmetric and even the stellar material on the axis is pushed sideways. More importantly, the head of the highly relativistic jet beam is also asymmetric and does not have a concave surface to “hold” the plug. A movie of these runs shows the plug forming, slipping off to the side, then forming again. The plug in these models has a much lower density and is not prominent. The presence or absence of this plug may have important implication for the production of short-hard gamma-ray bursts in the massive star models (Paper 1; Waxman & Mészáros, 2003).

### 4.3 Stability of the Jet

Given the ability to model jets in 3D, we undertook a study to test their survivability against non-radial instabilities. Three jets were introduced on the standard grid (Table 7) with inclinations to the radial of 3, 5, and 10 degrees (Table 6). These jets were made to precess with a period equal to $`2\mathrm{s}`$, i.e., a non-trivial fraction of the time it takes the jet to traverse the star and break out. The results are shown in Fig. 12.

With an angle of 3 degrees, the jet in Model 3P3 escapes the star with its relativistic flow at least partly intact. Though Fig. 12 shows Lorentz factors of only $`\mathrm{\Gamma }1020`$, the internal energy is still very large and some of these regions will attain $`\mathrm{\Gamma }180`$ if they expand freely. However, there is some “baryon-poisoning”. No clear line of sight exists along a radial line for the most energetic material even in the jet beam. Perhaps a jet with intermediate Lorentz factor would emerge. We have not yet followed these calculations to large radius (as in Figs. 2, 3, and 4) because the 3D grid was not large enough. Larger calculations are planned.

At larger angles Model 3P5 and especially 3P10 show the break up of the jet. Because there is no well focused highly relativistic jet beam, especially in Models 3P5 and 3P10, more baryon mass is mixed into the jet. The Lorentz factor of the jet will decrease due to dissipation. And it will be very difficult for these jets to make a common GRB. Again some sort of hard transient might be expected, especially after all the internal energy converts and $`\mathrm{\Gamma }`$ rises, but for Model 3P10 there will be no GRB of the common variety. In our simulations, we find that the critical angle for jet precession is about 3 degrees. One would expect that the constraint on the angle of precession will be reduced if the jet bears more power or is powered longer.

Is two seconds a reasonable period for the jet to precess? There are many uncertainties, but the gravitomagnetic precession period of a black hole surrounded by an accretion disk is estimated by Hartle, Thorne, & Price (1986) to be

$$T_{\mathrm{GM}}=900(M_{\mathrm{disk}}/0.001\mathrm{M}_{})^1(M_{\mathrm{bh}}/3\mathrm{M}_{})^{1/2}(r/300\mathrm{km})^{2.5}\mathrm{s},$$

where $`M_{\mathrm{disk}}`$ is the mass of the accretion disk, $`M_{\mathrm{bh}}`$ is the mass of the hole and $`r`$ is the radius of the disk. $`r=300\mathrm{km}`$ is a reasonable value for the radius of the disk, which depends on the distribution of angular momentum of the star (MacFadyen & Woosley 1999). The dependence of the mass in the disk on the alpha-viscosity has been explored in both numerical simulations (MacFadyen & Woosley 1999) and analytic calculations (Popham, Woosley & Fryer 1999). For currently favored value of the alpha-viscosity, $`0.1`$, the disk mass is about $`0.001\mathrm{M}_{}`$. For the disk mass to approach $`1\mathrm{M}_{}`$ the viscosity would need to be $`\alpha <0.001`$. Given the low mass of the accretion disk and the radius of the disk, $`300\mathrm{km}`$, precession periods as small as two seconds are unlikely, but the black hole may be kicked or instabilities deeper in the region of the star not modeled here could give the jet some non-radial component. If so, Fig. 12, shows an alternate way in which baryon-loaded jets could give rise to less relativistic mass ejected and slower moving jets far away from the star($`10\mathrm{\Gamma }100`$). Softer transients like XRFs and sub-luminous GRBs like GRB 980425 could result.

## 5 Discussion

Our calculations show, so long as the orientation of the jet does not waiver on time scales of order seconds by angles greater than 3 degrees, that a relativistic jet can traverse a Wolf-Rayet star while retaining sufficient energy and Lorentz factor to make a GRB. This conclusion is robust in three dimensions as well as two.

As it breaks out, the jet is surrounded by a cocoon of mildly-relativistic, energy-laden matter. As internal energy converts to kinetic, matter with $`10^{51}`$ to $`10^{52}\mathrm{erg}`$ of equivalent isotropic energy moving with Lorentz factor $`\mathrm{\Gamma }20`$ is ejected to angles about three times greater than the GRB. That is, whatever transient the cocoon gives rise to will be an order of magnitude more frequently observable in the universe, but two orders of magnitude less energetic than a GRB. Weaker transients can of course be obtained at still larger angles and there is considerable diversity in the models (Figs. 5, 6, and 7).

The isotropic equivalent rest mass to $`10^{51}`$ \- $`10^{52}\mathrm{erg}`$ with $`\mathrm{\Gamma }20`$ is $`3\times 10^5`$ to $`3\times 10^4\mathrm{M}_{}`$. This matter would give up its energy upon encountering $`1/\mathrm{\Gamma }`$ times its rest mass. For an assumed mass loss rate for the Wolf-Rayet progenitor of $`10^5\mathrm{M}_{}\mathrm{yr}^1`$ and speed, $`1000\mathrm{km}\mathrm{s}^1`$, the energy will be radiated at about $`3\times 10^{14}`$ to $`3\times 10^{15}\mathrm{cm}`$. Emission from this deceleration will have a duration $`R/(2\mathrm{\Gamma }^2c)`$ 10 - $`100\mathrm{s}`$. Variations of a factor of 10 are easily achievable by varying the mass loss rate or Lorentz factor.

Unlike the central jet, we do not see evidence for large scale variation in the Lorentz factor at lower latitude and so it would be hard to make the transient emission by internal shocks. We thus attribute the emission to external shocks and invoke a variety of Lorentz factors emitting within our light cone in order to explain the spread in wavelength.

Might there be observable counterparts to these large angle, low Lorentz factor explosions? Typically, they have too low a Lorentz factor to make common GRBs (Guetta, Spada, & Waxman, 2001). Tan, Matzner, & McKee (2001) present arguments however that, by external shock interaction with the progenitor wind, a hard transient of some sort should result.

Our model predicts a correlation between $`E_{\mathrm{iso}}`$, the Lorentz factor, and the angle between the polar axis and the observer (Fig. 5, 6 and 7). Roughly speaking, the GRB outflows have a narrow highly relativistic jet beam and a wide mildly relativistic jet wing. Recent observations and afterglow modeling support this non-uniform jet model (Berger et al., 2003; Zhang et al., 2003; Salmonson, 2003; Rossi, Lazzati, & Rees, 2002). If, as seems reasonable, the Lorentz factor is, in turn, correlated with the peak energy observed in the burst one expects a continuum of high energy transients spanning the range from X-ray afterglows (keV), to hard X-ray transients (tens of keV), to GRBs (hundreds of keV). Observations reviewed by Amati et al. (2002) show such a correlation for bursts with $`E_{\mathrm{peak}}`$ from 80 keV to over 1 MeV and Lamb (private communication) finds that the relation extends to the lower energies.

In particular, XRFs form a new class of X-ray transients having a duration of order minutes and properties that in ways resemble GRBs (Heise et al, 2001; Kippen et al., 2003). XRFs are probably many phenomena (Arefiev, Priedhorsky, & Borozdin, 2003) and it could be that, especially some of the longer ones, have alternate explanations (see introduction). However, we have felt for some time (Woosley, Eastman, & Schmidt, 1999; Woosley & MacFadyen, 1999; Woosley, 2000, 2001, Paper 1) that many XRFs are the off-axis emissions of GRBs, made in the lower energy wings of the principal jet. We have not calculated the expected spectra for any of our models, only Lorentz factors and energies. However, it is reasonable that matter moving with $`\mathrm{\Gamma }`$ 20 would make a transient softer than a GRB and harder than a few keV. Larger amounts of energy are emitted at a larger range of angles for slower but still relativistic ejecta ($`\mathrm{\Gamma }=5`$, e.g., Fig. 5).

If our model is valid, XRFs and GRBs should be continuous classes of the same basic phenomenon sharing many properties. They should have a similar spatial distribution to GRBs because they are essentially the same sources. However, because they are much less luminous, their log N-log S distribution would not exhibit the same roll over attributed in GRBs to seeing the “edge” of the distribution (e.g., Fishman & Meegan 1995). Their median redshift should be considerably smaller, certainly less than 1. This is consistent with the low redshifts inferred for the host galaxies of two XRFs by Bloom et al. (2003). XRFs may, most frequently, be seen in isolation and will be characterized by softer spectra, but there would also be an underlying XRF in every GRB since the emission of the mildly relativistic cocoon material is beamed to a larger angle that includes the poles. In some cases these XRFs might be seen as precursors or extended hard X-ray emission following a common GRB. They should be associated with supernovae. Indeed XRFs may more frequently serve as guideposts to jet-powered supernovae than GRBs, especially the nearby ones.

Though considerable variation is expected, our calculations (e.g., Fig. 5) suggest that XRFs are typically visible at angles about three times greater than GRBs and hence to ten times the solid angle. However, their energy, a few percent of GRBs, implies that a flux-limited sample could observe GRBs out to roughly ten times farther, implying that XRFs would be about 1% as frequent in the sample. The actual value is detector sensitive, but may be larger.

One should also keep in mind the possibility that XRFs do not accompany GRBs, in any direction, but are a result of baryon loading of the central jet. Three-dimensional instabilities (Fig. 12) may play an important role in this and are being explored.

This research has been supported by NASA (NAG5-8128, NAG5-12036, and MIT-292701) and the DOE Program for Scientific Discovery through Advanced Computing (SciDAC; DE-FC02-01ER41176). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. AH was supported, in part, by the Department of Energy under grant B341495 to the Center for Astrophysical Thermonuclear Flashes at the University of Chicago, and a Fermi Fellowship at the University of Chicago. We appreciate the helpful comments of the referee.


format: text

text:
# Epidemics and Dimensionality in Hierarchical Networks

## I Introduction

Epidemics of all kinds are costly. The influenza in the early 20th century and the outbreak of SARS in 2003 claimed precious lives, and the spreading of viruses on the internet has led to much economic loss when computer systems of governments and airlines were attacked. Recent breakdowns of major power networks in North America and in Europe also brought great loss and inconvenience to millions of people. A less damaging example is, perhaps, that of the spread of a rumor in a community. Epidemiological processes, together with percolation and problems on searchability, are the most interesting dynamical problems in a system of connected entities diekmann ; albert ; dorogovtsev ; dorogovtsev1 . There are two key factors in the study of epidemic spreading, namely how infection and recovery occur and how people or computers are connected. The latter is best described in terms of networks or graphs albert ; dorogovtsev ; dorogovtsev1  in which the nodes representing people, communities, power plants, or computers are connected by links representing the acquaintance among people or the connections between computers. The most important question is the extent to which a disease (virus) may propagate. Since infections take place through direct contact between infected (I) and susceptible (S) nodes, the underlying structure of the network understandably plays a determining role on the spread of diseases Satorras .

The Susceptible-Infected-Refractory (or recovered/removal) (SIR) model diekmann ; anderson ; murray  has been widely used in different forms for studying epidemiological processes such as the spread of influenza (see, for example, Ref. liu ) and SARS riley ; lipsitch , and also for the spread of rumors zanette ; zanette1 . In its simplest form, a node takes on one of three states, namely S, I or R. The R-state represents the case in which an infected node becomes recovered and immunized or, unfortunately dead. For computer networks, the R-state may represent recovery with the virus removed and anti-virus software installed. A node in the R-state is not infectious and will not become infected again. The underlying connectivity obviously plays a decisive role. For nodes connected randomly, a substantial proportion becomes infected and eventually recovered from an epidemic sudbury . However, when a regular network is gradually re-wired into a small-world network albert ; dorogovtsev ; dorogovtsev1 ; watts1 ; watts2  and then a random network, it was found that a threshold exists on the degree of disorderness below which the disease can only spread locally, but not globally zanette ; zanette1 .

Networks have fascinating structural and dynamical properties, and they have wide applications in physical, biological, social, economic, and computer systems albert ; dorogovtsev ; dorogovtsev1 ; strogatz ; newman . Physicists have made important contributions to the understanding of networks in recent years. Structurally, the small-world networks watts1 ; watts2  and the growing scale-free networks barabasi  have been extensively studied albert ; dorogovtsev ; dorogovtsev1 . Depending on the model, the degree distribution may show power law or exponential forms. As the underlying connectivity plays a decisive role in the studies of dynamical processes, a good network model that reflects the structural properties is therefore required for modelling epidemics or other dynamical problems in a connected population, or in systems in which the entities are grouped in a hierarchical fashion as in a population. The hierarchical model of searchable social networks proposed recently by Watts et al. watts3  captures the essential ingredients of a network model of connected population and has enjoyed much attention. The model is physically transparent in that it is motivated by the observation that individuals within a population are grouped according to their function into different dimensions, for example, their occupation, hobby, home district, etc. The model exhibits the phenomena of “six degrees of separation”, which was discovered by Travers and Milgram in the 1960’s travers  and has been tested in a world-wide experiment on the internet watts4 .

Recent works on networks on robustness of networks under attacks shargel  and ways to prevent attacks cohen  indicate that, again, the connectivity of a network is important. It is, therefore, important to study the dynamical aspects of the successful social network model of Watts et al. watts3 . In the present work, we present detailed results of numerical simulations on the extent of a spread of disease or virus in the hierarchical social network model of Watts et al. within the SIR epidemic dynamics. It is found that the spread may be local or global if the population can be partitioned in a unique way. However, when the partition of the population becomes multi-dimensional, the result is different from that of the searchability problem. It is found that the spread is always global.

The plan of the paper is as follows. In Sec.II, we motivate the hierarchical structure and describe how links are established between nodes in a hierarchical structure with a tunable degree of homophily. In Sec.III, we describe the implementation of the SIR dynamics and discuss the results on the spread of an epidemic in populations that can be classified into one and two hierarchies. Section IV summarizes our findings.

## II Network model of connected population

The model of Watts et al. was motivated by the general structure in the groupings of individuals in a society as shown in Fig.1 watts3 . The highest level can be regarded as a population of $`N`$ individuals or nodes. These $`N`$ nodes may then be classified or partitioned into $`b`$ groups. Each group can further be divided into $`b`$ subgroups and so on. After $`(l1)`$ divisions the structure has a total of $`l`$ levels and ends at a level where an individual belongs to a close functional group of size $`g`$, where $`g`$ typically is of the order of $`10^1`$ to $`10^2`$. Individuals belonging to the same lowest-level subgroup have the highest chance of becoming most similar or getting acquaintance. For example, all physicists in academia can be classified roughly by their research area, say condensed matter physics or nuclear physics, and further classified by their subfields, e.g., magnetism or superconductivity, and then further grouped by their specific research interests, e.g., topics grouped into the same session in a conference or with the same code in the PACS classification scheme. Obviously the divisions are usually not unique. Physicists may be grouped geographically based on the region their institution is located, and then the country, state, county, department, and research groups. A population of individuals or a collection of nodes may thus be characterized by $`H`$ hierarchies or dimensions, each of which takes on the structure shown in Fig.1. The structure of the network of $`N`$ nodes is then characterized by the parameters $`H`$, $`b`$, $`l`$, and $`g`$ with $`N=gb^{l1}`$, where $`g`$ is the average size of the lowest-level functional subgroups.

An important quantity within a hierarchy is the social distance $`x_{ij}`$ that measures the similarity between two nodes $`i`$ and $`j`$. For nodes belonging to the same lowest-level subgroup in a given hierarchy, $`x_{ij}=1`$, otherwise $`x_{ij}`$ is the number of levels from the lowest for which the two nodes belong to the same group. Hence the largest separation is $`l`$ in a given hierarchy. For $`H>1`$, an important geometrical feature is that for nodes $`i`$ and $`j`$ with $`x_{ij}=1`$ in one hierarchy and $`j`$ and $`k`$ with $`x_{jk}=1`$ in another, $`x_{ik}`$ is in general not small. For example, it is unlikely that the colleague in the office next door knows your collaborator in another country. Within the context of epidemics in a community, the hierarchies may be regarded as family, relatives, friends, and so on; or neighbors in the same apartment building, same neighborhood, same district and so on. While the model was motivated by the structure in groupings within a society, it is obvious that many well-structured systems, e.g. computer clusters, are connected in a similar fashion.

Computationally, a social network with a tunable degree of homophily can be constructed as follows watts2 . For a population of $`N`$ individuals with social structure as shown in Fig.1, the individuals are first assigned randomly into the lowest subgroups with an average size $`g`$ in each of the $`H`$ hierarchies. Links between individuals are established as follows. For $`H>1`$, a hierarchy is randomly chosen. Two nodes $`i`$ and $`j`$ are then selected randomly within the chosen hierarchy. A link, specifying that $`i`$ and $`j`$ are friends and hence become acquaintance, is established with a probability $`P(x_{ij})=\mathrm{exp}(\alpha x_{ij})/_{n=1}^l\mathrm{exp}(\alpha n)`$ depending on $`x_{ij}`$ of the chosen nodes in the chosen hierarchy. No duplicated links between $`i`$ and $`j`$ are allowed. Here, the parameter $`\alpha `$ is a measure of homophily of the system. The process is then repeated until a mean number of $`z=g1`$ links are established for each node (individual) in the system. This guarantees that for $`\alpha \mathrm{ln}b`$ and $`H=1`$, only links between nodes with small separation are probable and the individuals are connected only to those most similar in character, leading to isolated subgroups of nodes. For $`\alpha =\mathrm{ln}b`$, links between individuals with any social distance are equally probable and a random network results, with the notion of similarity between nodes loses its meaning. For intermediate values of $`\alpha `$, i.e., $`\alpha 1`$, the network shows small-world features. Here, we aim at studying the extent of an epidemic as a function of $`\alpha `$ for systems with $`H=1`$ and $`H=2`$, respectively.

## III Epidemic Modelling and Results

We model epidemics on the social network model by the standard SIR dynamics, which is implemented as follows anderson ; murray . Initially, all nodes are in the S-state and one node is randomly chosen for infection. At each timestep, a node $`i`$ is randomly chosen among all the infected (I) nodes. A neighbor (friend) $`j`$ is then selected randomly among all the neighbors of node $`i`$, i.e., those with a link connected to $`i`$. Note that the effects of dimensionality $`H`$ is built in through the construction of links. If node $`j`$ is susceptible, it becomes infected and the chosen node $`i`$ remains in the I-state; otherwise (i.e., node $`j`$ is either I or R) the state of node $`j`$ remains unchanged and the chosen node $`i`$ becomes recovered (R) at the end of the timestep. As time evolves, the number of R-nodes (S-nodes) increases (decreases); while the number of I-nodes increases initially and then eventually drops to zero. The number of R-nodes at the end of the epidemic is denoted by $`N_R`$.

We performed extensive numerical simulations to explore the effects of the structural parameter $`\alpha `$ and the dimensionality or number of hierarchies $`H`$. For each set of parameters, data are collected over $`10^3`$$`10^4`$ runs, each of which corresponds to a different realization of the network and a different initially infected node. We take $`b=2`$, $`g=10`$ and hence $`z=g1=9`$. Results for $`H=1`$ are shown in Fig.2. The structural effects on epidemics can be studied via the probability distribution function $`P(N_R)`$ for different values of $`\alpha `$ and different population sizes $`N`$. For $`\alpha =10\mathrm{ln}2`$, $`P(N_R)`$ is independent of $`N`$ and peaks at $`N_R8`$, with no network having $`N_R>25`$ (see Fig.2(a)). The independence on the population size $`N`$ implies that the spread is local, for cases of large $`\alpha `$. This is the case of “regular networks” consisting of isolated subgroups in which only local spread at the lowest level is possible. Each subgroup can be regarded as a small fully connected network. Previous studies showed that about 79% of nodes in a large fully-connected network eventually recovered sudbury . Using this result as an estimation gives $`0.79g8`$ for the size of a spread, in good agreement with the peak value observed. As only the mean $`g`$ is fixed, each cluster has a different $`g`$ leading to a distribution $`P(N_R)`$ that extends to $`N_R=25`$. Simulations using different values of $`g`$ reveal a corresponding shift in $`N_R`$ at which $`P(N_R)`$ shows a maximum. For the small-world regime corresponding to $`\alpha 1`$ (see Fig.2(b)), $`P(N_R)`$ is again $`N`$-independent and decays exponentially with runs in which the spread can be up to $`N_R10^3`$. Geometrically, there exist some links of longer social distance as $`\alpha `$ decreases, leading to larger but still localized spreads. As $`\alpha `$ decreases towards the random network limit of $`\alpha =ln2`$, $`P(N_R)`$ takes on the form of two disjoint parts (see Fig.2(c)): (i) one that decays exponentially near $`N_R0`$ (inset) characterizing the small fraction of runs with small-scale epidemics due to the existence of isolated clusters and (ii) a $`N`$-dependent gaussian distribution characterizing global spreading centered at a value $`\overline{N_R}`$ that scales linearly with $`N`$. Writing $`\overline{N_R}=a(g)N`$, we found $`a0.66`$ for $`g=10`$. Extensive simulations reveal that $`a(g)`$ approaches $`0.79`$, the value for fully-connected networks, as $`g`$ increases. For $`H=1`$, the consequence of an epidemic depends sensitively on the structural parameter $`\alpha `$.

Social networks usually correspond to $`H>1`$. For $`H>1`$, the characteristic features of an epidemic is very different from that in $`H=1`$. The most striking result is that for $`H1`$, an epidemic can always spread globally regardless of the degree of homophily (i.e., value of $`\alpha `$) in the network. Figure 3 shows the results for $`H=2`$. For fixed $`N`$ and $`z`$, $`P(N_R)`$ becomes independent of $`\alpha `$ and takes on a form identical to that in Fig.2(c). The majority of the weight is represented by the gaussian distribution peaked at $`\overline{N_R}a(g)N`$, with $`a(g=10)0.66`$ and increases for larger $`g`$. This result should be contrasted with that of searchability on the same social network watts3 , which occurs for $`H>1`$ and $`\alpha >0`$ with increasingly restricted constraints on $`H`$ and $`\alpha `$ as $`N`$ increases. For epidemics, $`H>1`$ alone is sufficient for a global spread of the disease and the result remains valid for different values of $`N`$. The independence on $`\alpha `$ can be understood by the way that the links are constructed. Although a larger $`H`$ implies fewer links per hierarchy, the chance of nodes $`i`$ and $`j`$ being similar with a shorter $`x_{ij}`$ in one of the hierarchies is higher and so is the chance that a link exists between nodes $`i`$ and $`j`$, even for large values of $`\alpha `$. In other words, nodes $`i`$ and $`j`$ unknown to each other in one hierarchy has a chance of becoming connected when additional dimension(s) is (are) added. Thus an increase in $`H`$ has the effect of promoting the spread of a disease in a population by permeating through different connected nodes in different hierarchies. This is consistent with daily experience that a flu may spread first from a parent to a kid (linked nodes in a family) and then between the kids in school (linked nodes in school) and so on. We have also studied the case of $`H=3`$ and obtained similar results.

To further illustrate the qualitative difference between $`H=1`$ and $`H=2`$ networks, we show the mean fraction of R-nodes, $`rN^1_{N_R=0}^NN_RP(N_R)`$, as a function of $`\alpha `$ in Fig.4. For $`H=1`$, the large $`\alpha `$ ($`\alpha 1`$) regime gives $`rg/N0`$ corresponding to the case in Fig.2(a), while the small $`\alpha `$ ($`\alpha <\mathrm{ln}b`$) regime gives an $`N`$-independent value of $`r0.66`$ corresponding to the case in Fig.2(c). The fraction $`r`$ drops sharply in a small intermediate range of $`\alpha `$ ($`\mathrm{ln}2<\alpha <2`$) corresponding to the small-world regime watts3 . This feature is similar to that observed in the spreading of rumors in a small-world construction zanette ; zanette1 . In this regime, links between nodes with large dissimilarity exist and hence $`N_R/N`$ increases as $`\alpha `$ decreases. One would expect that the transition becomes sharper as $`N`$ increases, as $`N_R`$ is restricted only by the connections and is insensitive to $`N`$ in the intermediate regime. For $`H=2`$, the behavior is qualitatively different in that $`r`$ stays at a higher level for all values of $`\alpha `$ and the transition disappears as a result of the effect of the extra hierarchy for hooking up individuals in a population.

## IV Summary

In summary, epidemiological processes are studied on a hierarchical social network within the SIR model. For population characterized by $`H>1`$ independent hierarchies, global spreading results regardless of the structural parameter $`\alpha `$ and hence the homophily of individuals forming a social circle. This result should be contrasted with that of searchability problems in that $`H>1`$ alone is sufficient for a global spread of a disease within the SIR dynamics. For $`H=1`$, a transition from global to local spread occurs as the population becomes decomposed into increasingly homophilous groups when the structural parameter $`\alpha `$ increases. Since social networks usually correspond to $`H>1`$, it is therefore very difficult to confine an epidemic, unless the “dimensionality” of the infected nodes can be reduced at least temporarily. One way to achieve this is to quarantine the infected nodes, which is known to be an effective way in handling infectious diseases. It will also be interesting to extend the present study to other models of epidemics such as the SIS Satorras  and SIRS kuperman  models.

###### Acknowledgements.

This work was partly supported by the National Natural Science Foundation of China under Grant No. 70471081 and No. 70371069, by a DFG grant under Grant No. TR 3000/3-3, and by a grant from the Research Grants Council of Hong Kong SAR Government under Grant No. CUHK4241/01P.

FIGURE CAPTIONS

Figure 1: Schematic diagram of the groupings of individuals in a hierarchical social network. A population of $`N`$ nodes are classified into $`b`$ groups. Each group is further divided into $`b`$ subgroups and so on. After ($`l1`$) divisions, the lowest-level subgroups consist of $`g`$ individuals. There $`x_{ij}`$ is the social distance between nodes $`i`$ and $`j`$. A population can, in general, be characterized by $`H`$ such hierarchies.

Figure 2: The probability distribution function $`P(N_R)`$ of the number of R-nodes at the end of an epidemic in a population characterized by $`H=1`$ hierarchy. The structural parameter takes on (a) $`\alpha =10`$, (b) $`\alpha =1`$, (c) $`\alpha =\mathrm{ln}2`$. Different population sizes $`N=1280`$ (circles), $`2560`$ (squares), and $`5120`$ (triangles) corresponding to $`l=8,9,10`$ are studied. Each data point represents an average over $`10^4`$ runs. The inset in (c) shows the portion of $`P(N_R)`$ for small $`N_R`$.

Figure 3: The probability distribution function $`P(N_R)`$ in a population of $`N=2560`$ ($`l=9`$) characterized by $`H=2`$ hierarchies. The structural parameter takes on $`\alpha =10`$ (circles), $`\alpha =1`$ (squares), and $`\alpha =\mathrm{ln}2`$ (triangles). Each data point represents an average over $`10^4`$ runs. $`P(N_R)`$ signifies a global spread of disease regardless of the structural parameter for $`H>1`$.

Figure 4: The fraction of R-nodes at the end of an epidemic as a function of the structural parameter $`\alpha `$ in a population of $`N=1280`$ (circles), $`2560`$ (squares), and $`5120`$ (triangles) characterized by (a) $`H=1`$ and (b) $`H=2`$ hierarchies. Each data point represents an average over $`10^3`$ runs. For $`H>1`$, global spreading occurs for all values of $`\alpha `$, while a transition from global to local spreading occurs for $`H=1`$.


text:
# A New Framework for the Assessment and Calibration of Medium Range Ensemble Temperature Forecasts

## 1 Introduction

A number of different methods have been used for the assessment and calibration of ensemble forecasts (for example see Atger (1999), Richardson (2000), Roulston and Smith (2002) and Wilson et al. (1999)). In many applications of ensemble forecasts the forecast is used to derive the probability of a certain outcome, such as temperature dropping below zero. In this context, the reliability diagram is an appropriate method for assessing reliability (see Anderson (1996), Eckel and Walters (1998), Talagrand et al. (1997) and Hamill (1997)) and the relative operating characteristic (ROC) (Mason (1982), Swets (1988), Mason and Graham (1999)) is an appropriate method to evaluate the resolution.

In other applications of ensemble forecasts, however, the forecast is interpreted as providing a mean and a distribution of future values of temperature. For example in the field of weather derivatives the calculation of the fair strike for a certain class of weather swap contract222linear swaps on a linear temperature index needs an estimate of the mean of the future temperatures, while the calculation of the fair premium for weather option contracts needs an estimate of the whole distribution of future temperatures (see Jewson and Caballero (2002) for details). Additionally, the assumption is often made that temperature is normally distributed since this allows the temperature forecast to be summarised succinctly using just the mean and the standard deviation. For such mean-and-distribution or mean-and-standard deviation based applications of ensemble forecasts the reliability diagram and the ROC are not particularly appropriate.

In this paper we present a new framework for the assessment and calibration of ensemble temperature forecasts based on analysis of the mean and standard deviation of the distribution of temperatures. The method has been developed to respond to the need for a simple and practical method for assessment and calibration that can be used by companies that make use of ensemble forecasts in the weather derivative market. We postulate a parametric model for the mean and standard deviation and fit the parameters of the model using the maximum likelihood method. This approach has a number of advantages relative to the assessment and calibration methods mentioned above. The model is simple, easy to interpret, and the entire ensemble distribution can be calibrated in one simple step. Also the model gives a clear indication of how many days of useful information there are in a forecast.

In section 2 we describe the data sets we use for this study. In section 3 we describe the statistical model that forms the basis for the method we propose. In section 4 we describe the results from fitting the model. In section 5 we discuss extensions to other distributions and in section 6 we summarise our results and draw some conclusions.

## 2 Data

We will base our analyses on one year of ensemble forecast data for the weather station at London’s Heathrow airport, WMO number 03772. The forecasts are predictions of the daily average temperature, and the target days of the forecasts run from 1st January 2002 to 31st December 2002. The forecast was produced from the ECMWF model (Molteni et al., 1996) and downscaled to the airport location using a simple interpolation routine prior to our analysis. There are 51 members in the ensemble. We will compare these forecasts to the quality controlled climate values of daily average temperature for the same location as reported by the UKMO.

Throughout this paper all equations and all values have had both the seasonal mean and the seasonal standard deviation removed. Removing the seasonal standard deviation removes most of the seasonality in the forecast error statistics, and justifies the use of non-seasonal parameters in the statistical models for temperature that we propose.

## 3 The Moment-based Ensemble Assessment and Calibration Model

For forecasts of temperature anomalies, it has long been recognized (see for example Leith (1974)) that the use of a final regression step between ensemble mean and observations can eliminate bias and minimise the mean square error (MSE). For our purposes we will write this regression step as:

$$T_iN(\alpha +\beta m_i,\sigma )$$
(1)
where $`T_i`$ is the observed temperature on day $`i`$, $`N(\mu ,\sigma )`$ represents a normal distribution with mean $`\mu `$ and standard deviation $`\sigma `$, $`m_i`$ is the forecast of the temperature (in our case, the ensemble mean) and $`\alpha `$, $`\beta `$ and $`\sigma `$ are free parameters. This regression model postulates that temperatures come from a normal distribution with mean given by $`\mu _i=\alpha +\beta m_i`$ and standard deviation given by $`\sigma `$. The values for $`\alpha `$, $`\beta `$ and $`\sigma `$ come from fitting the model, and this is usually done using least squares linear regression. One justification for the use of least squares linear regression is that for this particular model it is equivalent to finding the parameters that maximize the likelihood of the data given the model (see Press et al. (1992)), as long as we assume that the forecast errors are uncorrelated in time. We note that although the model in equation 1 postulates that the data come from a normal distribution, it can be applied in situations in which the data is not strictly normal, and in fact it is common (although perhaps bad) practice not to test for normality when doing such linear regressions.

One of the assumptions in this model is that the standard deviation of the forecast errors $`\sigma `$ is constant. However it is well documented that the size of forecast errors varies in time (Palmer and Tibaldi, 1988) and that there is a relationship between the ensemble spread and the size of forecast errors (Toth et al., 2000). It thus makes sense to attempt to generalize the model in equation 1 to a model that takes these temporal variations in $`\sigma `$ into account. We will do this using the model:

$$T_iN(\alpha +\beta m_i,\gamma +\delta s_i)$$
(2)
where the free parameter $`\sigma `$ has been replaced by a linear function of the ensemble spread $`s_i`$, and two new parameters $`\gamma `$ and $`\delta `$ have been introduced. Modelling the standard deviation as a linear function of the ensemble spread in this way allows for both time variation and the correction of biases in the predicted uncertainty.333We note that one could alternatively model the variance as a linear function of the spread squared.

The optimum parameters for this model can no longer be fitted using least squares linear regression. However, they can be fitted if we can identify a cost function that can be minimised or maximised by varying the parameters. There are various possibilities for such a cost function, but one of the most natural is the likelihood, defined as the probability density of the observations given the calibrated forecast. Maximising the likelihood is the standard way to fit parameters in statistics (see for instance textbooks such as Casella and Berger (2002) or Lehmann and Casella (1998)), and gives the most accurate possible estimates of the parameters for most statistical models.

As with the linear regression model, this model is also not restricted to cases in which temperature is exactly normally distributed: the assumption of the normal distribution merely provides a metric in which the likelihood can be calculated and the parameters fitted. This metric is most appropriate when the data is at least close to normally distributed. For cases when the data is not close to normal other distributions can be used, or the data can be transformed to normal.

There are a number of useful features of the model we present. These include:

* Once the parameters have been fitted to past historical data, calibration of future ensemble forecasts is easy since it just involves applying linear transformations to the ensemble mean and standard deviation. The calibrated values for the mean and the standard deviation can be used to define the whole forecast distribution, or can be used to shift and stretch the individual ensemble members, if individual ensemble members need to be preserved. In the latter case non-normality in the distribution of the original ensemble members will not be destroyed.
* The optimum values of the parameters in equation 2 have clear interpretation and give us useful information about the performance of the ensemble. For instance $`\alpha `$ identifies a bias in the mean, and $`\beta `$ represents a scaling of the forecast towards climatological values. In a perfect forecast, $`\alpha `$ would be zero and $`\beta `$ would be one. The spread parameters $`\gamma `$ and $`\delta `$ combine to optimize the prediction of uncertainty about the mean. The value of the ensemble spread $`s`$ varies in time because of the dependence of the growth rate of differences between ensemble members on the actual model state. The calibrated standard deviation value $`\sigma _i=\gamma +\delta s_i`$ additionally includes uncertainty due to model error. If the spread of the ensemble contains very little real information, $`\delta `$ will tend to be small, and $`\gamma `$ will tend to be large to compensate.
* It is very easy to calculate approximate uncertainty levels on the values of the parameters as part of the fitting procedure. This is done using the curvature of the log-likelihood at the maximum (see the above references on likelihood methods). These uncertainty levels give us a clear answer to the question of whether the ensemble forecast has useful skill at different lead times. For instance, once $`\beta `$ is not significantly different from zero we can say that the ensemble mean no longer contains useful information (at least not within this framework) and once $`\delta `$ is not significantly different from zero then we can say that the ensemble spread no longer contains useful information. This raises the interesting possibility that we might identify situations in which the mean may contain more days of useful information than the spread.
* It is often necessary to decide which of two forecasts is the more accurate. If two forecasts are both calibrated using equation 2 then the log-likelihood provides a natural way to compare the forecasts. Log-likelihood measures the ability of the forecast to represent the whole distribution of observed temperatures, and is a generalisation of mean square error. It can be presented in a number of ways such as log-likelihood or log-likelihood skill score.

Forecasts calibrated using equation 2 will not necessarily minimise MSE. Users interested purely in a single forecast that minimises MSE should thus calibrate using equation 1. However, users interested in predictions of uncertainty, or, equivalently, in the whole distribution of possible temperatures, should calibrate using equation 2. In practice we have found that the mean temperature prediction produced by equation 2 is close to that produced by equation 1, presumably because the fluctuations in uncertainty are not large.

## 4 Results

The optimum values for the parameters in equation 2 for our 1 year of forecast data and observations are shown in figure 1. In each case we show the approximate 95% sampling error confidence intervals around the optimum parameters. In some cases they are so narrow that they are hard to see in the graphs.

Looking at $`\alpha `$ we see that there is a small and roughly constant bias in the temperatures produced by the ensemble. Correction of the ensemble mean (or each ensemble member) using $`\alpha `$ would eliminate this bias, as long as the ensemble stays stationary.

The parameter $`\beta `$ is slightly less than 1 at all leads. This shows that the ensemble mean varies too much: either the ensemble mean, or each ensemble member, should be reduced by the factor $`\beta `$ towards the climatology. Such a damping factor is presumably required because the ensemble members are more correlated with each other than they are with the observations and because the ensemble is finite in size. Even at lead 10 $`\beta `$ is highly significantly different from zero, implying that the ensemble mean still contains useful predictive ability at that lead. If we allow ourselves to extrapolate the $`\beta `$ curve to longer leads by eye, it would seem likely that the ensemble mean would still contain useful predictive information even beyond that.

The fact that our values of $`\delta `$ are significantly different from zero out to the end of the forecast (just) shows that there is significant information in the ensemble spread too. However, in this case if we extrapolate to higher lead times by eye it seems unlikely that there would be any more skill in $`\delta `$. Since $`\delta `$ is below one and $`\gamma `$ is non-zero we see that the standard deviation of the ensemble is not an optimal estimate of the uncertainty of the prediction.

The $`\gamma +\delta s`$ transformation can change both the mean spread (the time mean of the standard deviation across the ensemble) and the variability of that spread (the standard deviation in time of the standard deviation across the ensemble). To measure the effect on the mean spread, figure 2 shows values of $`\frac{\gamma +\delta \overline{s}}{\overline{s}}`$ (where the overbar indicates the mean in time over the year of data) which shows the factor by which the transformation increases the mean spread. We see that at short lead times, the ensemble spread $`s`$ is far too small on average and the calibration increases the spread by factors of around 4 (at lead 0) and 2 (at lead 1). At longer lead times the ensemble spread is still too small on average by a factor of around 1.2. This underestimation of the spread from ensemble forecasts has been noted by a number of authors such as Ziehmann (2000) and Mullen and Buizza (2000). It is likely to be due to model error in the prediction model and due to the fact that the forecast is a prediction of a large scale flow while the observation is site-specific and hence affected by small-scale variability not represented in the model.

The size of the effect of the calibration on the *variability* (in time) of the spread is given by the value of $`\delta `$. Since $`\delta `$ is significantly different from one at all lead times beyond the first we conclude that the variability of the spread from the ensemble needs to be reduced to be optimum at those lead times. This could be because the variability of the ensemble spread is too large, or because the variability of the ensemble spread is not highly correlated with the real variability of skill.

We can see from the values of $`\delta `$ that the variability of the ensemble spread alone will overestimate the state-dependent predictability of this model by a large factor at long leads. A better estimate for the level of state-dependent predictability is given by the variability of the calibrated spread, which is smaller by the $`\delta `$ factor.

Figure 3 shows the ratio of the standard deviation of the ensemble spread to the mean ensemble spread at different lead times. We call this ratio the coefficient of variation of the spread (COVS). Figure 3 shows the COVS estimated from both the uncalibrated and the calibrated ensemble data. These values give an indication of how much extra information we get about the forecast uncertainty by using the (uncalibrated or calibrated) spread of the ensemble rather than using a level of uncertainty which is constant with time. The uncalibrated data suggests that variations in uncertainty that are 20% to 55% of the mean uncertainty are predictable using the ensemble spread. However, because the uncalibrated data both underestimates the total spread (the numerator in the COVS) and overestimates the predictable part of the variability of the spread (the denominator in the COVS) these values seem to be overestimates. The calibrated data suggests that variations in the uncertainty that are only 5% to 20% of the mean uncertainty are predictable using the ensemble spread.

## 5 Other distributions

In cases where the forecast errors are not close to normally distributed, one can use other distributions. For example in the case where the forecast errors show skew the skew-normal distribution $`SN`$ can be used. The skew-normal distribution is a generalisation of the normal distribution which has a third parameter, and includes the possibility of modelling skew (Azzalini, 1985). Suppressing the index $`i`$ for clarity we then have:

$$TSN(\alpha +\beta m,\gamma +\delta s,\zeta +\eta k)$$
(3)
where we have introduced the ensemble skew $`k`$ and two new parameters $`\zeta `$ and $`\eta `$.

The skew-normal model can be fitted using maximum likelihood methods exactly as for the normal distribution. One of the results from such a fitting process would be a clear indication as to whether the forecast being calibrated does or does not contain statistically significant information about the skew of observed temperatures (this question has been discussed by Denholm-Price and Mylne (2002)).

For extremely non-normal distributions for which even the skew-normal is not non-normal enough, non-parametric distributions may be more appropriate. A simple non-parametric method would be to use a kernel density, with a single free parameter for the width of each kernel (see Bowman and Azzalini (1997) for a description of kernel densities). Such a method would look a little like the method of Roulston and Smith (2003) even though it is justified in a completely different way.

## 6 Conclusions

We have described a simple parametric method for the assessment and calibration of ensemble temperature forecasts. The method consists of applying linear transformations to the mean and standard deviations from the ensemble. The parameters of the model can be fitted easily using the maximum likelihood method. The model has various advantages and disadvantages relative to other calibration models currently in use. The main disadvantage is that the model only works for forecast errors that are reasonably close to normally distributed, although extensions have been described that should overcome that limitation. The advantages of the model are that:

* the calibration of forecasts using the model is extremely simple
* the model is transparent and easy to understand
* the model separates skill in predictions of the mean and the spread
* calculating confidence intervals on parameters is easy
* the model gives a clear indication of how many days of useful skill there are in a forecast

We have applied the model to one year of site-specific ECMWF ensemble forecasts. We find that the forecasts have highly significant skill for predicting both the mean and the standard deviation out to 10 days. The forecasts underestimate the mean uncertainty, as has been reported in other studies. They also over-estimate the variability of the uncertainty. For these forecasts we estimate that the predictable part of the uncertainty is only between 5% and 20% of the mean uncertainty, depending on lead time. For some applications this variability in the uncertainty may be small enough that it can be ignored and one could make the simplifying assumption that the uncertainty is constant in time.

Further work includes:

* Developing algorithms that avoid having to make the assumption that the forecast error is uncorrelated in time.
* Out of sample testing of the calibrated forecasts, using both measures from within the framework (i.e. likelihood) and also other measures such as rank histograms, reliability diagrams and ROCs.

## 7 Acknowledgements

We would like to thank K. Mylne for providing us with the forecast data, Risk Management Solutions for providing us with the observational data, and S. Mason, F. Doblas-Reyes and D. Anderson for useful discussions on the subject of ensemble validation. This research was funded by the authors.


text:
# Nuclear matter incompressibility coefficient in relativistic and nonrelativistic microscopic models

## I Introduction

The nuclear matter incompressibility coefficient $`K_{nm}`$ plays an important role in understanding a wide variety of phenomena ranging from giant resonances in finite nuclei to heavy-ion collisions and supernova explosions. The measurement of the centroid energy of the isoscalar giant monopole resonance (ISGMR) provides a very sensitive method to determine the value of $`K_{nm}`$. Over the last couple of decades, attempts have been made to measure very accurately the value of the ISGMR centroid energy, $`E_0`$. The recent experimental data Youngblood99  for the $`E_0`$ in heavy nuclei are accurate enough to provide unambiguous information about $`K_{nm}`$. However, the theoretical scenario of the $`K_{nm}`$ still remains unclear.

In the past, the value of $`K_{nm}`$ was determined using a macroscopic approach which relies on the liquid drop type of expansion for the breathing mode restoring force. It was shown that equally good fits can be obtained with values of $`K_{nm}`$ ranging from $`100400`$ MeV Pearson91 ; Shlomo93 . In other words, the liquid drop approach can not constrain the value of $`K_{nm}`$ better than $`50\%`$. On the other hand, the microscopic determination of $`K_{nm}`$ using the Hartree-Fock based random phase approximation (HF-RPA) has undergone a significant improvement over time. The earlier nonrelativistic HF$``$RPA calculations carried out using the Skyrme interaction, which nicely reproduced the gross properties of nuclei (such as nuclear binding energy and charge radii) and the available data on isovector giant dipole resonance and isoscalar giant quadrupole resonance, yielded values of about 370 MeV for the $`K_{nm}`$. With these interactions, the ISGMR in $`{}_{}{}^{208}\mathrm{Pb}`$ was predicted to be located at an excitation energy of about 18 MeV. The discovery of ISGMR in $`{}_{}{}^{208}\mathrm{Pb}`$ at an excitation energy of 13.7 MeV Marty75  led to the modification of the Skyrme interaction. Until today, the nonrelativistic calculations with Skyrme Colo92 ; Hamamoto97  and Gogny Blaizot95  interactions predict a value of $`K_{nm}`$ in the range of $`210220`$ MeV. We remark here that the long standing problem of the conflicting results deduced for $`K_{nm}`$ from data on the isoscalar giant dipole resonance (ISGDR) and data on the ISGMR was explained recently by Shlomo and Sanzhur Shlomo02  as being due to a missing strength in the experimental data for the high energy region of the ISGDR.

The relativistic mean field based RPA (also referred to as RRPA) calculations, with the contribution from negative-energy sea not included, yielded for $`K_{nm}`$ a value in the range of $`280350`$ MeV Ma97 . Recent RRPA calculations Ma01 ; Vretenar03 , with the inclusion of negative-energy states in the response function, yield a value of $`K_{nm}`$ between $`250270`$ MeV. The discrepancy of about $`20\%`$ in the value of $`K_{nm}`$ obtained from relativistic and nonrelativistic models is quite significant in view of the accuracy of the experimental data currently available on the ISGMR centroid energies. In recent studies Giai01 ; Niksic02  it has been claimed that these pronounced differences are due to the model dependence of $`K_{nm}`$. On the other hand, it has been pointed out in Ref. Piekarewicz02  that the differences in the values of $`K_{nm}`$ obtained in the relativistic and the nonrelativistic models can be attributed, at least in part, to the differences in the density dependence of the symmetry energy in these models. However, in Ref. Piekarewicz02 , the analysis was restricted to a single nucleus $`{}_{}{}^{208}Pb`$ and the interaction parameters for the several families of the effective Lagrangian considered were fitted only to the empirical values of saturation density, binding energy per nucleon in symmetric nuclear matter and the charge radius of the $`{}_{}{}^{208}Pb`$ nucleus. It may be pointed out here that in Ref. Serr81  a reasonable value of $`E_0`$ for the $`{}_{}{}^{208}\mathrm{Pb}`$ nucleus was obtained using an effective force with $`K_{nm}=400`$ MeV. But, the same effective force overestimated the value of $`E_0`$ in case of the $`{}_{}{}^{90}\mathrm{Zr}`$ nucleus. Moreover, it has been suggested in Ref. Serr81  that a wide range of combinations of bulk, surface and asymmetry contributions to the finite nucleus compressibility can fit the energy of the ISGMR in medium to heavy nuclei. This implies that for a meaningful informations about the discrepancy between the relativistic RRPA and the nonrelativistic HF-RPA calculations, one must compare the results obtained from these models for several nuclei.

It is noteworthy that in a crude approximation, the uncertainty of about $`20\%`$ in the values of $`K_{nm}`$ is tantamount to an uncertainty of $`10\%`$ in the value of $`E_0`$. This is because, in a semi-classical approach, $`E_0\sqrt{K_{nm}}`$. We have shown very recently Agrawal03  that the calculated value of $`E_0`$ can deviate by about $`5\%`$ if the particle-hole space is quite limited and/or self-consistency is not properly maintained. We note that in the published literature, the calculated values of $`E_0`$ for $`{}_{}{}^{208}\mathrm{Pb}`$, obtained for the same interaction, differs by up to 0.3 MeV Giai01 ; Colo99 ; Colo01 . Therefore, appropriate comparison for the values of $`K_{nm}`$ obtained from the different models is possible only when all the calculations are performed with the same procedure and numerical accuracy.

In this work we take a close look at the issue of the model dependence of the nuclear matter incompressibility coefficient derived from the ISGMR centroid energy. Toward this purpose, we generate different parameter sets for the Skyrme interaction and perform highly accurate calculations for the ISGMR strength function for several nuclei using the HF based RPA approach. For the sake of true comparison, the calculations using different parameter sets of Skyrme interaction are performed following exactly the same numerical procedure. The values of the Skyrme parameters are obtained by a least square fit to exactly the same experimental data for the nuclear binding energies, charge radii and neutron radii as adopted in Ref. Lalazissis97  in determining the NL3 parameter set for an effective Lagrangian used in the RMF model. We find that the model dependence is rather weak and the differences in the values of $`K_{nm}`$ in the relativistic and the nonrelativistic models essentially arise from the differences in the nuclear matter properties, in particular, in the values of the symmetry energy coefficient ($`J`$), associated with these models.

## II Formalism

In a self-consistent HF-RPA calculation Bertsch75 , one starts by adopting a specific effective nucleon-nucleon interaction, $`V_{12}`$. In this work we shall use the Skyrme type interaction of the form Vautherin72 ,

$`V_{12}`$ $`=`$ $`t_0\left(1+x_0P_{12}^\sigma \right)\delta (𝐫_1𝐫_2)+{\displaystyle \frac{1}{2}}t_1\left(1+x_1P_{12}^\sigma \right)`$ (1)
$`\times \left[\stackrel{}{k}_{12}^2\delta (𝐫_1𝐫_2)+\delta (𝐫_1𝐫_2)\stackrel{}{k}_{12}^2\right]+t_2\left(1+x_2P_{12}^\sigma \right)\stackrel{}{k}_{12}\delta (𝐫_1𝐫_2)\stackrel{}{k}_{12}`$
$`+{\displaystyle \frac{1}{6}}t_3\left(1+x_3P_{12}^\sigma \right)\rho ^\alpha \left({\displaystyle \frac{𝐫_1+𝐫_2}{2}}\right)\delta (𝐫_1𝐫_2)`$
$`+iW_0\stackrel{}{k}_{12}\delta (𝐫_1𝐫_2)(\stackrel{}{\sigma _1}+\stackrel{}{\sigma _2})\times \stackrel{}{k}_{12}`$
where, $`P_{12}^\sigma `$ is the spin exchange operator,$`\stackrel{}{\sigma }_i`$ is the Pauli spin operator, $`\stackrel{}{k}_{12}=i(\stackrel{}{}_1\stackrel{}{}_2)/2`$ and $`\stackrel{}{k}_{12}=i(\stackrel{}{}_1\stackrel{}{}_2)/2.`$ Here, right and left arrows indicate that the momentum operators act on the right and on the left , respectively. The parameters of the Skyrme force are obtained by fitting the HF results to the experimental data for the bulk properties of finite nuclei. Once the HF equations are solved using an appropriate parameter set for the Skyrme interaction, then one obtains the RPA Green’s function Bertsch75

$$G=G_0(1+V_{ph}G_0)^1,$$
(2)
where, $`V_{ph}`$ is the particle-hole ($`ph`$) interaction consistent with $`V_{12}`$ and $`G_0`$ is the free $`ph`$ Green’s function . For the single-particle operator

$$F=\underset{i=1}{\overset{A}{}}f(𝐫_i),$$
(3)
the strength function is given by

$$S(E)=\underset{n}{}\left|0\left|F\right|n\right|^2\delta \left(EE_n\right)=\frac{1}{\pi }Im\left[Tr\left(fGf\right)\right].$$
(4)
The steps involved in the relativistic mean field based RPA calculations are analogous to those described above for the nonrelativistic HF-RPA approach. However, the nucleon-nucleon interaction in case of relativistic mean field models are generated through the exchange of various mesons. An effective Lagrangian which represent a system of interacting nucleons looks like Lalazissis97 ,

$``$ $`=`$ $`\overline{\psi }\left(\gamma \left(ig_\omega \omega g_\rho \stackrel{}{\rho }\stackrel{}{\tau }eA\right)mg_\sigma \sigma \right)\psi +{\displaystyle \frac{1}{2}}\left(\sigma \right)^2`$ (5)
$`U(\sigma ){\displaystyle \frac{1}{4}}\mathrm{\Omega }_{\mu \nu }\mathrm{\Omega }^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\omega ^2\omega ^2{\displaystyle \frac{1}{4}}\stackrel{}{R}_{\mu \nu }\stackrel{}{R}^{\mu \nu }`$
$`+{\displaystyle \frac{1}{2}}m_\rho ^2\stackrel{}{\rho }^2{\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }`$
which contains nucleons $`\psi `$ with mass $`m`$; $`\sigma `$, $`\omega `$, $`\rho `$ mesons; the electromagnetic fields; and non-linear self-interactions of the $`\sigma `$ field,

$$U(\sigma )=\frac{1}{2}m_\sigma ^2\sigma ^2+\frac{1}{3}g_2\sigma ^3+\frac{1}{4}g_3\sigma ^4.$$
(6)
The Lagrangian parameters are usually obtained, as in the case of nonrelativistic mean field calculations, by fitting procedure to some bulk properties of a set of spherical nuclei Reinhard89 . The values of various coupling constants and the meson masses appearing in Eqs. (5) and (6) for the most widely used parameter set NL3 are $`m_\sigma =508.194`$ MeV, $`m_\omega =782.501`$ MeV, $`m_\rho =763.000`$ MeV, $`g_\sigma =10.217`$, $`g_\omega =12.868`$, $`g_\rho =4.474`$, $`g_2=10.431`$ fm<sup>-1</sup> and $`g_3=28.885`$.

## III Results and discussions

### III.1 HF results

In the present work, for an appropriate comparison, we carry out a least square fit to the same experimental data for the nuclear binding energies, charge radii and neutron radii used in Ref. Lalazissis97  to obtain the NL3 parameter set. Furthermore, we deal with the centre of mass correction to the total binding energy, finite size effects of the proton and the Coulomb energy in the way similar to that employed in determining the NL3 parameter set in Ref. Lalazissis97 . It may be pointed out that pairing is not included in our HF calculations and instead of the 10 nuclei considered in Ref. Lalazissis97 , we consider 7 closed shell nuclei. The open shell nuclei $`{}_{}{}^{58}\mathrm{Ni}`$, $`{}_{}{}^{124}\mathrm{Sn}`$ and $`{}_{}{}^{214}\mathrm{Pb}`$ are excluded from our least square fit. We also ignore the proton and the neutron pairing gaps in case of $`{}_{}{}^{90}\mathrm{Zr}`$ and $`{}_{}{}^{116}\mathrm{Sn}`$ nuclei, respectively. However, we have verified that if we increase the error bars for the experimental data on these two nuclei in order to compensate for the pairing, we find that the values of the Skyrme interaction parameters remains practically unaltered.

Since, the main objective of this paper is to delineate the differences in the value of $`K_{nm}`$ predicted by the relativistic and nonrelativistic mean field based RPA calculations, we generate a Skyrme interaction having $`K_{nm}`$ and $`J`$ very close to those associated with the NL3 parameter set, i.e., 271.76 and 37.4 MeV, respectively. Also, as most of the calculation claiming the model dependence are restricted to the ISGMR centroid energy for the single heavy nucleus $`{}_{}{}^{208}\mathrm{Pb}`$, we generate a parameter set by demanding a very high accuracy for the root mean square charge radius of $`{}_{}{}^{208}\mathrm{Pb}`$. We denote this parameter set as SK272. In addition, we also obtain a parameter set SK255 having characteristics very much similar to the SK272 parameter set, but, $`K_{nm}`$ is taken to be 255 MeV.

In Table 1 we give the values of the Skyrme parameters SGII together with the new parameter sets SK272 and SK255. In Table 2 we compare the nuclear matter properties for the SK272 and SK255 interactions with the corresponding ones obtained from the NL3 and SGII interactions. The quantities $`\rho _0`$, $`m^{}/m`$ and $`L`$ in this table denote the saturation density, effective nucleon mass and the slope of the symmetry energy coefficient ($`L=3\rho _0dJ/d\rho _0`$), respectively. In column 3 of Table 3 we have given the experimental data for the total binding energy, $`E`$, charge radii, $`r_c`$, and neutron radii, $`r_n`$, for the nuclei used in Ref. Lalazissis97  for determining the NL3 parameter set and in our fit, together with the assumed error bars in percent. The values obtained from the parameter sets SK272 and SK255 are shown in columns 5 and 6, respectively. For the sake of comparison, in Table 3 we also give in columns 4 and 7 the results for the NL3 and SGII interactions, respectively. It is evident from this table that the quality of our fit to the experimental data is quite comparable to the results obtained with the NL3 and SGII interactions. We have also listed the values of $`\mathrm{\Delta }r=r_nr_p`$, the difference between rms radii for neutrons and protons, (not included in the fit). Experimental values for $`\mathrm{\Delta }r`$ are obtained from the data for the neutron and charge radii and by using $`r_p=\sqrt{r_c^20.64}`$. It is interesting to note that the values of $`\mathrm{\Delta }r`$ obtained for the SK272, SK255 and NL3 interactions are closer and are quite large compared with the corresponding results for the SGII interaction. This is consistent with the fact that the values of the slope of the symmetry energy ($`L`$) associated with the SK272, SK255 and NL3 interactions are significantly larger than that associated with the SGII interaction (see Table 2). Since, as is well known Friedman77 , the value of $`\mathrm{\Delta }r`$ is sensitive to the density dependent form adopted for the symmetry interaction.

### III.2 RPA results

We have demonstrated very recently Agrawal03  that the strength function for giant resonances are quite sensitive to the various numerical approximations made. By numerical approximations we essentially mean, the size of the box used for the discretization of the continuum, restriction imposed on the maximum energy for the particle-hole excitations ($`E_{ph}^{max}`$) , smearing parameter ($`\mathrm{\Gamma }/2`$) used to smear the strength function, etc. We have shown in Ref. Agrawal03  that in order to reproduce the results obtained in the continuum RPA calculation the size of the box should be consistent with the value used for the smearing width. For $`\mathrm{\Gamma }/2=1`$ MeV, one must use a large box of size 72 fm. In this work, we have used a box of size 90 fm and $`\mathrm{\Gamma }/2=1`$ MeV. From Ref. Agrawal03  we also note that in order to obtain a reliable value of the ISGMR centroid energy, $`E_0`$, accurate within 0.1 MeV, one must use $`E_{ph}^{max}>400`$ MeV. Here, the centroid energy is given by $`E_0=m_1/m_0`$, where $`m_0`$ and $`m_1`$ are the non-energy-weighted and energy-weighted sums of $`S(E)`$ of Eq. (4), respectively. In the present work, the lowest value of $`E_{ph}^{max}`$ we have used is higher than 500 MeV. In addition, the centroid energy depends strongly on the range adopted for the excitation energy interval of the giant resonance. Nevertheless, one often encounters in the published literature that the values of the centroid energies are given without any reference to the corresponding excitation energy range considered. For instance, in the case of the $`{}_{}{}^{208}\mathrm{Pb}`$ nucleus we find for the SGII parameter set that $`E_0=`$ 13.7, 13.9, 14.4 MeV for the excitation energy ranges $`040`$, $`060`$ and $`1040`$ MeV, respectively. These differences are quite significant, since, as pointed out earlier, a variation of $`5\%`$ in $`E_0`$ corresponds to a change in $`K_{nm}`$ by $`10\%`$. In what follows, we shall concentrate mainly on the results for $`E_0`$ obtained by integrating the strength function over the energy range $`060`$ MeV, since the RMF results presented in Ref. Giai01  for the NL3 parameter set were obtained using the same energy range and the strength function was smeared using $`\mathrm{\Gamma }/2=1`$ MeV MaP .

In Table 4 we give the results for the ISGMR centroid energy obtained using the parameter sets SK272 and SK255 and compare them with the RRPA results of Ref. Giai01  for the NL3 interaction. It is clear from this table that the parameter set SK272 yields values for $`E_0`$ which are higher by about $`46\%`$ than those obtained using the NL3 parameter set in the RRPA calculations. This implies that if we reduce the compressibility by about $`10\%`$, we can reproduce reasonably well the RRPA results for the NL3 interaction. For this reason we have generated another parameter set SK255 with $`K_{nm}=255`$ MeV, keeping $`J=37.4`$ MeV (see Table 2). In fact, we see that the differences between the values of the $`E_0`$ obtained from the parameter sets SK255 and NL3 are on the level of the uncertainty associated with the experimental data for $`E_0`$. We emphasize that though the values of $`K_{nm}`$ associated with the SGII and SK255 parameter sets differ by about 40 MeV, the values of $`E_0`$ for $`{}_{}{}^{208}\mathrm{Pb}`$ nucleus for these interactions are close within 0.3 MeV. Thus, by requiring a fixed value of $`E_0`$ we find that an increase in $`J`$ by $`10\%`$ leads to an increase in $`K_{nm}`$ by about $`5\%`$. We also compute the values of $`E_0`$ over the same energy range as used in experimental determination of the centroid energy Youngblood99 . It should be noted that the experimental values of the energy range, given in the Table 4, are more or less symmetric around the corresponding $`E_0`$. It can be seen from Table 4, that for the parameter set SK255, we obtain a good agreement with the experimental data for $`E_0`$, calculated over the experimental excitation energy range. We remark that, in our calculations for the $`{}_{}{}^{208}\mathrm{Pb}`$ nucleus with the SK272 and SK255 parameter sets, the peak energy for the isovector giant dipole resonance is 13.2 and 13.3 MeV respectively, which is in good agreement with the experimental value of $`13.3\pm 0.1`$ MeV Ritman93 .

### III.3 Conclusions

In summary, we have analyzed in detail the recent claim that the nuclear matter incompressibility coefficient $`K_{nm}`$ extracted from the ISGMR centroid energy calculated using the relativistic and nonrelativistic based RPA models differ by about $`20\%`$. For a meaningful comparison, we have generated parameter sets for the Skyrme interaction by a least square fitting procedure using exactly the same experimental data for the bulk properties of nuclei considered in Ref. Lalazissis97  for determining the NL3 parameterization of an effective Lagrangian used in the relativistic mean field models. Further, we also demanded in our fitting procedure that the values of $`K_{nm}`$, $`J`$ and the charge radius of the $`{}_{}{}^{208}\mathrm{Pb}`$ nucleus should be very close to the results obtained with the NL3 interaction. The parameter sets thus obtained were used to calculate the ISGMR centroid energy for several nuclei. For the parameter set SK272 ($`K_{nm}=272`$ MeV), the calculated values of $`E_0`$ are higher by about $`5\%`$ compared to the corresponding NL3 results. This implies that the difference in the value of $`K_{nm}`$ obtained in the relativistic and the nonrelativistic microscopic models could be at most $`10\%`$. In view of this, we generate another parameter set having $`K_{nm}=255`$ MeV. As expected, the parameter set associated with $`K_{nm}=255`$ MeV, yields for the ISGMR centroid energies values which are quite close to the NL3 results. Moreover, for the SK255 parameter set, we find a good agreement with experimental data for $`E_0`$ for all the nuclei considered, provided, the corresponding excitation energy ranges used in determining $`E_0`$ are the same as those used in obtaining the experimental data. We have thus shown that the difference in the values of $`K_{nm}`$ obtained in the relativistic and nonrelativistic models is mainly due to the differences in the values of the symmetry energy coefficient ($`J`$) and its slope ($`L`$) associated with these models.