text:
# Building up the globular cluster system of the Milky Way. The contribution of the Sagittarius galaxy

## 1 Introduction

The quest to find ordered structures in the system of satellites of the Milky Way galaxy dates back to 30 years ago when possible alignments of globular clusters and/or dwarf galaxies along wide streams were first noted (Hodge & Michie, 1969; Lynden-Bell, 1976; Kunkel & Demers, 1977; Kunkel, 1979; Lynden-Bell, 1982). The mounting consensus for scenarios in which the accretion of satellites has a major role in the formation of the outer halo of the Galaxy (Searle & Zinn, 1978; Zinn, 1993) prompted a new burst of such kind of studies since the mid ’90s up to the present day (Majewski, 1994; Lynden-Bell & Lynden-Bell, 1995; Fusi Pecci et al., 1995; Dinescu, Girard & van Altena, 1999; Palma, Majewski & Johnston, 2002). Despite the many interesting suggestions, none of the quoted studies was able to provide a conclusive proof of the reality of the alignments, mainly because of the overwhelming difficulty to assess the statistical significance of structures formed by inherently small numbers of objects.

However, recent theoretical and observational achievements may help us to look into the problem from a different and more fruitful perspective:

1. N-body simulations of the process of Galaxy assembly, starting from standard cosmological conditions (Cold Dark Matter - CDM), strongly suggest that the hierarchical merging of satellites is the main driver of galaxy formation (see Moore et al., 1999; Klypin et al., 1999, and references therein).
2. Independent of the assumed cosmology, it has been demonstrated that the accretion/disruption of satellites into the halo of a larger galaxy may leave long-lived relics in the form of streams of stellar (and/or dark matter) remnants that remain aligned to the orbital path of the parent satellite (Johnston, Hernquist & Bolte, 1996; Johnston et al., 1999; Ibata et al., 2001a; Ibata & Lewis, 1998; Mayer et al., 2001).
3. Convincing observational evidence of the clumpy and “filamentary” nature of the Galactic halo have been provided by many different groups (e.g Newberg et al., 2002; Vivas et al., 2001; Dohm-Palmer et al., 2001; Yanny et al., 2000; Ivezic et al., 2000; Ibata et al., 2001a; Helmi et al., 1999; Helmi, 2001; Majewski et al., 1999).

The in vivo example of a satellite accretion/disruption is provided by the Sagittarius dwarf Spheroidal galaxy (Sgr dSph; Ibata et al., 1994, 1997), which is currently merging with the Milky Way, and is carrying its own globular cluster system (i.e., M 54, Ter 8, Arp 2 and Ter 7, previously believed to be normal Galactic globulars). There is now clear observational evidence that the Sgr dSph is loosing stars under the strain of the Milky Way tidal field. These tidally-removed stars are found along a huge (and quite coherent) stream extending all over the sky (Sgr Stream, see Ibata et al., 2001a; Dohm-Palmer et al., 2001; Martínez-Delgado et al., 2001; Newberg et al., 2002, and references therein), tracing the orbit of the parent galaxy.

Ibata & Lewis (1998, hereafter IL98), and Ibata et al. (2001a) have simulated the evolution of the Sgr dSph over several orbital periods (P $`1`$ Gyr), computing the orbit of the galaxy as well as the phase-space distribution of the debris under different assumptions about the flattening of the CDM halo. The initial conditions of the simulations were based on the known position and radial velocity of Sgr dSph and on its proper motion as estimated by Ibata et al. (1997, 2001a). The orbit has a planar rosette structure, with the pole of the orbit located at \[$`\mathrm{}=90^{}`$, $`b=13^{}`$\] (i.e. a nearly polar orbit), and peri- and apo-Galactic distances of $`15\mathrm{kpc}`$ and $`60\mathrm{kpc}`$ respectively. The derived orbit has been successfully compared with the observed position of the Sgr Stream (Ibata et al., 2001a, b), providing also remarkable indications that the dark halo of the Milky Way is nearly spherical.

In this framework it is a tantalizing application to look for other halo globulars that may be correlated with the orbital path of the Sgr dwarf, and which could be lying in the Sgr Stream. In particular, we look for the phase-space coincidence of outer halo globulars with the computed orbit of the Sgr dSph from 1 Gyr ago up to the present day, searching for the most recent episodes of globular cluster loss, i.e. the ones whose traces are most likely to be still detectable.

## 2 Looking for structures

For our comparison we selected from the catalogue by Harris (1996) the 35 globular clusters in the range of galactocentric distance $`10\mathrm{kpc}R_{GC}40\mathrm{kpc}`$. Among these, 33 have also measured radial velocity $`V_r`$. For sake of brevity and clarity we will call this sample the Outer Halo Sample (OHS), in the following. With this selection we avoid the central part of the Galactic halo where it is less likely that ordered structures can survive for a long time, and we leave out of the sample the handful of clusters lying outside of $`R_{GC}60\mathrm{kpc}`$, a region that lies beyond the Sgr Stream according to the IL98 orbit111There is a 28 kpc wide gap in the radial distribution of Galactic globular clusters (see, e.g. Zinn, 1985, and references therein). There are only five clusters beyond $`R_{GC}=40`$ kpc, namely Pal 14, Eridanus, Pal 3, NGC 2419, Pal 4, and AM1. Their respective galactocentric distances are $`R_{GC}=`$ 65, 83, 90, 98, 99, and 117 kpc, much beyond the apogalacticon of the Sgr orbit. The adopted outer radial threshold ($`R40\mathrm{kpc}`$, quite similar to the one adopted by Palma, Majewski & Johnston, 2002, i.e., $`R_{GC}36\mathrm{kpc}`$) provides the selection of a homogeneous sample without significant gaps.. The adopted OHS does not include the known Sgr globulars, to avoid the detection of the obvious signal of their clustering around the center of the Sgr galaxy.

In Figure 1 we show the OHS clusters (small solid circles) and the Sgr orbit in the planes formed by the rectangular Galactocentric coordinates222This Galactocentric coordinate system is defined such that the origin lies at the Galactic Center; at the Solar position, $`(8,0,0)`$, the $`Y`$-axis points in the direction of Galactic rotation; while the $`Z`$-axis towards the North Galactic Pole. ($`X,Y,Z`$, in kpc) and in the $`R_{GC}`$ \[kpc\] vs. $`V_r`$ \[km/s\] plane. The large full circles are the known Sgr globulars, which we also show in the plots for completeness. Note that these clusters lie around the end of the orbit corresponding to the present time ($`t=0`$). We highlight (with encircled solid circles) six more clusters that lie remarkably close to the orbit in all the considered planes. These clusters are: Pal 12 (whose association to the Sgr Stream has been already established by Dinescu et al., 2000; Martínez-Delgado et al., 2002), NGC 4147, NGC 5634, NGC 5053, Pal 5 and Ter 3. Is this association real or could it be the mere occurrence of a chance alignment? Though chance alignments in the four-dimensional phase space (X,Y,Z,$`V_r`$) are not expected to be very likely, the key point is to quantify the probability that the observed structure could have originated from a statistical fluctuation. To do this we will compare the observed distribution - and its phase space distance to the Sgr orbit - with synthetic samples (having the same dimension as the OHS) extracted from a model representing an unstructured parent halo.

The most conservative comparison that can be made is with a model that closely resembles the observed radial and velocity distribution of the OHS. Figure 2 (upper panel) shows that the cumulative radial distribution of the OHS is well reproduced by a spherical halo model with a density distribution $`R^{1.6}`$. A Kolmogorov-Smirnov (KS) test shows that the probability that the OHS is drawn from the $`\mathrm{\Phi }R^{1.6}`$ model is $`90`$%. On the other hand, the probability that the same sample is drawn from the other two models shown for comparison ($`\mathrm{\Phi }R^{1.0}`$, and $`\mathrm{\Phi }R^{2.5}`$) is $`15`$%. Doubts may be cast on the appropriateness of a spherical model. It may be conceived that if the true parent halo is flattened, some excess of clustering of the observed points along an orbit with low inclination may artificially emerge in the comparison with a spherical model. This is clearly not the case, however, since the IL98 orbit is nearly polar, i.e. it is almost perpendicular to the Galactic Plane (see Figure 1). In the lower panel of Figure 2 it is shown that the observed distribution of radial velocity of the OHS is well reproduced by a Gaussian distribution with $`<V_r>=38\mathrm{km}\mathrm{s}^1`$ and $`\sigma _V=175\mathrm{km}\mathrm{s}^1`$. According to a KS test the probability that the observed sample is drawn from the model distribution is $`90`$%.

In the following simulations we extract all the synthetic samples from a spherical and isotropic model with $`\mathrm{\Phi }(R_{GC})R_{GC}^{1.6}`$ and with the Gaussian distribution of radial velocities shown in Figure 2. For each simulated cluster (as well as for all the OHS ones) we computed the spatial distance from the nearest point in the Sgr orbit ($`D_{orb}`$, in kpc) and the difference between their radial velocity and the one predicted from the computed orbit at that point ($`\mathrm{\Delta }V_r=V_r(obs)V_r(orb)`$).

In Figure 3, the $`D_{orb}`$ values of the selected clusters (large filled circles) are plotted against their $`\mathrm{\Delta }V_r`$. The encircled points are the six clusters highlighted in Figure 1. A sample of 10000 synthetic clusters (dots) extracted from the adopted model is also shown, for comparison, in the upper panel of Figure 3. The OHS clusters show a remarkable over-density toward the Sgr orbit, that lies in the origin of the axis in the considered plane. The dashed dotted lines enclose the points whose observed radial velocity is within $`\pm 60\mathrm{km}\mathrm{s}^1`$ of the velocity predicted by the IL98 orbit. Note that the expected velocity dispersion of the Sgr debris along the Sgr Stream is $`\sigma 60\mathrm{km}\mathrm{s}^1`$, according to Ibata et al. (2001b). The continuous vertical segments are placed at $`D_{orb}=`$6, 12, and 18 kpc.

The lower panel of Figure 3 is arranged in the same way, but in this case the dots represent the distribution of the points in the best-fitting Ibata et al. (2001a) simulation that retains a bound core to the present day. The boldface dots are the particles that remain bound to Sgr or were bound less than 3 Gyr ago, while the ordinary dots are particles that were already unbound at that time. Note the remarkable similarity with the distribution of the OHS clusters. In particular, the correlation between the particles that have flown in the Sgr Stream in recent times (less than 3 Gyr ago) and the OHS clusters with $`|\mathrm{\Delta }V_r|<60\mathrm{km}\mathrm{s}^1`$ at any $`D_{orb}`$ is quite striking.

We define $`N_{D<x}`$ as the number of synthetic or real clusters having $`60\mathrm{km}\mathrm{s}^1<\mathrm{\Delta }V_r<+60\mathrm{km}\mathrm{s}^1`$ and $`D_{orb}<x\mathrm{kpc}`$. For example, it can be seen from Figure 3 that there are five OHS clusters with $`|\mathrm{\Delta }V_r|<60\mathrm{km}\mathrm{s}^1`$ and $`D_{orb}<6\mathrm{kpc}`$, thus $`N_{D<6}=5`$. In the same way we find $`N_{D<12}=10`$ and $`N_{D<18}=14`$ from the observed sample. Now the question is what is the probability that a sample having $`N_{D<x}`$ greater or equal to the observed one ($`N_{D<x}^{obs}`$) is drawn from the assumed unstructured model? To answer this question we randomly extracted 10000 samples of 32 synthetic clusters from the assumed model, and for each of them we measured $`N_{D<6}`$, $`N_{D<12}`$, and $`N_{D<18}`$.

Figure 4 reports the distributions of $`N_{D<6}`$ (upper panel), $`N_{D<12}`$ (middle panel), and $`N_{D<18}`$ (lower panel) for the 10000 simulated samples. The respective observed values are indicated by a dashed line. In all of the considered cases the mean and modal $`N_{D<x}`$ values of the distribution are much smaller than the observed ones. $`N_{D<6}5`$ occurs for only 200 simulated samples in 10000 (2% of the cases), $`N_{D<12}10`$ occurs only for 42 samples in 10000 ($`0.4`$%), and $`N_{D<18}14`$ occurs only in 2 cases in 10000 ($`0.02`$%). It can be concluded that it is highly improbable that the observed clustering of the OHS globulars around the Sgr orbit occur by chance alignment. It is important to remark again that the model from which the simulated samples have been extracted has the same $`R_{GC}`$ and $`V_r`$ distribution of the observed sample. Thus the Sgr orbit appears to be a strongly preferred subset of the phase-space for the Galactic globulars in the considered range of $`R_{GC}`$.

Finally, we note that the simulated samples that have a significant probability of realization (e.g., $`P20`$%) have $`N_{D<6}3`$, $`N_{D<12}6`$, and $`N_{D<18}8`$, significantly lower than the observed values. Hence, the significance of the detected phase-space structure cannot be due to the actual correlation of just a couple of clusters with the Sgr orbit. According to the distributions shown in Fig. 4, the probability that $`N_{D<18}10`$ is $`P5`$%. Therefore it can be stated, with 95% confidence, that at least 4 real associations are needed to produce the observed signal of $`N_{D<18}=14`$.

The above test demonstrates that the observed phase-space clumping around the orbit of the Sgr dwarf galaxy is highly unlikely to have occurred by a chance coincidence, if the halo globular cluster population is indeed distributed according to the simple spherical model described above. To address this concern, we performed a second test, using “bootstrapped” artificial data sets constructed from the real OHS. The artificial samples were constructed by rotating the position of each real globular cluster by a random azimuthal angle about the Galactic center, and introducing a random flip perpendicular to the Galactic plane. This implicitly assumes, of course, that the globular cluster distribution is symmetric about the Galactic center, as well as above and below the Galactic plane. A slight complication arises from the fact that we do not know the proper motion of many of the clusters, so we cannot deduce the heliocentric radial velocity that the cluster should have at its new position. We therefore have to assume some model for the halo velocity distribution; we take a simple Gaussian model, and draw random realizations of the total space velocity $`\stackrel{}{V}`$ consistent with the observed heliocentric velocity $`V_r`$ (that is, we calculate the conditional probability of $`\stackrel{}{V}`$ given $`V_r`$). Having defined thereby the 3-dimensional velocity vector, we rotate the vector to the new position, and project it along the line of sight to the Sun (corrections for the peculiar motion of the Sun and for Galactic rotation are made). With a Gaussian model that has $`\sigma =175\mathrm{km}\mathrm{s}^1`$ and $`V_r=38\mathrm{km}\mathrm{s}^1`$, we find that out of 10000 artificial data sets, $`F(N_{D<6}5)=399`$, $`F(N_{d<12}10)=102`$, and $`F(N_{D<18}14)=18`$. For comparison, for a model with $`\sigma =110\mathrm{km}\mathrm{s}^1`$ (or $`\sigma =150\mathrm{km}\mathrm{s}^1`$) and $`V_r=0\mathrm{km}\mathrm{s}^1`$, the statistics are as follows: $`F(N_{D<6}5)=463(425)`$, $`F(N_{D<12}10)=140(121)`$, and $`F(N_{D<18}14)=35(28)`$, where the result in brackets refers to the higher velocity dispersion model. According to these tests using “bootstrapped” random data sets, the observed phase-space clumping is highly unlikely to have occurred by chance, in agreement with our previous test. The fact that the observed structure has a slightly lower statistical significance with the “bootstrapped” samples than was deduced from the previous test with its perfectly isotropic and isothermal random samples is not surprising. This is presumably a consequence of occasional chance reappearances of (part of) the structure of the real OHS in the “bootstrapped” samples.

As a final remark, we note that the region of the $`\mathrm{\Delta }V_rD_{orb}`$ plane with $`|\mathrm{\Delta }V_r|<60\mathrm{km}\mathrm{s}^1`$ and $`D_{orb}<18`$ kpc is also populated by particles that become unbound more than 3 Gyr ago. Thus it has to be considered the possibility that some of the OHS clusters found in that region drifted into the Sgr Stream in more ancient times. For instance, this can be the case for NGC 5634, as we argue elsewhere (Bellazzini, Ferraro & Ibata, 2002). On the other hand, the discussion of the individual OHS clusters is beyond the scope of the present analysis.

### 2.1 Proper Motions

Thanks to the painstaking effort of a few teams of astronomers, estimates of the proper motions (hereafter PMs) of 41 Galactic globular clusters are now available (see Dinescu, Girard & van Altena, 1999; Palma, Majewski & Johnston, 2002, for complete lists and references). Most of these measures have quite sizeable uncertainties, ranging from $`0.2`$ mas/yr to $`2`$ mas year. Moreover, subtle - and unaccounted - systematics may still affect them, as suggested by the large differences between independent estimates of the motion of the same cluster. In the list by Palma, Majewski & Johnston (2002) there are 14 cluster with more than one PM estimate. The average absolute difference among independent estimates for the same cluster are $`1.72\pm 1.29`$ mas/year in $`\mu _\alpha cos\delta `$ and $`2.90\pm 2.50`$ mas/year in the $`\mu _\delta `$ component. Differences much larger than the quoted errors (up to $`10`$ times) are common and a reasonable assumption of the minimum uncertainty of the proper motion estimates is probably $`1`$ mas/year in each component. Hence, while available PMs may be very valuable, for instance, to characterize the kind of orbit followed by a given globular, they are expected to have a modest constraining power in the application presented in this paper. Such constraining power is also lowered by the fact that PMs predicted by the model we are comparing with carry their own, non negligible, uncertainty (typically $`0.85`$ mas/yr in each component, which we estimate from the intrinsic dispersions in distance, position and velocity in the N-body stream models, assuming a 20% uncertainty in the Galactic circular velocity at $`50\mathrm{kpc}`$). Radial velocities, being unaffected by uncertainties in the distance, are more reliably predicted by the model and are easily, and much more accurately, measured. For these reasons we preferred to limit the bulk of our analysis to the the (X,Y,Z,$`V_r`$) phase-space.

Yet it may be interesting to compare the PMs predicted by the IL98 model for the clusters that are candidate members of the Sgr Stream (listed in Table 1, see below), with the observational estimates. The orientations (directions) of the PM vectors are expected to be less sensitive to systematics with respect to the actual moduli. In Fig. 5 we compare the predicted and observed directions of the PM vectors of the seven clusters listed in Tab. 1 for which PM estimates are available. Four of these seven clusters have two independent PM estimates: in these cases we report the comparison with both values (see caption). We stress that this kind of comparison cannot have a serious impact on our previous conlusions except if both the following conditions were simultaneously fulfilled: (a) if accurate PM estimates were available for all the considered clusters (or the large majority of them), and (b) if less than 3-4 clusters were found to have observed PMs in agreement with the predictions. This case would imply that the correlation among OHS clusters and the Sgr orbit we detect in the (X,Y,Z,$`V_r`$) phase-space is largely due to a very unlikely (but not impossible) chance alignment in four dimensions. On the other hand some (up to $`10`$, actually) of the clusters listed in Table 1 are expected not to belong to the detected structure.

In panel (a) the predicted proper motions include a correction for the Solar reflex motion (hereafter S.r.m.), assuming a Galactocentric distance of $`8\mathrm{kpc}`$, a circular velocity of $`220\mathrm{km}\mathrm{s}^1`$ and the peculiar motion of the Sun as derived by Dehnen & Binney (1998). The predicted PM directions of all the seven clusters agree with the observed values to within $`\pm 45^{}`$. The adoption of alternative PM estimates does not change significantly the observed correlation. However, the result of this comparison may be affected by the inclusion of the Solar reflex motion. Hence, in panel (b) we present the comparison between the plain (uncorrected for the S.r.m) model predictions and the observed PM directions after correction for the S.r.m. The result is very similar to that presented in panel (a) despite of the larger errors involved. The only cluster whose predicted and observed PM directions differ by more than $`45^{}`$ is Pal 12, e.g. the only one for which independent evidences of membership to the Sgr Stream are already available (Martínez-Delgado et al., 2002).

In conclusion, it is found that the available proper motion estimates are fully compatible with the results presented in §2. Given the high degree of correlation among predicted and observed PM directions shown in both panels of Fig. 5 one may also be tempted to draw firm conclusion on the actual membership of individual clusters. Nevertheless, in our view, this may still be an hazard. If one looks at the amplitude of the individual corrections it is found that (a) all the considered clusters have at least one component of the PM vector that is corrected by more than 50 % of the observed value, (b) 4 out of seven clusters have at least one component whose S.r.m correction is larger that 100 % of the observed values and, (c) corrections as large as 200 - 400 % are indeed applied in some cases. Given the remarkable uncertainties in the observed PMs and the uncertainties involved in the S.r.m correction it is easy to conclude that currently measured PMs are not good tracers of the actual motion of the considered clusters and are instead predominantly measures of the Solar reflex motion. These considerations fully support the approach followed in §2, e.g. using only positions and radial velocity measurements in the search for statistical structures in the globular cluster system of the Milky Way.

### 2.2 Comparison with previous analysis

The possible association of Galactic globulars with the Sgr dSph galaxy has been considered before (Irwin, 1999; Dinescu, Girard & van Altena, 1999; Dinescu et al., 2000; Siegel at al., 2001; Palma, Majewski & Johnston, 2002). All these studies strongly rely on PM estimates. Furthermore they are centered on testing the actual association of individual clusters to the Sgr galaxy and its remnants. On the other hand our approach and our results are statistical in nature. We find a statistically significant clustering of OHS clusters along the Sgr orbit but we have no possibility to check individual membership, except for the limited test shown in §2.1. We can only provide a list of possible members, ranked according to their distance to the Sgr orbit in the (X,Y,Z,$`V_r`$) phase-space (e.g., Tab. 1). Hence, the comparison with the above quoted studies can be made only in a broadly general sense.

Dinescu, Girard & van Altena (1999) performed a thorough study of the orbits of the Galactic globulars for which PM estimates are available. In doing this, they identify three metal-poor clusters having orbits typical of the thick disc (instead that of the halo, as expected) and they check the hypothesis that the peculiar orbits of these clusters (namely NGC 6254, NGC 6626 and NGC 6752) could be due to their former membership to the Sgr dSph, concluding that this possibility is unlikely. All these clusters have $`R_{GC}<10`$, and are thus are not included in the OHS considered here and we cannot comment further on the Dinescu, Girard & van Altena (1999) results.

Following the suggestion of Irwin (1999), Dinescu et al. (2000) modeled the (PM based) orbits of the Sgr dSph and of Pal 12 and concluded that the available observations are compatible with the possibility that the cluster was torn from the galaxy $`1.41.7`$ Gyr ago. Our analysis confirms this result.

Siegel at al. (2001) obtained an estimate of the PM of Pal 13 and compared the general characteristics of the integrated orbit of the cluster (total energy, angular momentum, eccentricity, perigalactic and apogalactic radii) to those of five satellite galaxies of the Milky Way, including Sgr. They conclude that a common origin for Pal 13 and Sgr is unlikely. However we note that (a) Sgr provides (by far) the best match to the orbital parameters of Pal 13 with respect to any other considered galaxy (see Siegel et al.’s Tab. 8), the maximum discrepancy being $`40`$% in angular momentum, and (b) the orbital parameters are provided without any uncertainty, as if the adopted PMs were perfectely known. In our view this shortcoming seriously weakens the conclusion by these authors. Here we adopt the same PM estimate obtained by Siegel at al. (2001), reaching the conclusion that Pal 13 is a possible member of the Sgr Stream by direct comparison with a model that includes the Sgr Stream itself.

Palma, Majewski & Johnston (2002) (hereafter PMJ02) performed a global search for structures made of Galactic satellites, adopting the technique of the “orbital poles”, originally introduced by Lynden-Bell & Lynden-Bell (1995). In this kind of analysis one searches for intersections among the great circles described by the familiy of possible poles of the orbit of any given satellite. As stated by PMJ02 “any set of objects having common origin and mantaining a common orbit, no matter how spread out in the sky, will have great circle pole families (GCPFs) that intersects at the same pair of antipodal points on the celestial sphere”. Adding the information from PMs, PMJ02 limit the valid pole family of each satellite to an arc of its great circle of orbital poles (arc segment pole family, ASPF). With this technique they select a set of candidate-former-members of the Sgr galaxy (Ter 7, Ter 8, Arp 2, M 54, M 53, NGC 5053, Pal 5, M 5, Pal 12 and NGC 6356) by requiring that their GCPF come within 5 deg from the Sgr ACPF and $`6R_{GC}36`$ Kpc. The first four clusters in their list are the well known present members of the galaxy. They rise doubts of the clear association of each of the other candidate members. For instance, they exclude NGC 5053 and M 53 because their metal content is lower than any of the known Sgr globulars333This is not completely true since M 53 has the same metallicity as Ter 8, i.e. $`[Fe/H]=1.99`$ (see Montegriffo et al., 1998; Da Costa & Armandroff, 1995). See also below, for further discussion about non-kynematic criteria of selection of candidate Sgr members.. With the same arguments adopted by Siegel at al. (2001) for Pal 13 they argue that Pal 5 and Pal 12 are unlikely members of the family, but they cannot rule out the possibility. Finally, they conclude that despite the counter-arguments presented they mantain their original list of candidates, with the exception of NGC 5053. Elsewhere in the paper, they also note that NGC 5466 has energy and angular momentum within 1-$`\sigma `$ of those of the Sgr galaxy. NGC 6356, M 5 and M 53 have $`R_{GC}<10`$ kpc, hence are not included in our sample. The possible association of NGC 5466 with the Sgr dSph is confirmed by our analysis. The possible association of NGC 5053, Pal 5 and Pal 12 is also confirmed by the present work. PMJ02 also suggested a possible association of NGC 4147 with a group of clusters whose orbits may be related to those of the Magellanic Clouds, but they state that the result is very uncertain. Furthermore we note that on the same basis the possible association of M 53 (elsewhere classified as a candidate for association with Sgr) with the Small Magellanic Cloud and Ursa Minor is also suggested by PMJ02. We conclude that there is no serious disagreement between the present analysis and the one by PMJ02.

The present contribution, however, provides the first proof that the orbit of Sgr is a preferential subset of phase-space for globular clusters inhabiting the outer halo of the Milky Way, and we place this result on a sound statistical basis, a result not accomplished by any of the previous studies. The success of the present analysis is due to the direct comparison of the positional and kinematical properties of the OHS clusters with a realistic model of the orbit of the Sgr dSph and its relics, that has been previously tested against independent observations (Ibata et al., 2001a, b).

Finally, we shall shortly comment on “non-kinematic” criteria that have been used to assess the association of a given cluster to the Sgr galaxy. Many authors (see, e.g. Dinescu, Girard & van Altena, 1999; Palma, Majewski & Johnston, 2002, and references therein) have discussed the likeliness of the membership of their candidates on the basis of the similarity of their metallicity, Horizontal Branch (HB) morphology and/or structural parameters with those of the four known Sgr clusters. In our opinion this approach may not be very useful. It is quite clear that Sgr was a much larger and complex system in the past and its present status (as well as its present GC system) may be not fully representative of its original range of properties. If, by chance, the cluster Ter 7 had been lost by the Sgr galaxy in the previous perigalactic passage, this kind of criteria would have erroneously classified it as a bad candidate for Sgr membership, since Ter 7 is much more metal rich ($`[Fe/H]0.5`$) than any of the other clusters that are still present in the main body of Sgr dSph ($`[Fe/H]<1.5`$). Moreover, the existence of population and/or metallicity gradients may have driven the preferential loss of metal poor populations (see Alard, 2001; Newberg et al., 2002; Bellazzini, Ferraro & Ibata, 2002, and references therein), making the actual low metallicity limit of the Sgr stellar population quite uncertain. The judgement of the likelihood of membership on the basis of the HB morphology is even less justified since it is well known that clusters of different morphologies do actually co-exist in real galaxies (see, e.g., Buonanno et al., 1998, and references therein, for the case for Fornax).

A much more reliable discriminant may be provided, in the future, by the detailed abundance patterns of the clusters. For instance, the behaviour of $`\alpha `$ elements as a function of $`[Fe/H]`$ is determined by the star formation history of the parent galaxy (see McWilliam, 1997, and references therein) that was (probably) not the same in the Sgr dSph and in the Galactic Halo. In this context it is interesting to note that Pal 12, Ru 106 (Brown, Wallerstein & Zucker, ) and Pal 5 (Smith, Sneden & Kraft, 2002) have been found to be less $`\alpha `$-enhanced than Galactic Halo globulars of similar metallicity.

For the above reasons, we rely only on phase-space parameters in our analisys, recalling that by searching the lost relics of the Sgr system we are trying the reconstruct its original properties, not the present day ones.

## 3 Conclusions

We have demonstrated that there is a coherent structure in the phase-space distribution of the outer halo globular cluster of the Milky Way, strongly correlated with the orbital path of the Sgr dSph galaxy. This correlation cannot have originated by chance, as a random realization of an unstructured halo, if such halo has the same distribution of galactocentric distance and radial velocity as the considered OHS clusters.

Several of the OHS globular clusters have spatial positions and radial velocities compatible with the hypothesis that they are following the same orbit as the Sgr dSph, i.e. they probably belong to the Sgr Stream. Furthermore the spread in phase space around the Sgr orbit of such clusters is similar to that predicted for the Sgr Stream population, according to the numerical simulations by Ibata et al. (2001a).

We conclude, with 95% confidence, that at least 4 (out of 32) of the OHS clusters are physically associated with the Sgr Stream and were former members of the Sgr dSph. It should be noted however, that the analysis presented in this contribution is best suited to identify clusters that became unbound relatively recently. In principle, more OHS globular clusters could be associated with the Sgr stream.

In Table 1 we report the list of the selected globulars having $`|\mathrm{\Delta }V_r|<60`$ Km/s and $`D_{orb}<18\mathrm{kpc}`$, in order of increasing $`D_{orb}`$. The proper motions predicted by the IL98 model are also reported in Table 1. Since the adopted analysis is statistical in nature it is clear that some of the clusters listed in Tab. 1 may be just occasional interlopers of the Sgr orbit (up to 10 over 14 in the worst case, according to our estimate). The comparison with the available proper motion estimates shown in §2.1 does not provide any strong indication in this regard. While the actual membership to the Sgr Stream of each individual cluster can be firmly established only with very accurate proper motion measures (as will be achieved by dedicated space missions like GAIA or SIM), it is reasonable to assume that the clusters with the smaller $`D_{orb}`$ are the ones for which the membership is more likely (see §2.1). We make the hypothesis that the six clusters more closely associated with the Sgr orbital path (i.e. those put in evidence in Figures 1 and 3) are former members of the Sgr dSph. In this case the number of globular clusters originally in the Sgr galaxy is $`N_{cl}=10`$. If we adopt for the Sgr dSph the same globular cluster specific frequency ($`S_N=N_{cl}10^{0.4(M_V+15)}`$ Harris & Van den Bergh, 1981) of the only other Galactic dSph that has a globular cluster system, i.e. Fornax $`S_N26`$, we obtain an estimate of the total absolute magnitude of Sgr before the occurrence of any tidal stripping, $`M_V14`$. Newberg et al. (2002) estimated that in the Sgr Stream there are as many stars as in the present undisrupted body of the galaxy, thus the ab initio total luminosity of the Sgr dSph was roughly two times the present value. Hence the total absolute magnitude at that time was $`M_V\mathrm{log}(2)14.1`$ (where $`M_V=13.4`$, from Mateo, 1998). The excellent agreement between the two independent estimates of the initial $`M_V`$ fully supports the plausibility of the proposed scenario.

According to the results presented, it emerges that the Sgr dSph was not only an important contributor of halo field stars (Ivezic et al., 2000; Yanny et al., 2000; Newberg et al., 2002) but it also had a significant role in the building-up of the globular clusters system of the Milky Way.

M.B. and F.R.F. acknowledge the financial support to this research by the italian Ministero dell’Universitá e della Ricerca Scientifica (MURST) through the grant p. 2001028879, assigned to the project Origin and Evolution of Stellar Populations in the Galactic Spheroid. This research has made use of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# Electromagnetic Wave Propagation in Media with Indefinite Permittivity and Permeability Tensors

## Abstract

We study the behavior of wave propagation in materials for which not all of the principle elements of the permeability and permittivity tensors have the same sign. We find that a wide variety of effects can be realized in such media, including negative refraction, near-field focusing and high impedance surface reflection. In particular a bi-layer of these materials can transfer a field distribution from one side to the other, including near-fields, without requiring internal exponentially growing waves.

The range of available electromagnetic material properties has been broadened by recent developments in structured media, notably Photonic Band Gap materials and metamaterials. These media have allowed the realization of solutions to Maxwell’s equations not available in naturally occurring materials, fueling the discovery of new physical phenomena and the development of devices. Photonic crystals, for example, can modify the radiative density of states associated with nearby electromagnetic sources. Optical effects such as superradiance John+Quang , enhanced or inhibited spontaneous emission Koenderink , ultrarefraction Kosaka  and even negative refraction notomi ; gralak  have been predicted for various photonic lattice configurations.

Photonic crystal effects typically occur when the wavelength is on the same order or smaller than the lattice constant of the crystal. Metamaterials, on the other hand, have unit cell dimensions much smaller than the wavelength of interest. A homogenization process can thus be applied, allowing the otherwise complicated composite medium to be described conveniently by a permittivity tensor ($`\epsilon `$) and a permeability tensor ($`\mu `$) Halevi , rather than by band diagrams.

In 2000, it was shown experimentally that a metamaterial composed of periodically positioned scattering elements, all conductors, could be interpreted as having simultaneously a negative effective $`\epsilon `$ and a negative effective $`\mu `$ smithPRL . Such a medium had been previously shown by Veselago to be consistent with Maxwell’s equations veselago , but had never been demonstrated in a naturally occurring material or compound. A medium with simultaneously isotropic and negative $`\epsilon `$ and $`\mu `$ supports propagating solutions whose phase and group velocities are antiparallel; equivalently, such a material can be rigorously described as having a negative index of refractionsmithkroll . An experimental observation of negative refraction was reported using a metamaterial composed of wires and split ring resonators deposited lithographically on circuit board material shelby .

The prospect of negative refractive materials has generated considerable interest, as this simply stated material condition suggests the possibility of extraordinary wave propagation phenomena, including near-field focusing perfectLens . So remarkable have been the claims surrounding negative refraction, that some researchers have been prompted to examine critically the achievability of negative refraction in existing metamaterials Valanju ; insaneGarcia . While such concerns might appear relevant in the context of frequency-dispersive materials, the interpretation of these structured materials as negative refractive has been entirely consistent with all aspects of reported experimental data smithPRL ; shelby , as well as with numerical simulations of both monochromatic and modulated beams kongBeam ; modulate . For the purposes of this Letter, we thus assume that the descriptions presented here will be applicable to real materials, although the extent to which this is true remains a topic of active pursuit.

Lindell et al. lindell  have shown that the property of negative refraction is not confined to materials with negative definite $`ϵ`$ and $`\mu `$, but can be expected to occur in certain classes of uniaxially anisotropic media. In addition to exhibiting negative refraction, some classes of anisotropic media have very unusual dispersion relations, characterized for example by hyperbolic dispersion curves. The range of electromagnetic wave propagation behavior expected for such materials greatly extends beyond that available in isotropic negative refractive materials, and motivates our pursuit here to find configurations that take advantage of the unique properties.

To simplify the proceeding analysis, we assume a linear material with no magnetoelectric coupling, so that the media can be fully described by a permittivity, $`\epsilon `$, and permeability, $`\mu `$, tensor. \[Bianisotropy effects have typically played a minor role in the overall response of the experimental metamaterials, and can be mitigated by design Marques .) We further assume that these tensors are simultaneously diagonalizable, having the form

$$\epsilon =\left(\begin{array}{ccc}\epsilon _x& 0& 0\\ 0& \epsilon _y& 0\\ 0& 0& \epsilon _z\end{array}\right)\mu =\left(\begin{array}{ccc}\mu _x& 0& 0\\ 0& \mu _y& 0\\ 0& 0& \mu _z\end{array}\right).$$
(1)
Metamaterials can be readily constructed that closely approximate these $`\epsilon `$ and $`\mu `$ tensors, with elements of either algebraic sign. In fact, the scattering elements comprising the metamaterials used to demonstrate negative refraction shelby  are appropriate building blocks for the classes of materials to be discussed here, (see Fig. 1.)

We are interested in an anisotropic medium in which not all of the principle components of the $`\epsilon `$ and $`\mu `$ tensors have the same sign. For brevity, we refer to such a medium as *indefinite*. We will consider layered media with surfaces normal to one of the principle axes, which we define to be the $`z`$-axis. We demonstrate our analysis using a plane wave with the electric field polarized along the $`y`$-axis, though it is generally possible to construct media that are polarization independent, or exhibit different classes of behavior for different polarizations.

$$𝐄=\widehat{𝐲}e^{i\left(k_xx+k_zz\omega t\right)}.$$
(2)
Plane wave solutions to Maxwell’s equations with this polarization have $`k_y=0`$ and satisfy:

$$k_z^2=\epsilon _y\mu _x\frac{\omega ^2}{c^2}\frac{\mu _x}{\mu _z}k_x^2.$$
(3)
Since we have no $`x`$ or $`y`$ oriented boundaries or interfaces, real exponential solutions, which result in field divergence when unbounded, are not allowed in those directions. Thus $`k_x`$ is restricted to be real. Also, since $`k_x`$ represents a variation transverse to the surfaces of our layered media, it is conserved across the layers, and naturally parameterizes the solutions.

In the absence of losses, the sign of $`k_z^2`$ can be used to distinguish the nature of the plane wave solutions. $`k_z^2>0`$ corresponds to real valued $`k_z`$ and propagating solutions. $`k_z^2<0`$ corresponds to imaginary $`k_z`$ and exponentially growing or decaying (evanescent) solutions. When $`\epsilon _y\mu _z>0`$, there will be a value of $`k_x`$ for which $`k_z^2=0`$. This value, which we denote $`k_c`$, is the cutoff wave vector separating propagating from evanescent solutions. From Eq. (3) this value is $`k_c=\frac{\omega }{c}\sqrt{\epsilon _y\mu _z}.`$ We identify four classes of media based on their cutoff properties:

|  | Media Conditions | | Propagation |
| --- | --- | --- | --- |
| Cutoff | $`\epsilon _y\mu _x>0`$ | $`\mu _x/\mu _z>0`$ | $`k_x<k_c`$ |
| Anti-Cutoff | $`\epsilon _y\mu _x<0`$ | $`\mu _x/\mu _z<0`$ | $`k_x>k_c`$ |
| Never Cutoff | $`\epsilon _y\mu _x>0`$ | $`\mu _x/\mu _z<0`$ | all real $`k_x`$ |
| Always Cutoff | $`\epsilon _y\mu _x<0`$ | $`\mu _x/\mu _z>0`$ | no real $`k_x`$ |

Note the analysis presented here is carried out at constant frequency, and that the term *cutoff* always refers to the transverse component of the wave vector, $`k_x`$, not the frequency, $`\omega `$. Iso-frequency contours, $`\omega \left(𝐤\right)=const`$, show the required relationship between $`k_x`$ and $`k_z`$ for plane wave solutions (Fig. 2).

To understand wave refraction at an interface between vacuum and an indefinite medium, we must first examine the general relationship between the directions of energy and phase velocity for waves propagating within an indefinite medium. We can accomplish this by calculating the group velocity, $`𝐯_g_𝐤\omega \left(𝐤\right)`$. $`𝐯_g`$ specifies the direction of energy flow for the plane wave, and is not necessarily parallel to the wave vector. $`_𝐤\omega \left(𝐤\right)`$ must lie normal to the iso-frequency contour, $`\omega \left(𝐤\right)=`$ constant, but may lie on either side, as illustrated in Fig. 2. To determine the correct direction, we utilize the dispersion relation in Eq. (3) to calculate $`_𝐤\omega \left(𝐤\right)`$. Performing an implicit differentiation of Eq. (3) leads to a result for the gradient that does not require square root branch selection, removing any sign confusion. To obtain physically meaningful results, a causal, dispersive response function, $`\xi \left(\omega \right)`$, must be used to represent the negative components of $`\epsilon `$ and $`\mu `$, since these components are necessarily dispersive landau10 . The response function should assume the desired (negative) value at the operating frequency, and satisfy the causality requirement that, $`\left(\xi \omega \right)/\omega 1`$ landau10 ; smithkroll . Combing this with the derivative of Eq. (3) determines which of the two possible normal directions applies, without specifying a specific functional form for the response function. Fig. 2 relates the direction of the group velocity to a given material property tensor sign structure.

Having calculated the energy flow direction, we can determine the refraction behavior of indefinite media by applying two rules: (1) the transverse component of the wave vector, $`k_x`$, is conserved across the interface, and (2) energy carried into the interface from free space must be carried away from the interface inside the media; i.e., the normal component of the group velocity, $`v_{gz}`$, must have the same sign on both sides of the interface. Fig. 2 shows typical refraction diagrams for the three types of media that support propagation.

Despite the interesting refraction properties associated with indefinite media, a finite slab of indefinite material will in general present a significant impedance mismatch to waves incident from vacuum. We can, however, increase the transparency at a given frequency by combining a layer of indefinite material with a second layer possessing the same dispersion and impedance but opposite refractive index. Such a compensated bilayer will have inter-layer relative impedance near unity, and will exhibit zero phase shift from one side to the other, for all plane wave components.

To illustrate the possibilities associated with compensated bilayers of indefinite media, we first recall that a motivating factor in the recent metamaterials effort has been the prospect of near-field focusing. Pendry showed theoretically that a planar slab with isotropic $`\epsilon =\mu =1`$, could act as a lens with resolution well beyond the diffraction limit. It is difficult, however, to realize significant sub-wavelength resolution with an isotropic negative index material, as the required exponential growth of the large $`k_x`$ field components across the negative index lens leads to extremely large field ratiosziolkowski . Sensitivity to material loss and other factors can significantly limit the sub-wavelength resolution.

By combining positive and negative refracting layers of *never cutoff* indefinite media, we can produce a compensated bilayer that accomplishes near-field focusing in a similar manner to the perfect lens, but with significant advantages. Fig. 2 indicates that for the same incident plane wave, the $`z`$-component of the transmitted wave vector is of opposite sign for these two materials. Combining appropriate lengths of these materials results in a composite indefinite medium with unit transfer function. We can see this quantitatively by considering the general expression for the transfer function of a bilayer.

$$T=8\left[\begin{array}{c}e^{i\left(\varphi +\psi \right)}\left(1Z_0\right)\left(1+Z_1\right)\left(1Z_2\right)+\\ e^{i\left(\varphi \psi \right)}\left(1Z_0\right)\left(1Z_1\right)\left(1+Z_2\right)+\\ e^{i\left(\varphi +\psi \right)}\left(1+Z_0\right)\left(1Z_1\right)\left(1Z_2\right)+\\ e^{i\left(\varphi \psi \right)}\left(1+Z_0\right)\left(1+Z_1\right)\left(1+Z_2\right)\end{array}\right]^1$$
(4)
The relative effective impedances are defined as

$$Z_0=\frac{q_{z1}}{\mu _{x1}k_z},Z_1=\frac{\mu _{x1}}{\mu _{x2}}\frac{q_{z2}}{q_{z1}},Z_2=\mu _{x2}\frac{k_z}{q_{z2}}.$$
(5)
where $`𝐤`$, $`𝐪__1`$ and $`𝐪_2`$ are the wave vectors in vacuum and the first and second layers of the bilayer respectively. The individual layer phase advance angles are defined as $`\varphi q_{z1}L_1`$ and $`\psi q_{z2}L_2`$, where $`L_1`$ is the thickness of the first layer and $`L_2`$ is the thickness of the second layer. If the signs of $`q_{z1}`$ and $`q_{z2}`$ are opposite as mentioned above, the phase advances across the two layers can be made equal and opposite, $`\varphi +\psi =0`$. If we further have the condition that the two layers are impedance matched to each other, $`Z_1=1`$, then Eq. (4), reduces to $`T=1`$. In the absence of loss, the material properties can be chosen so that this occurs for all values of the transverse wave vector, $`k_x`$. A transfer function with some realistic loss added is shown in Fig 3.

A proposed implementation is shown in Fig. 1. The elements shown in the top and bottom of the figure will implement media that focuses electric $`y`$-polarized and magnetic $`y`$-polarized waves respectively. Combining the two structures results in a bilayer that focuses both polarizations and is $`x`$-$`y`$ isotropic. The materials are formed from split ring resonators and wires with numerically and experimentally confirmed effective material properties shelby . Each split ring resonator orientation implements negative permeability along a single axis, as does each wire orientation for negative permittivity.

While compensated bilayers of indefinite media exhibit reduced impedance mismatch to free space and high transmission, uncompensated semi-infinite indefinite media exhibit unique high reflection properties. Aside from *cutoff* materials, all other classes of indefinite media have a reflection coefficient amplitude near unity for incident propagating waves. The phase of the reflection coefficient, however, varies, as illustrated here for positive refracting *anti-cutoff* media. The reflection coefficient for electric $`y`$-polarization is given by

$$\rho =\frac{\mu _xk_zq_z}{\mu _xk_z+q_z}$$
(6)
where $`𝐤`$ and $`𝐪`$ are the wave vectors in vacuum and the media respectively. For unit magnitude *anti-cutoff* material we have, from Eq. (3),

$$q_z^2=\frac{\omega ^2}{c^2}+k_x^2=k_z^2$$
(7)
Thus $`q_z=\pm ik_z`$. The correct sign for positive refracting media, $`+`$, is determined by the requirement that the fields must not diverge in the domain of the solution. We then have

$$\rho =\frac{1i}{1+i}=i$$
(8)
The magnitude of the reflection coefficient is unity, with a phase of -90 for propagating modes of all incident angles. An electric dipole antenna placed $`\lambda /8`$ away from the surface would thus be enhanced by interaction with this “mirror” surface. Customized reflecting surfaces are of practical interest, as they enhance the efficiency of nearby antennas, while at the same time providing shielding yablonovitch ; ma . In this case, the invariance of the reflection phase over all angles of incidence may be advantageous. Furthermore, surface modes, which correspond to poles in Eq. 6, are not supported on this interface. Surface modes require evanescent solutions on both sides of the interface. It is possible, as shown in Fig. 2, to have no overlap of evanescent solutions (dashed line) between *cutoff* and *anti-cutoff* media. In many communications applications the energy lost to the excitation of surface modes is undesirable, as it represents loss of signal.

In conclusion, we have begun to explore the properties of media with indefinite $`\epsilon `$ and $`\mu `$ tensors. Indefinite media are governed by hyperbolic dispersion relations, previously found in more exotic situations, such as relativistic moving media kongBook . Consideration of layered structures has led to useful and interesting reflection and refraction behavior, including a new mechanism for sub-diffraction focusing. We note that neither the analysis nor the fabrication of these media is complicated, and thus anticipate other researchers will quickly assimilate the principles and design structures with unique and technologically relevant properties.

We thank Claudio Parazzoli (Phantom Works, Boeing) for motivating this work. This work was supported by DARPA through grants from ONR (Contract No. N00014-00-1-0632) and AFOSR (Contract No. 78535A/416250/440000) and a grant from AFOSR (Contract Number F49620-01-1-0440).


text:
# 1 SS433/W50 Review

## 1 SS433/W50 Review

SS 433 is the famous galactic microquasar known by its two-sided precessing mildly relativistic (0.26c) jets (Margon 1984). The question whether the compact object is a neutron star or a black hole remains unanswered, in spite of its discovery over 2 decades ago. ASCA observations revealed Doppler blue- and red-shifted emission lines (Kotani et al. 1996) indicating a jets’ scale of $``$ 1013 cm. Recent Chandra observations with the gratings (Marshall et al. 2002) gave a mass outflow rate of 1.5$`\times `$10-7 M yr-1 and a kinetic power of 3.2$`\times `$1038 erg s-1 (at an assumed distance of 4.85 kpc). The latter values are small compared to previous estimates of $``$ 10-6 M yr-1 and 1039-1041 erg s-1 (possibly because SS 433 was in a low state during the Chandra observations).

SS 433 is near the center of the supernova remnant (SNR) W50. Observations in the radio (Dubner et al. 1998), millimeter (Durouchoux et al. 2000; Sood et al. 2002 in preparation), infrared (Fuchs, Koch-Miramond, & Abraham 2001), optical (Mazeh et al. 1983), and X-rays (Safi-Harb & Petre 1999 and refs therein) show that the SS 433 jets are interacting with the surounding inhomogeneous medium, causing the unusual morphology and the X-ray lobes of W50.

SS 433 was observed by the RXTE at several occasions. Here we present the preliminary results of our spectral and timing study. A more detailed analysis and interpretation of our results will be presented elsewhere. We also refer the reader to a companion paper (Kotani et al. this volume) highlighting the multiwavelength 2001 campaign, and to several other SS 433 papers in this volume highlighting the most recent multi-wavelength observations: Namiki et al. for recent Chandra observations, Fabrika et al. for the optical observations, Fuchs et al. for the infrared observations suggesting a Wolf-Rayet star origin, Blundell et al. and Paragi et al. for the evidence of an equatorial outflow, Migliari et al. for thermal reheating in the jet, Trushkin et al. for the 6-day modulation in the quite radio emission, and Chakrabarti et al. for the theoretical implications.

## 2 Radio and X-ray anti-correlation

RXTE observed SS 433 in 1996 ($`\varphi `$=0.83–1.12; $`\psi `$=0.78–0.81), 1998 ($`\varphi `$=0.01–1.0; $`\psi `$=0.1–0.2), and 2001 ($`\varphi `$=0.1–1.3; $`\psi `$=0.29–0.38); where $`\varphi `$ and $`\psi `$ are the optical and precession phases, respectively. A radio flare occurred mid March 1998 (around MJD 50890) and in November 2001 (MJD 52215 and 52235), as shown in Fig. 1 by the arrows. Luckily the 1998 flare (detected with the NRAO GBI) was covered by two short RXTE observations, and the November 2001 flare detected with the RATAN 600-m telescope (Kotani & Trushkin 2001) was monitored by RXTE for 2 weeks (except for November 18, 2001).

The PCA spectra are well fitted with a thermal bremsstrahlung model (or a power law model at occasions, see below), modified by interstellar absorption, and with two Fe-lines: A broad line and a narrow line. In Fig. 1, we show the hard X-ray flux variations and its anti-correlation with the radio flares that occurred in 1998 (left) and 2001 (right). In Fig. 2, we show the variations of the spectral parameters, $`kT`$ and $`\mathrm{\Gamma }`$, fixing $`N_H`$ at 0.7$`\times `$1022 cm-2. We find that 1) the X-ray flux decreases during the flare, 2) the spectrum softens at the onset of the flare then hardens shortly after, and 3) the spectrum switches from thermal to power law right after the flare in 1998. The narrow line variations (believed to originate from the jets) are also consistent with the redshifts or blueshifts determined from optical spectroscopy (Kotani et al. 2002).

## 3 Fast time variability

Previous observations revealed no variability on time-scales shorter than 300 s. A search for variability in the 0.5-300 s range in the ROSAT data (Safi-Harb 1997, see Fig. 3, left) using the method of autocorrelation function shows flickering around 3–10 s. However, this variability does not appear consistently in all the observation segments, and the time scale varies from an observation segment to the other. Using the PCA Nov. 2001 data, we also found fast time variability (Fig.3, right), which we confirmed using previous archival PCA observations. The time scale of 50–100 s indicates a length scale $``$1012cm. Fast time variability is expected when the source enters a highly non-stationary regime expected from super-critical accretion into the surface of a neutron star.

## 4 Progress and the future

We have presented the first evidence for the hard X-ray/radio anti-correlation, and for fast time variability. The spectral behavior mimics the behavior seen in other microquasars (Mirabel & Rodriguez 1999). The short time scales indicate that we are, for the first time, probing the high-energy emission regions closer to the compact object. Chakrabarti et al. (2002) show that non-steady shocks in sub-Keplerian accretion flow provides the basic timescale of the ejection interval.

SS 433 remains unique compared to other microquasars: 1) the ratio of its X-ray flux to the jets power is much smaller (10-4 compared to 0.1 for GRS 1915+105), 2) it’s the only microquasar with large scale X-ray and radio lobes resulting from the interaction between the jets and the surrounding medium/SNR (see however the discovery by Corbel et al., this volume, of large scale jets from a black hole), and 3) it’s the only source confirmed to have precessing jets with Doppler shifted emission lines revealing a complex system with baryonic jet matter. The current picture that best explains the properties of this enigmatic object is a binary system embedded in an expanding thick disk that is fed by the wind from a super-Eddington accretion rate (see Gies et al. 2002). The companion star is most likely a Wolf-Rayet star dumping some 10-4 M yr-1 into the compact object. This would explain the large infrared flux, the low X-ray luminosity, and the difficulty to detect any pulsations.

Several burning questions still need to be addressed with further observations: 1) the discrepancy in the distance derived for SS 433 and W50 needs also to be resolved to derive better estimates for the mass outflow and the jets power; 2) the nature of the compact object remains ambiguous. While the majority of published papers point to an underlying black hole (e.g. Zwitter & Calvani 1989), a later estimate by D’Odorico et al. (1991) indicates a neutron star origin. Interestingly, even in the latter scenario, the neutron star’s mass is unusually low. If the estimated mass is indeed 0.8$`\pm `$0.1 $`M_{}`$, then SS 433 would join the growing class of low-mass ‘neutron’ stars (Gondek-Rosinska, Kluzniak, & Stergioulas 2002).

SS 433 remains unique. While other objects have been claimed to be SS 433-like (see reviews in this volume by Safi-Harb, Fender, and the study of CI-Cam and V4641 Sgr by Rupen), none of them is an SS 433–twin. SNRs offer a promising laboratory to search for binary compact objects (potential SS 433’s). Simultaneous monitoring of SS 433 in the radio, infrared, high energy X-rays and gamma-rays is needed to study in more detail the spectral and timing behavior discussed here, in comparison with other microquasars. Most importantly, a reliable estimate of the mass function will hopefully solve a 25-year old puzzle.

## Acknowledgments

SSH thanks NSERC for the support and the organizers for the invitation.


text:
# LHC(CMS) SUSY discovery potential for nonuniversal gaugino and squark masses and the determination of the effective SUSY scale

## Abstract

We review the results on the investigation of the LHC(CMS) SUSY discovery potential for the case of nonuniversal gaugino and squark masses.

1 Institute for High Energy Physics, Protvino, Russia

One of the LHC supergoals is the discovery of the supersymmetry. In particular, it is very important to investigate a possibility to discover strongly interacting superparticles (squarks and gluino). In Ref. (see also Refs.) the LHC squark and gluino discovery potential has been investigated within the minimal SUGRA-MSSM framework  where all sparticle masses are determined mainly by two parameters: $`m_0`$ (common squark and slepton mass at GUT scale) and $`m_{\frac{1}{2}}`$ (common gaugino mass at GUT scale). The signature used for the search for squarks and gluino at LHC is $`(n0)isolatedleptons+(m2)jets+E_T^{miss}`$ events. The conclusion of Ref.  is that LHC is able to detect squarks and gluino with masses up to $`(22.5)TeV`$.

Despite the simplicity of the SUGRA-MSSM framework it is a very particular model. The mass formulae for sparticles in SUGRA-MSSM model are derived under the assumption that at GUT scale ($`M_{GUT}210^{16}`$ GeV) soft supersymmetry breaking terms are universal. However, in general, we can expect that real sparticle masses can differ in a drastic way from sparticle masses pattern of SUGRA-MSSM model due to many reasons, see for instance refs. . Therefore, it is more appropriate to investigate the LHC SUSY discovery potential in a model-independent way for general case of nonuniversal gaugino and squark masses.

In this report we review our previous results - on the investigation of the LHC(CMS) SUSY discovery potential for the case of nonuniversal squark and gaugino masses. We also give some new results on the determination of the effective SUSY scale for models with nonuniversal gaugino masses.

The decays of squarks and gluinos depend on the relation among squark and gluino masses. For $`m_{\stackrel{~}{q}}>m_{\stackrel{~}{g}}`$ squarks decay mainly into gluino and quarks and gluino decays mainly into a quark-antiquark pair and a gaugino. For $`m_{\stackrel{~}{q}}<m_{\stackrel{~}{g}}`$ gluino decays mainly into squarks and quarks, whereas squarks decay mainly into quarks and a gaugino. The next-to-lightest gauginos have several leptonic decay modes giving a lepton and missing energy, for instance:

$$\stackrel{~}{\chi }_1^\pm \stackrel{~}{\chi }_1^0+l^\pm +\nu ,$$
(1)
$$\stackrel{~}{\chi }_2^0\stackrel{~}{\chi }_1^0+l^\pm l^{}$$
(2)
As a result of chargino and second neutralino leptonic decays, besides the classical signature $`(n2,3,4)jetsplusE_T^{miss}`$ signatures such as $`(n2,3,4)jetsplus(m1leptons)plusE_T^{miss}`$ with leptons and jets in final state arise.

In our concrete calculations all SUSY processes, with full particle spectrum, couplings, production cross section and decays have been generated with ISASUSY code . The Standard Model backgrounds have been generated by PYTHIA code . Our simulations have been made with parametrized detector responses based on detailed detector simulations. Namely, the CMS(Compact Muon Solenoid) detector fast simulation code CMSJET  has been used. It appears that the following SM processes give the main contribution to the background:

$`W+jets,Z+jets,t\overline{t},WZ,ZZ,b\overline{b}andQCD(22)reactions`$ .

As it has been mentioned previously in our study we have considered as signatures $`(nl)`$ jets plus $`(mk)`$ isolated leptons plus $`E_T^{miss}`$, where $`l=2,3,4`$ and $`k=0,1,2,3,4`$. More detailed information on the used cuts can be found in Refs.-.

For squark and gluino pair production we have considered 4 different kinematic regions:

* A. $`m_{\stackrel{~}{g}}m_{\stackrel{~}{q}}`$
* B. $`m_{\stackrel{~}{q}}m_{\stackrel{~}{g}}`$
* C. $`m_{\stackrel{~}{q}}m_{\stackrel{~}{g}}`$, $`m_{\stackrel{~}{q}}>m_{\stackrel{~}{g}}`$
* D. First 2 squark generations are heavy $`m_{\stackrel{~}{q}_{1,2}}1TeV`$. Only 3 rd squark generation is relatively light.

We have found that the LHC(CMS) SUSY discovery potential depends rather strongly on the relation among $`m_{\stackrel{~}{\chi }_1^0}`$ and $`min(m_{\stackrel{~}{g}},m_{\stackrel{~}{q}})`$ and it decreases with the increase of the LSP mass $`m_{\stackrel{~}{\chi }_1^0}`$. For LSP mass $`m_{\stackrel{~}{\chi }_1^0}`$ close to $`min(m_{\stackrel{~}{g}},m_{\stackrel{~}{q}})`$ it is possible to detect SUSY for $`min(m_{\stackrel{~}{g}},m_{\stackrel{~}{q}})(1.21.5)TeV`$ (cases A,B,C). For the case D it is possible to detect SUSY with $`m_{\stackrel{~}{q}_{1,2}}800GeV`$. Our results are presented in Figs.(1-5).

Gluinos and squarks are strongly produced at the LHC, so it is important to find a variable that measures the produced mass for events with missing transverse energy. In Ref. 111See also recent Ref. on this subject. a variable that works rather well within SUGRA-MSSM model has been defined as the sum of the missing energy and the $`p_T`$’s of the first four jets,

$$M_{eff}=E_T^{miss}+p_{T,1}+p_{T,2}+p_{T,3}+p_{T,4}$$
(3)
As it has been mentioned in Ref. in SUGRA-MSSM model the distribution of SUSY cross section $`\frac{d\sigma _{SUSY}}{dM_{eff}}`$ has maximum at some point $`M_{SUSY}M_{peak}`$ and the value of $`M_{SUSY}`$ coincides within 10 percent accuracy with $`min(m_{\stackrel{~}{g}},m_{\stackrel{~}{q}})`$, namely

$$M_{SUSY}M_{peak}min(m_{\stackrel{~}{g}},m_{\stackrel{~}{q}})$$
(4)
We have investigated the dependence of $`M_{SUSY}`$ on squark, gluino and LSP masses in general case of arbitrary relations among LSP, squark and gluino masses. Our results are presented in Figs.(6-8). As it follows from our results in general case for LSP mass $`m_{\stackrel{~}{\chi }_1^0}m_{\stackrel{~}{g}}`$ there is no linear dependence of $`M_{SUSY}`$ on gluino or squark masses. The value of $`M_{SUSY}`$ is nonlinear and nontrivial function on gluino, squark and LSP masses. So, in general case it would be rather nontrivial to determine more or less accurately the squark and gluino masses with the help of $`M_{SUSY}`$.

In conclusion we would like to stress that even for the most difficult case when LSP mass is rather heavy and close to squark or gluino masses the LHC is able to discover SUSY for the most interesting from the theoretical point of view case when sparticle masses are lighter than $`1TeV`$.

This work has been supported by CERN-INTAS Grant 99-0377 and RFFI Grant N02-01-00601.


text:
# Temperature inversion in granular fluids under gravity

## Abstract

We study, via hydrodynamic equations, the granular temperature profile of a granular fluid under gravity and subjected to energy injection from a base. It is found that there exists a turn-up in the granular temperature and that, far from the base, it increases linearly with height. We show that this phenomenon, observed previously in experiments and computer simulations, is a direct consequence of the heat flux law, different form Fourier’s, in granular fluids. The positive granular temperature gradient is proportional to gravity and a transport coefficient $`\mu _0`$, relating the heat flux to the density gradients, that is characteristic of granular systems. Our results provide a method to compute the value $`\mu _0`$ for different restitution coefficients. The theoretical predictions are verified by means of molecular dynamics simulations, and the value of $`\mu _0`$ is computed for the two dimensional inelastic hard sphere model. We provide, also, a boundary condition for the temperature field that is consistent with the modified Fourier’s law.

In a granular system a small amount of energy is lost at each collision between particles. This energy loss at the microscopic level leads to macroscopic dynamics quite different from those of systems with particles undergoing elastic collisions. One of the surprising consequences of the dissipative dynamics is the granular temperature inversion in fluidized granular systems under gravity. The granular temperature —the temperature from now on— is defined proportional to the mean kinetic energy per particle in the reference frame of the fluid. It has been observed that an open system subjected to a permanent energy injection from the base, exhibits a minimum of temperature at some distance from the boundary, and from there on temperature increases with height. In spite of the minimum the energy flux continues always pointing upwards. This fact confirms the existence of a heat flux law, different from Fourier’s law, in these granular fluids. The heat flux is defined, for granular fluids, as the kinetic energy flux in the reference frame of the fluid. Nevertheless, the fact that the kinetic energy per particle, that is the granular temperature, can increase with increasing height from the base is unexpected, given the inelastic nature of collisions.

This temperature turn-up can be indirectly observed in experiments. For example, in  the authors report a temperature increase 30 times greater than expected due to the slow convective movement they observe in their experiments. Their conclusion is that this effect may not be attributed to convection currents, and suggest that it could be due to density gradients. The temperature minimum has also been observed in computer simulations  and a recent article accounts for this minimum with an hydrodynamic approach .

In this article, we deduce the presence of the minimum. Instead of solving the whole hydrodynamic equations for dilute granular gases as in , we show, using simple hydrodynamic arguments, that the temperature minimum has its origin in the modification of the Fourier law for the heat flux in granular fluids . Also, unlike in , these simple hydrodynamic arguments predict, far from the base, a linear increase of temperature with height with a slope that can be deduced directly from the theory. This linear dependence is verified using molecular dynamic simulations of a simple granular model.

Consistently with these results, we establish the appropriate boundary condition for the temperature of an unbounded system, and for a system bounded by an upper adiabatic wall. Finally we provide a method to measure experimentally the transport coefficient associated with the non-Fourier heat flux law. All these results are compared with molecular dynamics simulations for a bi-dimensional system.

Hydrodynamic equations (that is, closed equations for a reduced number of macroscopic fields), resembling those of elastic systems, have been derived for the description of granular fluids . Even though the validity of such description might be controversial , hydrodynamic equations together with appropriate constitutive relations have shown to be valid in quasi-elastic low-density stationary states.

The difference between the elastic and the granular hydrodynamic equations is an extra sink term in the energy balance equation which accounts for the energy dissipation. Also, constitutive relations are quite different in granular materials and in elastic systems. Fourier’s Law, as mentioned, does not hold for these dissipative systems, but rather the heat flux follows a law in the form

$$𝐉=kT\mu n,$$
(1)
$`𝐉`$ being the heat flux, $`k`$ and $`\mu `$ the transport coefficients, $`T`$ the granular temperature, and $`n`$ the number density. For small dissipation and low density, kinetic theory predicts that $`k`$ and $`\mu `$ must be of the form

$`k`$ $`=`$ $`k_0\sqrt{T},`$
$`\mu `$ $`=`$ $`\mu _0{\displaystyle \frac{T^{3/2}}{n}},`$ (2)
where $`k_0`$ and $`\mu _0`$ depend only on the dissipative coefficient and $`\mu _0`$ vanishes in the elastic limit. In quasi-elastic systems, $`\mu _0k_0`$ .

Experimentally, fluidized stationary states can be achieved by means of a vibrating base with some amplitude and frequency. In the low amplitude and high frequency limit, quasi-elastic systems reach non-equilibrium stationary states in which the vibrating base plays the role of a stationary boundary. A thermal wall is a useful theoretical model to describe this kind of boundary: each time a particle collides with the wall, it comes out with a velocity sorted out from a Maxwellian distribution.

While the results presented here should not depend on the particular mechanism of the energy injection, for simplicity, we will focus on an open system fluidized by a thermal base at the bottom. The system is kept bounded by a gravitational field pointing downwards with intensity $`g`$.

Such a system may exhibit different regimes. If density, gravity and dissipation are not too large, the system reaches a stationary state: the conductive regime. This state is characterized by a vanishing velocity field and density, temperature, and heat flux fields depending only on the vertical coordinate $`z`$. For increasing density, dissipation and/or gravity the system starts developing different instabilities such as Rayleigh Benard-like convection and solidification .

In the conductive regime, the energy dissipation at collisions induces a vertical temperature gradient. Energy is injected through the base and dissipated in the volume, with the corresponding decrease of the temperature with height. A snapshot of a typical configuration in the conductive regime is shown in Fig. 1.

Close to the base, we can state that the temperature decreases with height due to the high dissipation rate. Depending on the gravity acceleration and the inelasticity , the negative temperature gradient, is sometimes accompanied by a density increase which, at some extreme conditions, can lead to clustering and solidification.

Far from the base, on the other hand, we know that the density decreases by the action of gravity. Indeed, if the height is large enough ($`z45`$ in Fig. 1) we find just a few particles moving mainly in free flight, so that we can safely consider that there is a height above which the pressure is given by the ideal gas equation of state, $`p=nT`$.

The momentum balance equation states that the pressure obeys the barometric law $`\frac{dp}{dz}=ng`$, so that together with the ideal gas expression for $`p`$ we obtain,

$$n\frac{dT}{dz}+T\frac{dn}{dz}=ng.$$
(3)
Inserting the expression for the density gradient derived from (3) into (1) and using the dilute limits for the transport coefficients (2), the heat flux far from the base is given by

$$J=\left(k_0\mu _0\right)\sqrt{T}\frac{dT}{dz}+\mu _0g\sqrt{T}.$$
(4)
As the density goes to zero with increasing height, the heat flux must consequently vanish as there are basically no particles to transport energy. It can be assumed then that there is a height where $`J\mu _0g\sqrt{T}`$. Above this height the heat flux in Eq. (4) can be neglected leading to the following relation for the temperature gradient as $`z`$ goes to infinity

$$\frac{dT}{dz}=\frac{\mu _0g}{k_0\mu _0}.$$
(5)
That is, the temperature gradient is a constant at high positions so that, far from the heating base, we can predict that the temperature field increases linearly with height with a slope proportional to $`g`$ and $`\mu _0`$.

Furthermore, as the temperature near the base decreases with height, it can be deduced, by continuity, that the temperature field has a minimum. Above the minimum, the temperature increases until it reaches the linear asymptotic dependence given by Eq. (5).

This prediction is fully corroborated in our MD simulations, as observed in Fig. 2 where the measured temperature profile has been plotted for two different systems. In both of them we clearly observe the minimum and the further linear increase of the temperature.

In our MD simulations, the inelastic hard sphere model (IHS) in two dimensions is used as the microscopic model for the particles. Grains are modeled as smooth hard disks that dissipate energy at each collision through a constant dissipative coefficient $`q`$, related to the normal restitution coefficient $`r=12q`$, so that $`q=0`$ represents the elastic case. The simulated systems are composed of $`N`$ inelastic hard disks. In the horizontal direction $`x`$, we set periodic boundary conditions. The box width is $`L`$. At the bottom we imposed a thermal wall at temperature $`T_0`$. Due to the absence of an intrinsic energy scale in the IHS model, we fix the temperature $`T_0`$ to unity. Results for other temperatures are obtained by simple dimensional analysis. Also, units are chosen such that the disks diameters $`\sigma `$ and masses $`m`$ are set to one. In these units the dimensionless gravity is $`\widehat{g}=mg\sigma /T_0`$. The system is completely defined giving the values of $`N`$, $`L`$, $`\widehat{g}`$, and the dissipative coefficient $`q`$. For each value of the parameters we thermalize the system after which we measure the macroscopic fields, averaging over a long simulation time.

In Fig. 2 we present the results for: system A (dashed line) with parameters $`N=1120`$, $`L=89.44`$, $`q=0.01`$, and $`\widehat{g}=0.02`$; and system B (solid curve) with $`N=560`$, $`L=44.72`$, $`q=0.01`$, and $`\widehat{g}=0.03`$. The snapshot presented in Fig. 1 corresponds to the system A. Comparing both figures it can be observed that the temperature minimum (at $`z45\sigma `$) lies close to the free surface of the system but still into the bulk. The system size, under which 99% of the particles lie, is $`60\sigma `$ so that the observed linear dependence of the temperature remains up to twice the system size, proving to be the real asymptotic behavior. Similar behavior is observed in system B, that has an effective height of $`z=45\sigma `$ and the temperature minimum is located at $`z=35\sigma `$. For higher positions, the density is very low implying large statistical errors.

To test any possible finite size dependence on the temperature slope, we performed simulations for different number of particles $`N`$, keeping fixed the ratio $`N/L`$. We found that even if the profiles $`n(z)`$ and $`T(z)`$ may suffer variations, the asymptotic temperature slope remains independent of $`N`$ for $`N280`$.

We could deduce from these results that the temperature could achieve arbitrarily large values as height increases. However, our hydrodynamic description is limited because of the zero density, or infinite Knudsen number, limit at these heights. Our prediction at high positions is, nevertheless, physically meaningful as the combined effect of the temperature increase (linear) and the density decrease (almost exponential) is a net decrease of the energy density with height. So, even if the average kinetic energy per particle increases, the total energy of the system remains always finite. In summary, the prediction of the linear regime is valid up to some height after which the starting equations (hydrodynamics) are not valid, but anyway give physically meaningful predictions.

A system bounded by an upper adiabatic wall (that is, grains are reflected elasticaly) is a similar, although not identical, problem to that of the unbounded system. At an adiabatic boundary, the heat flux is exactly zero so that Eq. (5) is the exact expression for the temperature gradient at the wall if the density at the wall is low enough. The temperature gradient at the elastic wall is then positive, provided $`\mu _0k_0`$. This ensures the existence of a minimum, as the temperature gradient is negative at the base. The minimum of temperature induced by an upper elastic wall has been observed in one dimensional simulations .

It is clear that these results concern closely the boundary condition for the temperature field that must be used to solve the hydrodynamic equations in granular systems. In an unbounded system under gravity, Eq. (5) must be imposed as the boundary condition for the temperature at infinity. The same boundary condition applies to any adiabatic wall, contrary to the usual condition $`dT/dz=0`$ that has traditionally been used in computations of granular systems.

Those experiments where the stationary temperature inversion has been observed , have the normal lack of statistics at high positions and the linear increase of temperature can not be easily detected, specially if this is not the final purpose of the experiment. We want to stress that the measurement of the temperature slope, either experimentally or in simulations, provides a method to obtain the transport coefficient $`\mu _0`$ for small dissipation coefficients.

As an example, we have calculated the value of $`\mu _0`$, at the lowest order in the dissipative coefficient $`q`$, for the two dimensional inelastic hard sphere model. We consider systems with $`N=560`$ and $`L=44.72`$. In order to avoid instabilities –especially the convection described in Ref. – systems with low dissipative coefficients and low gravity values have to be considered. Our simulations have been done with $`\widehat{g}=0.01`$ and $`\widehat{g}=0.02`$, and $`0.002q0.02`$.

For these small values of $`q`$, we can use $`(k_0\mu _0)k_0`$ to first order in $`q`$ . Also, in the same order of approximation, $`k_0`$ can be set to that of the elastic hard disk system, $`k_0=2/\sqrt{\pi }`$, neglecting the linear correction in $`q`$. Then, $`\mu _0`$ is given by

$$\mu _0=\frac{2}{\sqrt{\pi }g}\frac{dT}{dy},$$
(6)
expression that is valid only up to order $`q`$.

In Fig. 3, we present results for $`\mu _0`$ obtained in the simulations. It can be observed that the two sets, with different values of $`\widehat{g}`$, collapse into the same curve, confirming the linearity with $`\widehat{g}`$ of the temperature gradient.

As expected, $`\mu _0`$ goes to zero in the elastic limit. The computed values of the two series are fitted to a single quadratic function in $`q`$ where it is imposed that $`\mu _0`$ must vanish for $`q=0`$. The obtained fit is

$$\mu _0=0.64q+106q^2$$
(7)
noting that only the linear term in the last fit is consistent with the approximation done. In Ref.  the authors computed $`\mu _0`$ using a different method from that described in here, with the result $`\mu _0=0.7q`$ which agrees with the result found in the present paper. Nevertheless, we want to point out that the present method for computing $`\mu _0`$ uses the temperature profile over a wide range (the region where $`T`$ increases linearly) instead of using the information over just one point as in Ref. . Therefore this method of computing $`\mu _0`$ is more accurate in the dilute case.

In the highly fluidized regime, the particular energy injection mechanism becomes irrelevant far from the energetic boundary. This implies that even in a system on a vibrating base, the method presented here should be fully valid. Also, even though our predictions are valid in the dilute case, we expect that in denser regimes there should be also a temperature turn-up and a temperature increase with height (though maybe not linear). In this case, there are explicit density dependences in the equation of state and transport coefficients, giving rise to a more complex differential equation for $`T`$ compared to (5). Nevertheless, as long as there is a non-vanishing coefficient $`\mu `$, the equation is not trivial and a temperature gradient is predicted.

In summary, we have shown that the granular temperature turn-up observed in vibrofluidized granular systems has its origin in the modification of the Fourier law for the granular heat flux. The vanishing energy flux far from the base together with the density gradient produced by the gravity imposes a positive and constant temperature gradient. Careful measurements of this gradient provide a method to compute the transport coefficient $`\mu _0`$, either experimentally or in simulations.

R.S. acknowledges the grant from MIDEPLAN, Chile, the Programa de inserción de científicos chilenos of Fundación Andes, and the support of the FONDECYT project 1010416. R.R. acknowledges Marie Curie Individual Fellowships from the European Community.