text:
# Formation and evolution of bars in disc galaxies

## 1. Introduction

Bars can form spontaneously in galactic discs (e.g. Miller, Prendergast & Quirk 1970, Hohl 1971). They then evolve over a large number of rotations, changing their length, strength, shape and angular frequency. During this period they can drive spiral or ring formation, push gas to the center-most parts of the galaxy and thus trigger starbursts and activity in the nucleus, and interact with the outer disc, bulge and/or halo by exchanging energy and angular momentum with them. I will discuss a few of these processes, relying – due to the strongly nonlinear nature of bars – mainly on results of $`N`$-body simulations.

## 2. Simulations

The results presented in sections 3 to 6 come from a preliminary analysis of well over 100 simulations run on our GRAPE-5 systems (Kawai et al. 2000). Six of these have been discussed by Athanassoula & Misiriotis (2002, hereafter AM) and Athanassoula (2002, hereafter A02), where more information on the numerical techniques can be found. Initially the disc is always exponential radially, and follows an $`sech^2`$ law with scalelength $`z_0`$ in the vertical direction. Its initial $`Q`$ value (hereafter $`Q_{init}`$) does not vary with radius. The disc mass ($`M_d`$) and initial scalelength ($`h`$) will be used as units of mass and length in the simulations. The halo is initially isotropic and non-rotating, and follows the density law

$$\rho _h(r)=\frac{M_h}{2\pi ^{3/2}}\frac{\alpha }{r_c}\frac{exp(r^2/r_c^2)}{r^2+\gamma ^2},$$
(1)
where $`\alpha `$ is a normalisation constant, $`M_h`$ is the mass of the halo and $`\gamma `$ and $`r_c`$ are scalelengths. This functional forms contains three free parameters and thus allows a large flexibility. For relatively large values of $`\gamma `$ the halo profile resembles those used by observers in rotation curve decompositions, and is thus observationally motivated. On the contrary, for sufficiently small values of $`\gamma `$ the halo rotation curve rises very steeply, thus mimicking a cusp in the halo distribution. Indeed, the core can be chosen sufficiently small to be hidden by the softening length. In all the simulations discussed here, except for the ones mentioned in the end of section 6, $`M_h`$ = 5 and $`r_c`$ = 10.

The total number of particles per simulation varies between 1 and 1.5 million. I adopted a softening of 0.0625 and a time-step of 0.015625. I assessed the numerical robustness of my results by trying in a few cases double the number of particles, different values of the softening and time step, as well as direct summation (with half the number of particles) and a non-GRAPE tree code. Various calibrations are possible. AM proposed to set the unit of length equal to 3.5 kpc and the unit of mass to 5 $`\times `$ $`10^{10}`$ $`M_{}`$. In this case $`t`$ = 500 corresponds to 7 Gyrs. Alternatively, one could use a considerably smaller value in kpc for the disc scalelength, arguing that this is the initial disc scalelength and that this length will increase with time. Thus opting for a 2 kpc scalelength would set $`t`$ = 500 to roughly 3 Gyrs. Time, in particular, may be difficult to calibrate, since numerical parameters, like the softening and the number of particles, can influence the rate at which the bar grows in the initial stages of the simulation. Thus the clock may tick differently, at least for certain parts of the evolution, in simulations than in real galaxies. For these reasons, here we will stay with computer units, allowing the reader the freedom to convert them to astronomical units according to his/her preferences and needs.

The two parameters in the initial conditions influencing most the evolution are the central concentration of the halo and $`Q_{init}`$. The former is given by the parameter $`\gamma `$, which can be thought of as the core radius of the halo.

## 3. The effect of halo central concentration

Figure 1111In the third row of panels, giving the face-on view, I used density values for the isocontours such as to show the structure in the bar region, but leaving out the disc region at larger radii. compares three simulations. The one illustrated in the left column of panels has $`\gamma `$ = 5 and its inner parts are disc dominated from the onset (upper left panel). I will hereafter refer to it as model MD, for massive disc. The middle and right panels refer to simulations with $`\gamma `$ = 0.5. Here the halo is much more centrally concentrated and contributes to the total circular velocity curve somewhat more than the disc up to roughly two disc scalelengths, and considerably more at larger radii (middle and left upper panels). I will refer to these two simulations as MH, for massive halo, and RH, for rigid halo, respectively. Indeed in simulation MH, as in MD, the halo is live, i.e. composed of particles, while in simulation RH it is rigid, i.e. is described by a force and potential imposed on the disc particles. All three simulations start off with $`Q_{init}`$ = 1.2.

The differences between the results of simulations MD and MH are very important. The bar in model MH is longer and stronger than that in model MD. Its isodensities are more rectangular-like and, seen edge-on with the bar seen side-on, it displays an ‘X’ shape, while for the same orientation MD shows a boxy structure. Both shapes have been observed in edge-on galaxies. Seen end-on, MD again shows a boxy structure, while MH shows a spheroidal feature of considerable vertical extent. If such a feature was observed in an edge-on galaxy, it would be erroneously classified as a sizeable bulge component. Thus a number of observed bulges could in fact be bars seen end-on. A final, very important difference can be seen by comparing the relative Fourier components of the two models (AM). Model MH has a considerably larger $`m`$ = 2 component. More important, it has also considerable $`m`$ = 4, 6 and 8 components, very much in agreement with what is observed in early type barred spirals (e.g. Ohta 1996). On the contrary the Fourier components of model MD are more reminiscent of those of late type galaxies. All these differences are surprising and most of them are contrary to previous wisdom, since the stronger bar is found in the more halo dominated galaxy.

## 4. Evolution of some basic quantities

The evolution does not stop after the initial phases of bar formation. I will here describe the evolution of some basic properties, namely the strength and pattern speed of the bar and the angular momentum and mass distribution in the disc and halo components.

### 4.1. Pattern speed

Due to the interaction with the halo and the outer disc, the bar gradually slows down (see e.g. Tremaine & Weinberg 1984, Weinberg 1985, Combes et al. 1990, Little & Carlberg 1991, Hernquist & Weinberg 1992, Athanassoula 1996, Debattista & Sellwood 1998 and 2000). The rate at which this happens depends crucially on the two parameters which mainly influence the evolution, namely $`\gamma `$ and $`Q_{init}`$.

Simulations with more centrally concentrated haloes, i.e. of MH-type, show a considerably larger slowdown rate than simulations of MD-type. This is clearly seen in Fig. 2, where I plot the slowdown of the bar between times 500 and 600 as a function of $`\gamma `$. In particular for the smallest values of $`\gamma `$, i.e. 0.01 and 0.1, the slowdown is non negligible. Even for larger values, however, there is a definite trend of decreasing slowdown rate with increasing $`\gamma `$ (right panel). This is in good qualitative agreement with the results of Debattista & Sellwood (1998, 2000), who showed that in cases with massive discs there is an acceptable slowdown of the bar, contrary to simulations with massive halo components. A quantitative comparison will be given elsewhere.

The halo mass distribution, however, is not the only parameter influencing the slowdown of the bar. The left panel of Fig. 3 shows the evolution of the pattern speed with time for simulations which have initially identical live haloes and identical mass distributions in the disc component, but different $`Q_{init}`$. It is clear that $`Q_{init}`$ influences strongly the bar slowdown, in the sense that bars in initially colder discs show a much stronger slowdown rate than bars in initially hotter discs. This is further explicited in the right panel, where I give the slowdown of the bar between times 500 and 600 as a function of $`Q_{init}`$. The solid line is simply a least squares fit to the data to guide the eye and does not imply that the decrease should be linear.

Thus the slowdown rate is a function of at least two parameters, $`Q_{init}`$ and $`\gamma `$, and a particular value of this rate can be achieved not only for a single ($`Q_{init}`$, $`\gamma `$) pair, but for a sequence of such pairs. One should, therefore, refrain from using this rate to obtain limits on the values of one of the two parameters, unless the value of the other is relatively well constrained. This should be kept in mind when using arguments based on the bar slowdown in the ongoing debate about whether discs are maximum or sub-maximum, since little is known about $`Q_{init}`$. One could argue that discs were initially very cold, as they were composed mainly of gas, which formed stars with low velocity dispersion. On the other hand, it could also be argued that discs had initially very inhomogeneous mass distributions with big lumps of matter, whose interactions could heat the initial disc to relatively high temperatures. This, together with the inflow of fresh cool gas from the halo, should then determine $`Q_{init}`$.

If the $`Q_{init}`$ is sufficiently large, then the bar will not slow down substantially. Furthermore, if $`Q_{init}`$ is intermediate, then the increase in bar length (see §4.2) will compensate for the slowdown, in the sense that $`r_L/a`$ (where $`r_L`$ the corotation radius and $`a`$ the bar length) will remain roughly constant during the evolution. This will be discussed in detail elsewhere (Athanassoula 2002, in preparation). Here let us just note that, since our present knowledge does not allow us to have a clear view of the situation, we should refrain from drawing conclusions, if these are sensitive on the value of $`Q_{init}`$. Finally, the slowdown rate could also be dependent on the resolution of the simulations. However, comparisons of simulations with different resolutions should necessarily involve only simulations with the same or similar $`Q_{init}`$, if these are to give meaningful information on the effect of resolution. This issue will be further discussed elsewhere (Athanassoula 2002, in preparation).

### 4.2. Bar strength and length

The strength of the bar and its evolution with time also depend crucially on $`\gamma `$ and $`Q_{init}`$. The left panel of Fig. 4 compares the evolution of the bar strength for the three simulations already discussed in section 3. The bar strength is here defined as the total $`m`$ = 2 component of the mass, i.e. integrated over the surface of the disc, divided by the corresponding total $`m`$ = 0 component. Note that in the initial phases of the evolution the bar grows faster in simulation MD than in MH, in good agreement with the results of the 2D simulations by Athanassoula & Sellwood (1986). At later stages, however, the bar continues growing slowly but steadily in simulation MH, compared to hardly, if at all, in simulation MD. Thus, eventually, the bar in MH becomes considerably stronger than in MD.

The length of the bar in model MD does not show much evolution with time. On the contrary, for model MH it grows considerably with time.

### 4.3. Angular momentum and mass distribution

During the evolution, the disc gives angular momentum to the halo (see also Debattista & Sellwood 2000) and, as a result, the latter obtains a net rotation in the same sense as the disc. This effect is stronger for MH-type models than for MD-type ones.

As can be seen by comparing the two upper panels of Fig. 1, the disc mass radial profile evolves considerably, becoming much more centrally concentrated. This effect is stronger for initially colder discs. On the other hand the mass distribution in the halo component shows little evolution with time (Fig. 1).

## 5. Maximum versus sub-maximum discs

Observed rotation curves and velocity fields provide useful information on the gravitational forces that act in galactic discs. It is, however, not clear what fraction of this force is due to the disc and the halo components, respectively. This has been the subject of animated debates, some pieces of evidence arguing that discs are maximum, i.e. that they dominate in the inner regions, and others that their contribution is substantially lower (e.g. e.g. Athanassoula, Bosma & Papaioannou 1987; Bosma 1999, 2000 and these proceedings; Bottema 1993; Courteau & Rix 1999; Kranz, Slyz & Rix 2001; Sellwood 1999; Weiner, Sellwood & Williams 2001).

Sackett (1997) and Bosma (2000) give a simple, straightforward criterion, allowing us to distinguish maximum from sub-maximum discs. Consider the ratio $`S=V_{d,max}/V_{tot}`$, where $`V_{d,max}`$ is the circular velocity due to the disc component and $`V_{tot}`$ is the total circular velocity, both calculated at a radius equal to 2.2 disc scalelengths. According to Sackett (1997) this ratio has to be at least 0.75 for the disc to be considered maximum. Of course in the case of strongly barred galaxies the velocity field is non-axisymmetric and one should consider azimuthally averaged rotation curves, or “circular velocity” curves. Furthermore, in the case of strongly barred galaxies it is not easy to define a disc scalelength, so it is better to calculate $`S`$ at the radius at which the disc rotation curve is maximum, which is a well defined radius and is roughly equal to 2.2 disc scalelengths in the case of an axisymmetric exponential disc. After this small adjustment, we can apply the above criterion to our simulations.

As can be seen from Fig. 1, the MD disc starts off as maximum, and its importance increases further with time. On the other hand simulation MH starts as sub-maximum, with a value of $`S`$ not far from the value of 0.63 advocated by Bottema (1993). The right panel of Fig. 4 plots the evolution of $`S`$ with time and shows clearly that after roughly time 200 the disc becomes maximum, and then stays so, at a roughly constant value of $`S`$. It is interesting to note that the rather abrupt jump of the disc from sub-maximum to maximum happens right after the abrupt increase of the bar strength, i.e. right after bar formation.

In the above example the disc started as sub-maximum and ended as maximum, due to the rearrangement of the disc material. Yet this is not always the case. In cases with strong bulges, small $`\gamma `$ and/or large $`Q_{init}`$, the value of $`S`$ may not rise sufficiently for the disc to become maximum. The limiting value of $`\gamma `$ for which the disc will become maximum after bar formation depends on the values of the remaining parameters of the simulation. Interpolating from a set of simulations with $`Q_{init}`$ = 1.2, $`z_0`$ = 0.2 and $`M_h`$ = 5, we find the limiting value of $`\gamma `$ to be around 0.4. On the other hand, for initially very cold discs, with $`Q_{init}`$ = 0.1, the limiting value of $`\gamma `$ is much smaller, of the order of 0.13. This last value corresponds to a very fast rising rotation curve, reminiscent of those obtained with a cuspy halo density profile. In fact, since this value of $`\gamma `$ is only slightly larger than twice the softening length, a smaller distance would not be well resolved with the present simulations.

We can thus conclude that whether a disc ends as maximum or as sub-maximum depends on the central concentration of the halo and on the initial velocity dispersion in the disc. For highly centrally concentrated but non-cuspy haloes ($`\gamma `$ = 0.5), as those favoured by observers for early type galaxies, an initially sub-maximum disc will become maximum after the bar has formed. This holds for all such initial conditions that I tried, except if $`Q_{init}`$ is larger than 1.7. On the other hand, if the haloes are cuspy, then the threshold $`Q_{init}`$ is considerably smaller. It should also be noted that the exact position of this threshold may well depend on the initial $`z_0`$ and on the halo-to-disc mass ratio – which here has been taken to be equal to 5 – and perhaps also on numerical parameters like the softening. Thus, although the above discussion argues that it is quite likely that discs of strongly barred galaxies are maximum, a considerable amount of work is still necessary in order to produce a rigid criterion.

## 6. The role of the halo

The strongest differences in Fig. 1 are between the middle (model MH) and right (model RH) columns of panels. Indeed, the face-on view of model MH shows a very strong bar and the edge-on one a very strong peanut, or rather an ‘X’-type feature. On the contrary, model RH shows at most a slight oval distortion in the inner parts of the face-on view and not even a boxiness in the edge-on views. Similarly, the rotation curves in the second row of panels reveal that the disc material for model MH has considerably concentrated towards the center, while there is only a very small such effect for RH. Since the two models are initially identical in everything, except that the halo of RH is rigid, while that of MH is live, it must be the halo response that is responsible for the big differences.

In order to understand the halo response we need first to understand its orbital structure. We thus need to calculate, at a given time $`t`$, the principal frequencies of the halo and disc orbits, i.e. the values of the angular frequency $`\mathrm{\Omega }`$, the epicyclic frequency $`\kappa `$ and the vertical frequency $`\kappa _z`$ for each orbit (A02). For this, I calculate the total potential and the bar pattern speed at this time. In this potential, now considered as stationary in a frame of reference rotating with the measured pattern speed, I calculate the orbits of 100 000 disc and 100 000 halo particles and from these their basic frequencies. This reveales that a considerable number of halo stars are in resonance with the bar. Fig. 5 plots the number of particles (orbits), $`N_R`$, that have a frequency ratio $`R_f=(\mathrm{\Omega }\mathrm{\Omega }_p)/\kappa `$ within a bin of a given width centered on a value of this ratio, plotted as a function of $`R_f`$. The distribution is far from uniform and shows high peaks at the main resonances, both for the disc and the halo components. Disc particles at inner Lindblad resonance will emit energy and angular momentum (Lynden-Bell & Kalnajs 1972) which can be absorbed either by disc particles at corotation and outer Lindblad resonance, or by halo resonant particles. Since the bar is a negative angular momentum perturbation, the halo, by taking angular momentum from it, helps it grow. Thus haloes, instead of stabilising discs, as was previously thought, can help bars grow. This new instability mechanism (A02) can account for strong bars in cases with strong halo components, and explains why the bar in model MH is stronger than that in model MD.

However, one should not extrapolate this behaviour further to ever stronger halo components and conclude that the stronger the halo the stronger the bar should be. There will obviously be a point where the disc will be comparatively too light to emit all the angular momentum the halo can absorb. For simulations with $`M_h/M_d`$ = 5 the discs proved to be bar unstable for all the values of $`\gamma `$ I tried. This, however, is not always the case for smaller disc-to-halo mass ratios. Thus no bar formed for $`M_h/M_d`$ = 10 and $`\gamma `$ = 0.5 at least up to $`t`$ = 900, and presumably much longer. This shows that, although bars can grow in sub-maximum discs, they can not grow in galaxies with an arbitrarily small disc contribution.

## 7. Evolution in the presence of gas

The evolution of a barred galaxy containing a considerable amount of gas in the disc can differ substantially from that of a galaxy with no gaseous component.

One of the differences is that the pattern speed may stay constant, or even increase with time, contrarily to what was seen in the purely stellar cases (Friedli & Benz 1993, Heller & Shlosman 1994, Berentzen et al. 1998). Unfortunately not much quantitative information from appropriate $`N`$-body simulations can be found in the literature, and a study of the evolution of the pattern speed for different gas masses and distributions would certainly be most useful. The problem is particularly interesting since in the exchange of energy and angular momentum, discussed in the previous section, a third partner has been added to the stellar disc and the halo components, namely the gas.

As a response to the bar forcing, the gas concentrates in two very narrow and elongated regions near the leading edges of the bar, which are in fact the loci of shocks. Depending on the bar parameters, and in particular its strength, they are either straight, or curved with their concave side towards the bar major axis (Athanassoula 1992). As a result of these shocks, the gas flows inwards and accumulates in an area near the center. The size of this central area is of the order of a kpc, and is determined by the extent of the largest orbits of the $`x_2`$ family.

In order to come yet further inwards, to the nucleus of the galaxy, the gas has to shed a yet larger fraction of its angular momentum. Several possibilities have been so far put forward. In the first one (Heller & Shlosman 1994) when the central disc becomes gas dominated it is unstable and breaks into clumps. Via their collisions, these clumps may loose angular momentum and spiral to the center. A second possibility is that the $`x_2`$ family does not exist, or that it has a very small extent, in which case the central area is very small. It has, however, been argued that the density distribution in barred galaxies is such as to permit an $`x_2`$ family of sizeable extent (e.g. Teuben et al. 1986; Athanassoula 1991, 1992; Athanassoula & Bureau 2000). A similar effect can be achieved in the case of a very high sound speed, in which case the shock loci are very near the center (Englmeier & Gerhard 1997, Patsis & Athanassoula 2000). At somewhat smaller, but still quite high, sound speeds the distance between the main shock loci and the center may be bridged by a spiral, which can also transfer gas to the center (Maciejewski et al. 2002). As a final alternative let me mention nested bars, which can be either gaseous, as initially suggested by Shlosman, Frank & Begelman (1989), or stellar (Erwin, these proceedings). Recent hydrodynamical simulations, however, argue that secondary stellar bars are unlikely to increase the mass inflow rate into the galactic nucleus (Heller & Shlosman 2002, Maciejewski et al. 2002). Further work on this subject is necessary in order to elucidate the properties and the role of secondary bars.

## 8. Bar destruction and regeneration

There are at least two ways by which bars could be destroyed, or, by which their amplitude could be very severely diminished. One can be assimilated to a self-destruction, since the culprit is the central concentration formed by the gas which is pushed towards the center of the galaxy by the bar itself. In the second method the culprit is a companion or small infalling galaxy.

As mentioned in the previous section, the gas is pushed inwards by the bar. Once it has reached the center, it may create a sufficient central concentration to destroy the bar that drove it there (Friedli & Benz 1993, Berentzen et al. 1998). Indeed, orbital calculations by Hasan & Norman (1990) and by Hasan, Pfenniger & Norman (1993) have shown that a sufficiently strong central concentration can make the $`x_1`$ orbits unstable. Norman, Sellwood & Hasan (1996) showed that a massive core, of mass 5% or more of the combined disc and bulge mass, can destroy a bar. Nevertheless, the mass of black holes in disc systems seems to be smaller, by an order of magnitude or more, than what is required by Norman et al. and thus may not be sufficient to destroy bars.

In purely stellar simulations, vertical impacts by a massive companion in the bar area result in a considerable decrease of the bar amplitude, but destruction was so far found only in cases in which the disc thickened so considerably that the galaxy could not be called a disc galaxy anymore (Athanassoula 1996, 1999). On the other hand, destruction seems to be possible for such impacts if a considerable gaseous component is present (Berentzen, Athanassoula, Heller et al., in preparation). If the companion is initially in a near-circular orbit, then it can easily destroy the bar (Athanassoula 1996, 1999), even in simulations with no gaseous component.

Can a disc form a new bar, once its first one has been destroyed? In other words, can more than one bar form consecutively in the same disc? It should be noted that after the demise of the first bar, a disc is inhospitable for a new bar, since it is left with a considerable central concentration and a high velocity dispersion, both factors inhibiting bar growth. Nevertheless, fresh infalling gas may cool the disc and make it again unstable. A first attempt at making a second bar in this way (Sellwood & Moore 1999) proved successful, so further work on this subject should follow.

#### Acknowledgments.

I would like to thank A. Bosma, A. Misiriotis, M. Tagger and F. Masset for stimulating discussions and A. Misiriotis and J. C. Lambert for their collaboration on the software calculating the orbital frequencies. Part of this paper was written while I was visiting INAOE. I would like to thank the ECOS-Nord and the ANUIES for financing this trip and INAOE for their kind hospitality.

## References

Athanassoula, E. 1991, in Dynamics of Disk Galaxies, ed. B. Sundelius (Göteborg: Chalmers Univ. Technology), 149

Athanassoula, E. 1992, MNRAS, 259, 345

Athanassoula, E. 1996, in ASP Conf. Ser. Vol. 91, Barred Galaxies, eds. R. Buta, D. Crocker and B. Elmegreen, (San Francisco: ASP), 309

Athanassoula, E. 1999, in ASP Conf. Ser. Vol. 160, Astrophysical discs, eds. J. A. Sellwood and J. Goodman, (San Francisco: ASP), 160, 351

Athanassoula, E. 2002, ApJ, 569, L83 (A02)

Athanassoula, E., Bosma, A., & Papaioannou, S. 1987, A&A, 179, 40

Athanassoula, E., & Bureau, M. 2000, ApJ, 522, 699

Athanassoula, E., & Misiriotis, A. 2002, MNRAS, 330, 35 (AM)

Athanassoula, E., & Sellwood, J. A. 1986, MNRAS, 221, 213

Berentzen, I., Heller, C. H., Shlosman, I., & Fricke K. J. 1998, MNRAS, 300, 49

Bosma, A. 1999, in ASP conference series 182, Galaxy Dynamics, eds. D. R. Merritt, M. Valluri & J. A. Sellwood, 339

Bosma, A. 2000, in ASP conference Series 197, ¿From the Early Universe to the Present, eds. F. Combes, G. A. Mamon & V. Charmandaris, (San Francisco: ASP), 91

Bottema, R. 1993, A&A, 275, 16

Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A&A, 233, 82

Courteau, S. ,& Rix, H. W. 1999, ApJ, 513, 561

Debattista, V. P., & Sellwood, J. A. 1998, ApJ, 493, L5

Debattista, V. P., & Sellwood, J. A. 2000, ApJ, 543, 704

Englmaier, P., & Gerhard, O. 1997, MNRAS, 287, 57

Friedli, D., & Benz W. 1993, A&A, 268, 65

Hasan, H., & Norman, C. 1990, ApJ, 361, 69

Hasan, H., Pfenniger, D., & Norman, C. 1993, ApJ, 409, 91

Heller, C. H., & Shlosman, I. 1994, ApJ, 424, 84

Heller, C. H., & Shlosman, I. 2002, ApJ, 565, 921

Hernquist, L., & Weinberg. M. D. 1992, ApJ, 400, 80

Hohl, F. 1971, ApJ, 168, 343

Kawai, A., Fukushige, T., Makino, J., & Taiji, M. 2000, PASJ, 52, 659

Kranz, T., Slyz, A. ,& Rix, H. W. 2001, ApJ, 562, 164

Little, B., & Carlberg, R. G. 1991, MNRAS, 250, 161

Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1

Maciejewski, W., Teuben, P. J., Sparke L. S., & Stone, J. M. 2002, MNRAS, 329, 502

Miller, R. H., Prendergast, K. H., & Quirk, W. J. 1970, ApJ161, 903

Ohta, K. 1996, in ASP conference Series Vol 91, Barred Galaxies, eds. R. Buta, D. A. Crocker, B. G. Elmegreen, (San Francisco: ASP), 37

Patsis, P., & Athanassoula, E. 2000, A&A, 358, 45

Sackett, P. D. 1997, ApJ, 483, 103

Sellwood, J. A. 1999, in ASP conference series 182, Galaxy Dynamics, eds. D. R. Merritt, M. Valluri & J. A. Sellwood, (San Francisco: ASP), 351

Sellwood, J. A., & Moore E. M. 1999, ApJ, 510, 125

Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature, 338, 45

Teuben, P., Sanders, R. H., Atherton, P. D., & van Albada, G. D. 1986, MNRAS, 221, 1

Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729

Weinberg, M. D. 1985, MNRAS, 213, 451

Weiner, B. J., Sellwood, J. A., & Williams, T. B. 2001, ApJ, 546, 931


format: text

text:
# Determination of 𝛼_𝑠 from Gross-Llewellyn Smith sum rule by accounting for infrared renormalon*footnote **footnote * Partly based on invited talk given by G.C. in WG3 of NuFact02, July 1-6, 2002, Imperial College, London

##

The infrared (IR) renormalons are singularities at positive values of the Borel variable $`b`$ in the Borel transforms of QCD physical observables. These renormalons represent the contributions of loop integrations in the low energy regimes . They represent obstacles for and lead to ambiguities in the Borel integration. This integration then gives discontinuities (cuts) of the observables in the complex coupling parameter plane. It is thus important to eliminate these unphysical cuts from the Borel integration. This approach was proposed in  and applied to the GLS sum rule in . Here we recapitulate the method and the results of , and emphasize some physical aspects. The GLS sum rule is the quantity

$`M_3(Q^2)3\left(1\mathrm{\Delta }(Q^2)\right)`$ $`=`$ $`{\displaystyle _0^1}𝑑xF_3(x;Q^2)\left[1+𝒪\left({\displaystyle \frac{m_N^2}{Q^2}}\right)\right],`$ (1)
where $`F_3=(F_3^{\nu p}+F_3^{\nu n})/2`$ is the nonsinglet structure function in the $`\nu N`$ DIS, $`m_N`$ is the nucleon mass. Theoretical knowledge about the GLS $`\mathrm{\Delta }(Q^2)`$ is the following. The truncated perturbation series (TPS) up to $`a^3`$ \[$`a(Q)\alpha _s(Q)/\pi `$\] is known

$`\mathrm{\Delta }(Q^2)`$ $``$ $`a(1+w_1a+w_2a^2),`$ (2)
as is also the structure of the first IR renormalon singularity of its Borel transform

$`\mathrm{BT}(b)`$ $`=`$ $`{\displaystyle \frac{C}{(1b)^{1+\nu }}}\left[1+\kappa _1(1b)+\mathrm{}\right].`$ (3)
Here $`\nu =(\beta _1/\beta _0\gamma _2)/\beta _0`$, where $`\beta _0`$ and $`\beta _1`$ are the one- and two-loop QCD beta function coefficients and $`\gamma _2`$ is the one-loop coefficient of the anomalous dimension of the twist-four ($`d=2`$) operator appearing in the Operator Product Expansion (OPE) for $`\mathrm{\Delta }(Q^2)`$. When the number of active quark flavors is $`n_f=3`$, we have $`\beta _0=9/4`$, $`\beta _1=4`$, and $`\gamma _2=8/9`$ (the latter value obtained in ), giving $`\nu =32/81`$. The TPS of $`\mathrm{BT}(b)`$ in powers of $`b`$ is known to $`b^2`$, as a consequence of (2). A consequence of the IR singularity (3) at $`b1`$ is that the Borel integral $`\mathrm{\Delta }_{\mathrm{BI}}(z)`$, which for the complex coupling $`z=|z|\mathrm{exp}(i\varphi )`$ is defined as

$$\mathrm{\Delta }_{\mathrm{BI}}(z)=\frac{1}{\beta _0}_0^{+\mathrm{}\mathrm{exp}(i\varphi )}𝑑b\mathrm{exp}\left[\frac{b}{\beta _0z}\right]\mathrm{BT}(b),(b=|b|e^{i\varphi }),$$
(4)
has a discontinuity (cut) at the positive real axis $`z=a(Q)0`$. Namely, for $`z`$ near the positive real axis, it can be shown from (4) that

$$\mathrm{\Delta }_{\mathrm{BI}}(z)|_{z=a(Q)\pm i\epsilon }=\frac{1}{\beta _0}_{\pm i\epsilon }^{\mathrm{}\pm i\epsilon }𝑑b\mathrm{exp}\left[\frac{b}{\beta _0a(Q)}\right]\mathrm{BT}(b).$$
(5)
This quantity is not even real; the discontinuity shows up in its imaginary part

$`\mathrm{Im}\mathrm{\Delta }_{\mathrm{BI}}(a(Q)\pm i\epsilon )={\displaystyle \frac{1}{\beta _0}}\mathrm{Im}{\displaystyle _{\pm i\epsilon }^{\mathrm{}\pm i\epsilon }}𝑑b\mathrm{exp}\left[{\displaystyle \frac{b}{\beta _0a(Q)}}\right]{\displaystyle \frac{C}{(1b)^{1+\nu }}}\left[1+\kappa _1(1b)+\mathrm{}\right]`$ (6)
$`=C\beta _{0}^{}{}_{}{}^{1\nu }\left[e^{1/(\beta _0a(Q))}\right]a(Q)^\nu \mathrm{\Gamma }(\nu )\mathrm{sin}(\pi \nu )\left[1+\kappa _1\beta _0a\nu +\mathrm{}\right].`$ (7)
Because $`(ai\epsilon )^\nu =a^\nu [\mathrm{cos}(\pi \nu )\pm i\mathrm{sin}(\pi \nu )]`$, (7) leads to

$$\mathrm{\Delta }_{\mathrm{BI}}(z)=\frac{C\mathrm{\Gamma }(\nu )}{\beta _0^{1+\nu }}\left[e^{1/(\beta _0a(Q))}\right]\left[(z)^\nu +\kappa _1(\nu )\beta _0(z)^{\nu +1}+\mathrm{}\right]+\stackrel{~}{\mathrm{\Delta }}_{\mathrm{BI}}(z),$$
(8)
where $`\stackrel{~}{\mathrm{\Delta }}_{\mathrm{BI}}(z)`$ is a function without cuts in the complex $`z`$-plane since the first expression on the right-hand side absorbs the renormalon-induced cut at $`z0`$. This expression, when $`z=a(Q)`$ is real, has the $`Q`$-dependence $`Q^2\alpha _s(Q)^{\gamma _2/\beta _0}[1+𝒪(\alpha _s)]`$, the same as the twist-four term in the OPE for the GLS sum rule . The Im part of this expression, i.e. expression (7), can be identified as the Im contribution from the leading IR renormalon (3). The premise of the method is that the full first expression on the right-hand side of (8), i.e., including its real part, represents the unphysical cut contribution which is to be eliminated. This leads to the final result for the resummed GLS value

$$\mathrm{\Delta }(a(Q))=\stackrel{~}{\mathrm{\Delta }}_{\mathrm{BI}}(a(Q))=\left[\mathrm{Re}\mathrm{cot}(\pi \nu )\mathrm{Im}\right]_{\pm i\epsilon }^{\mathrm{}\pm i\epsilon }𝑑b\mathrm{exp}\left[\frac{b}{\beta _0a(Q)}\right]\frac{R(b)}{(1b)^{1+\nu }}.$$
(9)
In the Borel integration here, the exactly known IR renormalon singularity has been factored out explicitly. Function $`R(b)=(1b)^{1+\nu }\mathrm{BT}(b)`$, whose TPS up to $`b^2`$ is known exactly due to (2), has much weaker singularity at $`b=1`$ than $`\mathrm{BT}(b)`$. The function $`R(b)`$ has, in principle, singularities from other renormalons, at $`b=1,\pm 2,\mathrm{}`$ To accelerate the convergence in the TPS for $`R`$ further, we introduce in integral (9) the change of variables $`bw`$ according to the conformal mapping

$`bw(b)`$ $`=`$ $`\left(\sqrt{1+b}\sqrt{1b/2}\right)\left(\sqrt{1+b}+\sqrt{1b/2}\right)^1,`$ (10)
which maps all these renormalon locations onto the unit circle $`|w|=1`$, except for the leading IR renormalon $`w(b=1)=1/3`$ whose structure has been factored out. Similar applications of conformal mappings have been applied in other methods in Refs. . Further, expansion of $`R(b(w))`$ in powers of $`w`$ suggests that the $`\mathrm{N}^3\mathrm{LO}`$ coefficient (at $`a^4`$) in the GLS TPS (2), in $`\overline{\mathrm{MS}}`$ and with renormalization scale $`\mu =Q`$, is $`w_3=158\pm 30`$, somewhat higher than the TPS PMS prediction $`w_3130`$ . Using $`w_3=158\pm 30`$, we use for $`R(b(w))`$ in (9) the TPS up to $`w^4`$.

Expression (9) represents the massless part of the GLS sum rule. For $`Q^22`$-$`3\mathrm{GeV}^2`$, where reasonable experimental data are available, the contributions of the massive $`c`$-quarks  represent corrections contributing about $`3`$-$`4\%`$ to $`\mathrm{\Delta }(Q^2)`$. Nuclear effects of the iron target  contribute about $`1\%`$ to $`\mathrm{\Delta }(Q^2)`$ \[$`\mathrm{\Delta }_{\mathrm{Fe}}(Q^2)410^3\mathrm{GeV}^2/Q^2`$\]. The “light-by-light” contribution in the GLS TPS (2) (it is $`a^3`$) is kept apart, it contributes to $`\mathrm{\Delta }(Q^2)`$ less than $`0.5\%`$.

We work in the $`\overline{\mathrm{MS}}`$ renormalization scheme. The renormalization scale $`\mu `$ was fixed according to the principle of minimal sensitivity (PMS), obtaining $`\mu ^23.3Q^2`$. Experimentally extracted values of GLS $`\mathrm{\Delta }(Q^2)`$ at $`Q^2=3.16\mathrm{GeV}^2`$ are given by the Fermilab CCFR Collaboration

$$\mathrm{\Delta }(Q^2=3.16\mathrm{GeV}^2)=0.149\pm 0.039,$$
(11)
where the nuclear corrections were subtracted out. When we confront (11) with the theoretical GLS expression \[Eq. (9) with the conformal mapping (10); plus $`c`$-quark corrections; plus “light-by-light;”\], the obtained QCD $`\alpha _s`$ in $`\overline{\mathrm{MS}}`$ is

$$\alpha _s(Q)=0.305_{0.074}^{+0.132}(\mathrm{exp})\pm 0.006(\mathrm{th}),\alpha _\mathrm{s}(\mathrm{M}_\mathrm{Z})=0.1167_{0.0118}^{+0.0128}(\mathrm{exp})\pm 0.0008(\mathrm{th}).$$
(12)
Using the central CCFR data value of $`\mathrm{\Delta }(Q^2)`$ at $`Q^2=2.00\mathrm{GeV}^2`$, we obtain almost the same central value: $`\alpha _s(M_Z)=0.1166`$. The values of $`\alpha _s(M_Z)`$ obtained are in agreement with the present world average. The experimental uncertainties are large. The theoretical uncertainties are dominated by the ambiguity with respect to the change of the renormalization scheme.

The term proportional to $`\mathrm{cot}(\pi \nu )`$ in (9), which eliminates the unphysical cut, was important numerically as it decreases the value of $`\mathrm{\Delta }(Q^2)`$ by about $`10\%`$, thus increasing the extracted value of $`\alpha _s(M_Z)`$ by about $`0.003`$. The CCFR Collaboration , by evaluating the NNLO TPS of $`\mathrm{\Delta }(Q^2)`$, corrected by a term $`1/Q^2`$ as estimated by QCD sum rules, obtained from the central values of the data $`\alpha _s(M_Z)0.113`$-$`0.114`$. An earlier similar (TPS-PMS) evaluation  gave $`\alpha _s(M_Z)0.115`$.

###### Acknowledgements.

The work of C.C. and G.C. was supported by DGIP (UTFSM) and FONDECYT (Chile) Grant No. 1010094, respectively. The work of T. L. and K.S.J. was supported in part by the BK21 Core Program.


text:
# New insights in high-resolution spectroscopy: a wide theoretical library of R=500 000 stellar spectra

## 1. Introduction

New-generation spectrographs, at the major ground-based telescopes, have begun to pour in the hard-disks of astronomers’ computers an increasing mass of high-quality spectroscopic data. Inverse resolutions as high as $`\lambda /\mathrm{\Delta }\lambda =\mathrm{100\hspace{0.17em}000}`$ can now be easily attained, at least for the brightest ($`V<15`$ mag) objects in the sky, and this is pushing the observation of local and extragalactic stellar systems to an ever unimagined resolution level. The outstanding performance of instruments like UVES and VIRMOS at the ESO Very Large Telescope, or HIRES and SARG at the Keck Observatory and Telescopio Nazionale Galileo, respectively, urges therefore theoretical tools of comparable accuracy level in order to consistently match and analyse such a huge amount of observational data.

In this framework, and to help filling the gap, we undertook a long-term project aimed at providing the community with a systematic theoretical library of high-resolution stellar spectra (virtually the largest sample currently available in the literature) in the optical range ($`\lambda =35007000`$ Å) and at an inverse resolution of $`R=\mathrm{500\hspace{0.17em}000}`$.

## 2. The optical grid of synthetic spectra

This first set of spectral energy distributions relied on the SYNTHE code of Kurucz (1993), using ATLAS 9 model atmospheres (Kurucz 1995) as input. Computations have been carried out at the Brera Observatory in Milan (Italy). The input models follow the classical approximations of steady-state, homogeneous, LTE, plane-parallel layers, with a microturbulence velocity of 2 km/s and a mixing-length value $`\mathrm{}/H_p=1.25`$. They make use of an improved treatment of convective overshooting in the calculation of the transfer equations (Castelli, Gratton, & Kurucz 1997).

Over 46 million absorption lines are accounted for in the code, including all atomic elements at different ionization states and the most important diatomic molecules. All line data were extracted from the Kurucz’ (1992) original database, with the major improvement of the TiO contribution, for which we adopted the new list of lines computed by Schwenke (1998).

The wavelength interval spans the whole optical range, from the Balmer break to H$`\alpha `$. All the passbands of the Lick/IDS spectrophotometric indices (Worthey et al. 1994) are included, as well as the high-resolution index set defined by Rose (1994). Spectra are sampled at variable wavelength step maintaining a constant resolution ($`\mathrm{\Delta }\lambda /\lambda =2\times 10^6`$).

The library extends to effective temperatures as hot as $`T_{\mathrm{eff}}=\mathrm{50\hspace{0.17em}000}`$ K, and provides a suitable match to O $``$ K spectral types. M stars, cooler than $`T_{\mathrm{eff}}<4000`$ K, are however still lacking in our dataset given the missing contribution of tri-atomic molecular opacity in the input model atmospheres (Kurucz 1992).

Gravity spans a wide range across the H-R diagram ($`5\mathrm{log}g0`$), namely from dwarf (MK class V) to supergiant stars (MK class I), while metallicity accounts for stars at the two extreme edges of the Galaxy \[Fe/H\] distribution with $`3.0<[\mathrm{M}/\mathrm{H}]<+0.3`$ as boundary limits. A total of 832 synthetic spectral energy distributions (SEDs) have been computed, as summarized in Fig. 1 and Table 1. An illustrative example of the model output, sampling the spectral emission around the atomic Mg $`b`$ triplet, is shown in Fig. 2 and 3, with varying physical parameters.

A complementary project, which will extend our analysis to the ultraviolet interval, between 850–4750 Å with an inverse resolution $`R=\mathrm{50\hspace{0.17em}000}`$, is in progress at INAOE and will be completed soon (Rodríguez-Merino 2002). The extended dataset consists of over 1000 SEDs covering the effective temperature range between 3500 and 50 000 K at surface gravity $`\mathrm{log}g`$ = 1.0 $``$ 5.0 and chemical compositions $`[\mathrm{M}/\mathrm{H}]=+0.5,\mathrm{\hspace{0.17em}0.0},0.5`$ and $`1.5`$ (Rodríguez-Merino et al. 2001). In Fig. 4 we display a subset of spectra computed for solar chemical composition, $`\mathrm{log}g=5`$ and different effective temperatures.

## 3. The solar spectrum validation

A comparison with the solar spectrum represents a first mandatory test in order to assess the degree of reliability of our theoretical library. We present in Fig. 5 the match of the Solar Flux Atlas by Kurucz et al. (1984) and the corresponding synthetic solar spectrum, which we computed at similar resolution (i.e. $`R=\mathrm{522\hspace{0.17em}000}`$). The input model atmosphere is from the ATLAS 9 and assumes $`T_{\mathrm{eff}}=5777`$ K, $`\mathrm{log}g=4.4377`$, a microturbulent velocity of 1.0 km/s, and chemical abundance according to Anders & Grevesse (1989).

The figure is made up by a dense sequence of positive and negative spikes, which mainly track the residual flux differences in the absorption line profiles. While positive spikes mainly originate from observed lines that are not present in the synthetic model, the negative peaks come from theoretical absorption lines that are too weak or undetected in the observed spectrum. Note, as a stricking feature in the plot, that both positive and negative spikes tend to decrease at longer wavelength, and no systematic drift is present in the data distribution. Fig. 5 gives also a direct measure of the intrinsic uncertainty in the input physics of the models; excluding the telluric bands, clearly affecting the residual distribution about 6300 and 6900 Å, the rms value in the plot is $`\pm 0.086`$.

## References

Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197

Castelli, F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841

Kurucz, R. L. 1992, in IAU Symp. 149, The Stellar Populations of Galaxies, ed. B. Barbuy & A. Renzini, (Dordrecht: Kluwer), 225

Kurucz, R. L. 1993, CD-ROM No. 18, SYNTHE spectrum synthesis programs and line data

Kurucz, R. L. 1995, CD-ROM No. 13, ATLAS9 Stellar Atmosphere Programs and 2 km/s Grid, revised

Kurucz, R. L., Furenlid, I., Brault, L., & Testerman, L. 1984, National Solar Observatory Atlas, Sunspot, New Mexico: National Solar Observatory

Rodríguez-Merino, L.H. 2002, Ph. D. Thesis, in preparation

Rodríguez-Merino, L. H., Chávez, M., Buzzoni, A., Bertone, E. 2001, in proc. of New Quests in Stellar Astrophysics: A Link between Stars and Cosmology, ed. M. Chavez, A. Bressan, A. Buzzoni & D. Mayya, (Dordrecht: Kluwer), 39

Rose, J. A. 1994, AJ, 107, 206

Schwenke, D. W. 1998, Faraday Discuss., 109, 321

Worthey, G., Faber, S. M., Gonzáles, J. J., & Burnstein, D. 1994, ApJS, 94, 687


text:
# Precision 𝑊-pair physics with the YFSWW3 and KoralW Monte Carlos

## 1 Introduction

In this note we give a short summary of our recent results and developments in the area of precision $`W`$-pair physics performed with the help of the Monte Carlo (MC) event generators KoralW  and YFSWW3 . We will concentrate mostly on the LEP2 physics; however, some results for the energy range of future linear colliders will be also discussed.

The LEP2 data analysis is gradually approaching its final, highest precision stage . The forecast accuracy of 1% for the W-pair total cross section, 30 MeV for $`W`$ mass or 0.01 for anomalous coupling $`\lambda `$ as an example of anomalous vector-boson triple-gauge couplings (TGCs), will be an impressive achievement. Also in the area of single-$`W`$ production the experimental precision for the total cross section is expected to be as high as 7%.

Such a good experimental performance immediately puts demands on the theoretical calculations – they should be accurate to at least half of the experimental precision. All kinds of theoretical corrections have been extensively studied throughout the last years and we will not attempt to cover this vast subject here. Any interested reader will find a review in the Reports . We will only mention here that some of the most difficult corrections turned out to be the $`𝒪(\alpha )`$ electroweak (EW) ones. These corrections have been included in the YFSWW3 code in the so-called “double-pole approximation” based on the expansion in the $`\mathrm{\Gamma }_W/M_W`$ parameter and keeping only the highest (doubly) resonant contributions, all done in a gauge-invariant way . Apart from YFSWW3, only one other MC code, the RacoonWW , is able to calculate these corrections.

The available studies in the literature of theoretical precision have been restricted to the total cross section , lacking any discussion on the other basic observables such as $`M_W`$ or TGCs. The attempts to naively extrapolate the error of $`\sigma _{WW}`$ to the other observables (in the form of “error bands” for example) have failed, as other observables do not depend so much on the overall normalization, but rather on the shape of some differential distributions. For that reason, we decided to study both $`M_W`$  and TGCs  in detail. The TGC study has been done in collaboration with R. Brunelière, A. Denner, S. Dittmaier, S. Jézéquel, M. Roth and D. Wackeroth.

We have also done a number of upgrades and extensions of our MC codes. Among them the most important ones are: (1) the “merge without merging” of KoralW and YFSWW3 into a one MC code on an event-per-event basis, with the help of the novel technique of Concurrent Monte Carlo , and (2) the improvements in KoralW in the domain of single-$`W`$ processes , necessary to meet the experimental precision in this area.

## 2 KoralW&YFSWW3 – the Concurrent Monte Carlo for $`W`$ physics

We have developed two Monte Carlo codes for the $`W`$-pair physics: KoralW and YFSWW3. The first of them, KoralW, includes the complete, fully massive matrix element for all four-fermion processes, generated by the GRACE package  along with two independent and efficient multi-channel algorithms for the four-fermion phase-space integration . The YFSWW3 includes the complete $`𝒪(\alpha )`$ electroweak corrections to the $`W`$-pair production process in the double pole approximation, based on the on-shell calculations of Refs. . Apart from that, the two programs share a number of features such as the initial state radiation (ISR) in Yennie–Frautschi–Suura exponentiation up to $`𝒪(\alpha ^3)`$ LL with multiple photons with non-zero transverse momenta, Coulomb correction or anomalous TGCs. Also the semi-analytical program KorWan for $`W`$-pair production and decay is included in both packages.

To meet the high precision requirements of the LEP experiments the corrections of both codes need to be available simultaneously, on an event-per-event basis. The most obvious way of achieving it, by merging the two codes at the level of source codes, is highly impractical. The already big codes (few hundred thousand lines) would become even larger; it would not be possible to maintain them separately, or the merging would have to be repeated after each new release. Therefore we proposed a different approach: the Concurrent Monte Carlo (CMC) . In CMC both codes remain completely independent, but they run simultaneously and communicate in real time through the UNIX/LINUX facility called “named pipes” or “FIFO files”. From the technical point of view, there is no need to modify any of the codes (apart from the changes necessary to accommodate the physics requirements) as the FIFO files behave just like ordinary ASCII files. A sample of such a setup is shown in Fig. 1,

where we present the main operation mode of CMC KoralW&YFSWW3 in which the events are constructed with full background four-fermion corrections and the $`𝒪(\alpha )`$ EW corrections.

The concept of CMC has a number of advantages: (1) it provides a natural mechanism of “encapsulation” for any codes written in any languages without any need to modify their source codes; (2) it has a “lego block” structure with many possible arrangements of involved codes (for example for reprocessing data by two KoralW codes running concurrently but initialized with different input parameters ); (3) it is naturally suited to use multiprocessor computers efficiently.

## 3 Theoretical precision of <br>$`M_W`$ measurement

How did we proceed with the estimation of the theoretical error of the $`M_W`$ measurement? We set-up a simplified framework based on one-dimensional fits to the invariant mass of one of the $`W`$’s. These fits were performed with the help of the semi-analytical code KorWan , based on the CC03 matrix element convoluted with the ISR and with the screened Coulomb correction, both matching exactly the corresponding corrections in KoralW and YFSWW3 MCs. In the first step, we checked that the fitting procedure works properly and we reproduce in the fit the value of the $`W`$ mass used in the MC simulation. This introductory exercise shows also that if one is interested only in the relative shifts of $`M_W`$ due to various effects, then any of the fitting functions (FFs) that we used is equally good.

In the next step we switched on and off various corrections in the MC sample and in FFs to see how they influence the $`M_W`$. We found out that the most complicated corrections: four-fermion background and EW non-leading $`𝒪(\alpha )`$ ones are both of the size of 1 MeV.

In the third step, we compared the results of the fits to the KoralW&YFSWW3 data with the fits to the data generated by the RacoonWW MC code, taken from Ref. . Both these codes include $`𝒪(\alpha )`$ EW and four-fermion corrections, but they differ in many aspects of the actual implementation of these corrections. Therefore such a comparison provides in particular a good test on the technical precision of the calculations. The fitted value of $`M_W`$ obtained from corresponding data from these two codes turned out to differ by less than 3 MeV.

Finally, we estimated the size of the possible missing corrections in KoralW&YFSWW3. This has been done for missing parts of ISR ($`<1`$ MeV), FSR ($`2`$ MeV) as well as non-factorizable ($`<2`$ MeV) ones.

All the above work brought us to the following conclusion: the electroweak theoretical precision of $`M_W`$ reconstruction through KoralW&YFSWW3 at LEP2 energies is $``$ 5 MeV. This estimate includes a “safety factor” of 2, due to the simplified 1-dimensional framework that we used. The above estimate is in good agreement with the small differences between the YFSWW3 and RacoonWW MCs. This value is well below the expected final LEP2 experimental precision of $``$ 30 MeV.

## 4 Theoretical precision of anomalous $`\lambda `$ coupling measurement

This study of $`\lambda `$ coupling has been done in collaboration with R. Brunelière, A. Denner, S. Dittmaier, S. Jézéquel, M. Roth and D. Wackeroth.

The measurement of anomalous TGCs is primarily based on the distribution of the cosine of the $`W`$ boson production angle $`\theta _W`$. The most general parametrization of the anomalous TGCs requires fourteen complex parameters. In our analysis we concentrated on the simplest scenario of only one anomalous TGC, the $`\lambda =\lambda _\gamma =\lambda _Z`$. We adopted essentially the same strategy as for the $`M_W`$ analysis, i.e. we performed 1-dimensional fits to the $`\mathrm{cos}\theta _W`$ distributions.

As in the $`M_W`$ case, we began by a simple fit of the semi-analytical FF from KorWan. It turned out, however, that the size of the various effects obtained from such a fit differ significantly, depending on the acceptances applied to the MC data. Therefore we had to find a way of incorporating more realistic experimental acceptances into FFs. We did this by using the Monte Carlo Parametric Fitting Function (MPFF), constructed directly as a parametrization of MC data, in the form of a polynomial in $`\lambda `$ and $`\mathrm{cos}\theta _W`$. In particular we used the real-life ALEPH acceptances. The performed fits to MPFFs thus show that the $`𝒪(\alpha )`$ NL EW corrections lead to a $`\lambda `$ shift of about 0.01 and that the differences between KoralW&YFSWW3 and RacoonWW are very small, below 0.001. In the next step we also looked at what happens when we extend our analysis to the different channels, other than $`\mu \nu _\mu qq`$. This has been done with the help of “transfer matrices”, constructed on the basis of the full ALEPH detector simulation. These matrices have been applied to our MC $`\mathrm{cos}\theta _W`$ distributions. As a result we showed that the effect of $`𝒪(\alpha )`$ NL EW corrections on $`\lambda `$ is similar (within 15%) in all channels, including the difficult four-quark and tau-semileptonic ones. In the last step we have estimated the size of missing corrections in both MC codes.

As a result we came to the conclusion that the EW theoretical uncertainty in $`\lambda `$ of the KoralW&YFSWW3 and RacoonWW MC codes is $`0.005`$. As before this number includes a “safety factor” of 2 due to the 1-dimensional fit and missing higher-order effects.

## 5 Improvements of KoralW for the single-$`W`$ processes

Apart from the $`W`$-pair production, another interesting process involving the $`W`$ boson is the single-$`W`$ process. It is usually defined as the $`e\overline{e}e\overline{\nu }_ef^{}\overline{f}`$ with outgoing $`e`$ close to the beam line (say, $`\mathrm{cos}\theta _e0.95`$) . In such a configuration the $`t`$-channel graphs dominate. As a consequence, the photonic radiation based on initial state only, sufficient for $`W`$-pair production, has to be corrected to include the screening of the electric charge (ECS) due to the dipole-like structure formed by initial- and final-state electrons. This screening is responsible for the strong suppression of photons with large transverse momenta and for the change of scale of the “leading logarithms” from $`s`$ to $`|t|`$.

The KoralW code, based on ISR, has missed the ECS effect. Therefore we have proposed  a simple way of reintroducing it by just adding an additional weight, without any modification of the original code. The weight results from an educated guess rather than the rigorous derivation from the first principles. It is constructed as a ratio of one-photon differential distributions of a coherent sum of IS and FS radiation over an analogous incoherent IS and FS sum. These distributions are exact in the soft photon limit and have the simple form of $`(p_a/kp_ap_b/kp_b)^2`$. In order to maintain the independence of the program of the dummy infra-red cut-off $`ϵ_{IR}`$ we also introduced the corresponding virtual weight. As in the real-emission case, we have constructed it as an ansatz that fulfills the requirement of cancelling the spurious $`ϵ_{IR}`$ dependence. The overall precision of this ECS ansatz we estimate at 2%, because of some missing nonleading terms.

In order to meet the experimental requirements of LEP2, we had to introduce also the effect of running the QED coupling constant in the $`e`$ vertex from the $`M_W^2`$ to $`t`$ scale. We have done this in the form of an overall factor $`(\alpha (t)/\alpha _{G_\mu })^2`$. Such a naive ansatz has a precision better than 2%, see Ref. .

In summary, we estimate at 3% the overall precision of the proposed improvements in KoralW in the single-$`W`$ region of phase space. This is sufficient for the expected LEP2 experimental precision of about 7%.

We present in Table 1 some numerical results for the single-$`W`$ total cross sections in the $`e^+e^{}e^{}\overline{\nu }_eu\overline{d}`$ process with the following cuts: $`\mathrm{cos}\theta _e0.997`$, $`M_{q\overline{q}}45`$ GeV. As can be seen, the effects of ECS and running of QED coupling cancel each other in the total cross section in a wide energy range, from the LEP2 up to future linear colliders area. This cancellation is however purely accidental and would not be true for differential distributions.

## 6 Summary

In this short note we presented recent results and developments in precision $`W`$-pair physics related to the KoralW and YFSWW3 MC programs:

(1) We established that the theoretical precisions, so far omitted in the literature, of fundamental LEP2 observables other than total cross section are: 5 MeV for $`M_W`$ (LEP2 experimental precision $``$ 30 MeV) and 0.005 for anomalous $`\lambda `$ coupling (LEP2 experimental precision $``$ 0.01). These precisions fully meet the LEP2 experimental requirements.

(2) We proposed a new way of combining independent MC codes on an event-per-event basis, by running them in parallel and providing the real-time data exchange between them by means of the UNIX/LINUX facility of “named pipes”. We used this scheme to create a Concurrent MC KoralW&YFSWW3 that generates $`W`$-pair events with both four-fermion background and $`𝒪(\alpha )`$ corrections.

(3) We presented the upgrade of KoralW in the kinematical domain of single-$`W`$ targeted at better emulation of the $`t`$-channel-dominated photonic radiation. Thanks to this upgrade, the overall precision of KoralW for this type of processes at LEP2 has become 3% (LEP2 experimental precision $``$ 7%).

## Acknowledgements

This work has been supported in part by the Polish Government grants KBN 5P03B09320, 2P03B00122, the European Commission 5-th framework contract HPRN-CT-2000-00149, the NATO Grant PST.CLG.977751, the US DoE contract DE-FG05-91ER40627 and the Polish-French Collaboration within IN2P3 through LAPP Annecy.


text:
# TEST OF GENERAL RELATIVITY: 1995-2002 MEASUREMENT OF FRAME-DRAGGING

## 1 Gravitomagnetic phenomena on test gyroscopes, test particles, clocks and photons

In this paper we review some phenomena arising in the vicinities of a rotating body and some proposals and recent measurements of frame-dragging and Lense-Thirring effect obtained by laser-ranged satellites. For a general review of tests and measurements of general relativistic and gravitational effects we refer to .

Einstein’s general theory of relativity  predicts the occurrence of peculiar phenomena on test gyroscopes, test particles, clocks and photons in the vicinity of a mass-energy current and thus in the vicinity of a spinning body due to its rotation.

In section (2) we describe how the orbit of a test-particle is influenced by the spin of a central body and the direct measurement of the orbital perturbations of laser ranged satellites due to the Earth spin, i.e. the Lense-Thirring effect . In 1995-2002 the Lense-Thirring effect was measured with about 20$`\%`$ accuracy using the LAGEOS and LAGEOS II satellites , see section (2).

Small test-gyroscopes, that determine the axes of a local, freely falling, inertial frame, where “locally” the gravitational field is “unobservable”, rotate with respect to “distant stars” due to the spin of a body. This effect is described in  and should be measured with about 1$`\%`$ accuracy by the Gravity Probe-B experiment .

However, not only test particles and gyroscopes are affected by the spin of the central object but also photons and clocks. In this section we review some gravitomagnetic phenomena on clocks and photons. A photon co-rotating around a spinning body takes less time to return to a “fixed point” (with respect to distant stars) than a photon rotating in the opposite direction . Since light rays are used to synchronize clocks, the different travel-time of co-rotating and counter-rotating photons implies the impossibility of synchronization of clocks all around a closed path around a spinning body; the behavior of light rays and the behavior of clocks around a spinning body are intimately connected. In several papers the “frame-dragging clock effect” around a spinning body has been estimated and some space experiments have been proposed to test it . Thus, when a clock, co-rotating very slowly (using rockets) around a spinning body and at a constant distance from it, returns to its starting point, it finds itself advanced relative to a clock kept there at “rest” (with respect to “distant stars”, see above). Similarly a clock, counter-rotating arbitrarily slowly and at a constant distance around the spinning body, finds itself retarded relative to the clock at rest at its starting point . For example, when a clock that co-rotates very slowly around the spinning Earth, at $`r`$ 6000 km altitude, returns to its starting point, it finds itself advanced relative to a clock kept there at “rest” (with respect to “distant stars”) by $`\mathrm{\Delta }\tau \frac{4\pi J_{}}{r}5\times 10^{17}`$ $`sec`$, where $`J_{}\mathrm{\hspace{0.17em}145}cm^2`$ is the Earth angular momentum. Similarly, a clock, that counter-rotates very slowly around the spinning Earth, finds itself retarded relative to a clock kept there at “rest” by the same amount. Then, the difference between the time read by the two clocks when they meet again after a whole revolution is about $`10^{16}`$ .

However, Einstein’s gravitational theory predicts peculiar phenomena also inside a rotating shell . In reference  we derive the time-delay in travel-time of photons due to the spin of a body both outside a rotating body and inside a rotating shell. We then show that this time-delay by the spin of an astrophysical object might be detected in different images of the same source by gravitational lensing.

Since here we are only interested to analyze the time delay due to spin, we chose a simple configuration where source, lens and observer are aligned and we use quasi-Cartesian coordinates . We then get, for a photon with impact parameter $`b`$ traveling on the equatorial plane of the source: $`\mathrm{\Delta }T=4M\mathrm{ln}\left(2\frac{\overline{z}}{b}\right)+\frac{4J}{b}`$. In this expression $`\overline{z}b`$ is the distance of source and observer from the lens, the first term is the standard Shapiro time delay and the second term is the gravitomagnetic time-delay due to the spin of the deflecting body.

Let us now give the order of magnitude of the time delay due to the spin of some astrophysical sources. For the sun, by considering two light rays on the equatorial plane of the Sun, grazing the Sun on opposite sides, the relative gravitomagnetic time delay is $`\mathrm{\Delta }T_{rel}^J=3.3510^{12}`$ sec. For the lensing galaxy of the Einstein cross, by assuming a simple model for rotation and shape of the central object (see  and references therein), we then get: $`\mathrm{\Delta }T_{12}^J8hr`$. As a third example we consider the relative time delay of photons due to the spin of a typical cluster of galaxies. Depending on the geometry of the system and on the path followed by the photons, we then find relative time delays ranging from a few minutes to several days. Then, at least in principle, one could detect the time delay due to the spin of a lensing galaxy by removing the larger quadrupole-moment time delay by a method described in ; of course, as in the case of the Sun, one should be able to accurately enough model and remove all the other delays, due to other physical effects, from the observed time delays between the images.

Let us now analyze the time delay in the travel time of photons propagating inside a shell of mass $`M`$ rotating with angular velocity $`\omega `$. Inside a spinning shell it is in general not possible to synchronize clocks all around a closed path. Indeed, if we consider a clock co-rotating very slowly on the equatorial plane along a circular path with radius $`r`$, when back to its starting point it is advanced with respect to a clock kept there at rest (in respect to distant stars). The difference between the time read by the co-rotating clock and the clock at rest is equal to: $`\delta T=\frac{8M}{3R}\pi \omega r^2`$. For a shell with finite thickness we just integrate this expression from the smaller radius to the larger one.

If we now consider a photon traveling with an impact parameter $`r`$ on the equatorial plane of a galaxy; the time delay due to the rotation of the external mass for every infinitesimal shell with mass $`dm=4\pi \rho R^2dR^{}`$ and radius $`R^{}\left|r\right|`$, is : $`\mathrm{\Delta }t_{dm}=\frac{8dm}{3}\omega r\frac{\sqrt{R^2r^2}}{R^{}}`$. This is the time delay due to the spin of the external thin shell. By integrating this expression from $`\left|r\right|`$ to the external shell radius $`R`$, we have:

$$\mathrm{\Delta }T=\frac{32\pi }{3}\omega r_{\left|r\right|}^R\rho R^{}\sqrt{R^2r^2}𝑑R^{}$$
(1)
This is the time delay due to the spin of the whole rotating mass of the external shell. From this formula we can easily calculate the relative time delay between two photons traveling on the equatorial plane of a rotating shell, with impact parameters $`r_1`$ and $`r_2`$.

Let us give the time delay corresponding to some astrophysical configurations. In the case of the ”Einstein cross”, after some calculations, based on a standard model for the lensing galaxy (see ref.  and references therein), the order of magnitude of the relative time delay of two photons traveling at a distance of $`r_1650pc`$ and $`r_2650pc`$ from the center, using (1) in the case $`r_1r_2`$, is: $`\mathrm{\Delta }T20`$ min. If the lensing galaxy is inside a rotating cluster, or super-cluster, to get an order of magnitude of the time delay due to the spin of the mass rotating around the deflecting galaxy, we use typical super-cluster parameters (see  and references therein). If the galaxy is in the center of the cluster and light rays have impact parameters $`r_115kpc`$ and $`r_215kpc`$ (of the order of the Milky Way radius), the time delay, applying formula (1) in the case $`r_1r_2`$ and $`\rho =constant`$, is: $`\mathrm{\Delta }t1day`$.

Finally, if the lensing galaxy is not in the center of the cluster but at a distance $`r=aR`$ from the center, with $`0a1`$ and $`R`$ radius of the cluster, by integrating (1) between $`r=aR`$ and $`R`$, when $`r_1r_2`$ we have:

$$\mathrm{\Delta }T=\frac{32\pi }{9}\omega (r_1r_2)\rho (1a^2)^{1/2}(14a^2)R^3$$
(2)
Thus, if the lensing galaxy is at a distance of 10 Mpc from the center of the cluster, the relative time delay due to the spin of the external rotating mass between two photons with $`(r_1r_2)30kpc`$, is: $`\mathrm{\Delta }T0.9day`$. Since the present measurement uncertainty in the lensing time delay is of the order of 0.5 day , the spin time delay might already be observable.

In conclusion, in ref.  we have derived and studied the ”spin-time-delay” in the travel time of photons propagating near a rotating body, or inside a rotating shell due to the angular momentum. We found that there may be an appreciable time delay due to the spin of the body, or shell, thus spin-time-delay must be taken into account in the modeling of relative time delays of the images of a source observed at a far point by gravitational lensing. This effect is due to the propagation of the photons in opposite directions with respect to the direction of the spin of the body, or shell. If other time-delays can be accurately enough modeled and removed from the observations , one could directly measure the spin-time-delay due to the gravitomagnetic field of the lensing body. We have analyzed the relative time-delay in the gravitational lensing images caused by a typical rotating galaxy, or cluster of galaxies. We have then analyzed the relative spin-time-delay when the path of photons is inside a galaxy, a cluster, or super-cluster of galaxies rotating around the deflecting body; this effect should be large enough to be detected at Earth. The measurement of the spin-time-delay, due to the angular momentum of the external massive rotating shell, might be a further observable for the determination of the total mass-energy of the external body, i.e. of the dark matter of galaxies, clusters and super-clusters of galaxies. Indeed, by measuring the spin-time-delay one can determine the total angular momentum of the rotating body and thus, by estimating the contribution of the visible part, one can determine its dark-matter contribution. The estimates presented in this paper are preliminary because we need to apply the spin-time-delay to some particular, known, gravitational-lensing images. Furthermore, we need to estimate the size and the possibility of modeling other sources of time-delay. Nevertheless, we conclude that, depending on the geometry of the astrophysical system considered, the relative spin-time-delay can be a quite large effect.

## 2 Measurement of gravitomagnetism with laser ranged satellites

In ref.  (see also ref. ) we describe the LARES experiment to measure the Lense-Thirring effect with relative accuracy of about 1$`\%`$. This laser ranged satellite, by detecting its perigee rate, would also test the foundations and other basic phenomena of general relativity and gravitational interaction. Indeed, LARES would improve the bounds on hypothetical long-range gravitational forces and the bounds on deviations from the inverse square law for very weak-field gravity; LARES would improve, by about two orders of magnitude, the accuracy in testing the equivalence principle and would provide an improved measurement in the field of Earth of the PPN (Parametrized-Post-Newtonian) parameters $`\alpha _1,\beta `$ and $`\gamma `$ .

In this section we describe the 1995-2002 measurements of the Lense-Thirring effect obtained using LAGEOS and LAGEOS II.

The measurement of distances has always been a fundamental issue in astronomy, engineering, and science in general. So far, laser ranging has been the most accurate technique for measuring the distances to the moon and artificial satellites. Short-duration laser pulses are emitted from lasers on Earth, aimed at the target through a telescope, and then reflected back by optical cube-corner retroreflectors on the moon or an artificial satellite , such as LAGEOS. By measuring the total round-trip travel time, one can determine the distance to a retroreflector on the moon with an accuracy of about 2 cm and to the LAGEOS satellites with a few millimeters accuracy.

Our detection and measurement of the Lense-Thirring effect was obtained by using the satellite laser ranging data of the satellites LAGEOS (LAser GEOdynamics Satellite, NASA) and LAGEOS II (NASA and ASI, the Italian Space Agency) and the Earth gravitational field models, JGM-3 and EGM-96. The LAGEOS satellites are heavy brass and aluminum satellites, of about 406 kg weight, completely passive and covered with retroreflectors, orbiting at an altitude of about 6,000 km above the surface of Earth. LAGEOS, launched in 1976 by NASA, and LAGEOS II, launched by NASA and ASI in 1992, have an essentially identical structure but they have different orbits. The semimajor axis of LAGEOS is $`a\mathrm{\hspace{0.17em}12},270`$ km, the period $`P\mathrm{\hspace{0.17em}3.758}`$ hr, the eccentricity $`e\mathrm{\hspace{0.17em}0.004}`$, and the inclination $`I\mathrm{\hspace{0.17em}109.9}^{}`$. The semimajor axis of LAGEOS II is $`a_{II}\mathrm{\hspace{0.17em}12},163`$ km, the eccentricity $`e_{II}\mathrm{\hspace{0.17em}0.014}`$, and the inclination $`I_{II}\mathrm{\hspace{0.17em}52.65}^{}`$.

We analyzed the laser-ranging data adopting the IERS (International Earth Rotation Service) conventions  in our modeling, however, in the 1998 analysis, we used the static and tidal EGM-96 model . Error analysis of the LAGEOS orbits indicated that the EGM-96 errors can only contribute periodic root-sum-square errors of 2 to 4 mm radially, and in all three directions they do not exceed 10 to 17 mm. The initial positions and velocities of the LAGEOS satellites were adjusted for each 15-day batch of data, along with variations in their reflectivities. Solar radiation pressure, Earth albedo, and anisotropic thermal effects were also modeled . In modeling the thermal effects, the orientation of the satellite spin axis was obtained from ref. . Lunar, solar, and planetary perturbations were also included in the equations of motion, formulated according to Einstein’s general theory of relativity with the exception of the Lense-Thirring effect, which was purposely set to zero. All of the tracking station coordinates were adjusted (accounting for tectonic motions) except for those defining the TRF terrestrial reference frame. Polar motion was also adjusted, and Earth’s rotation was modeled from the very long baseline interferometry-based series SPACE96 . We analyzed the orbits of the LAGEOS satellites using the orbital analysis and data reduction software GEODYN II .

The node and perigee of LAGEOS and LAGEOS II are dragged by the Earth’s angular momentum. From the Lense-Thirring formula , we get $`\dot{\mathrm{\Omega }}_I^{LenseThirring}`$ $`31`$ mas/yr and $`\dot{\mathrm{\Omega }}_{II}^{LenseThirring}31.5`$, where mas is a millisecond of arc. The argument of pericenter (perigee in our analysis), $`\omega `$, also has a Lense-Thirring drag , thus, we get for LAGEOS: $`\dot{\omega }_I^{LenseThirring}`$ $`32`$ mas/yr, and for LAGEOS II: $`\dot{\omega }_{II}^{LenseThirring}`$ $`57`$ mas/yr . The nodal precessions of LAGEOS and LAGEOS II can be determined with an accuracy of the order of 1 mas/year. Over our total observational period of about 4 years, we obtained a root mean square (RMS) of the node residuals of about 4 mas for LAGEOS and about 7 mas for LAGEOS II . For the perigee, the observable quantity is the product $`ea\dot{\omega }`$, where $`e`$ is the orbital eccentricity of the satellite. Thus, the perigee precession $`\dot{\omega }`$ for LAGEOS is difficult to measure because its orbital eccentricity $`e`$ is $`4\times 10^3`$. The orbit of LAGEOS II is more eccentric, with $`e`$ 0.014, and the Lense-Thirring drag of the perigee of LAGEOS II is almost twice as large in magnitude as that of LAGEOS. Over about 4 years, we obtained a root mean square of the residuals of the LAGEOS II perigee of about 25 mas , whereas the total Lense-Thirring effect on the perigee is, over 4 years, $`228`$ mas.

To precisely quantify and measure the gravitomagnetic effects we have introduced the parameter $`\mu `$ that is by definition 1 in general relativity  and zero in Newtonian theory.

The main error in this measurement is due to the uncertainties in the Earth’s even zonal harmonics and their time variations. The unmodeled orbital effects due to the harmonics of lower order are comparable to, or larger than, the Lense-Thirring effect. However, by analyzing both the JGM-3 and the EGM-96 models with their uncertainties in the even zonal harmonic coefficients and by calculating the secular effects of these uncertainties on the orbital elements of LAGEOS and LAGEOS II, we find  that the main sources of error in the determination of the Lense-Thirring effect are concentrated in the first two even zonal harmonics, $`J_2`$ and $`J_4`$. We can, however, use the three observable quantities $`\dot{\mathrm{\Omega }}_I`$, $`\dot{\mathrm{\Omega }}_{II}`$ and $`\dot{\omega }_{II}`$ to determine $`\mu `$ , thereby avoiding the two largest sources of error, those arising from the uncertainties in $`J_2`$ and $`J_4`$. We do this by solving the system of the three equations for the three residuals $`\delta \dot{\mathrm{\Omega }}_I`$, $`\delta \dot{\mathrm{\Omega }}_{II}`$ and $`\delta \dot{\omega }_{II}`$ in the three unknowns $`\mu `$, $`J_2`$ and $`J_4`$, obtaining:

$$\delta \dot{\mathrm{\Omega }}_{LAGEOSI}^{Exp}+c_1\delta \dot{\mathrm{\Omega }}_{LAGEOSII}^{Exp}+c_2\delta \dot{\omega }_{LAGEOSII}^{Exp}=$$

$$=\mu (\mathrm{\hspace{0.17em}31}+\mathrm{\hspace{0.17em}31.5}c_1\mathrm{\hspace{0.17em}57}c_2)mas/yr+othererrors\mu (60.2mas/yr),$$
(3)
where $`c_1=0.295`$ and $`c_2=0.35`$. Equation (3) for $`\mu `$ does not depend on $`J_2`$ and $`J_4`$ nor on their uncertainties; thus, the value of $`\mu `$ that we obtain is unaffected by the largest errors, which are due to $`\delta J_2`$ and $`\delta J_4`$, and is sensitive only to the smaller errors due to $`\delta J_{2n}`$, with $`2n6`$.

Similarly, regarding tidal, secular, and seasonal changes in the geopotential coefficients, the main effects on the nodes and perigee of LAGEOS and LAGEOS II, caused by tidal and other time variations in Earth’s gravitational field , are due to changes in $`J_2`$ and $`J_4`$. However, the tidal errors in $`J_2`$ and $`J_4`$ and the errors resulting from other unmodeled, medium and long period, time variations in $`J_2`$ and $`J_4`$, including their secular and seasonal variations, are eliminated by our combination of residuals of nodes and perigee. In particular, most of the errors resulting from the 18.6- and 9.3-year tides, associated with the lunar node, are eliminated in our measurement. An extensive discussion of the various error sources that can affect our result is given in , only a brief discussion of the error sources is given below.

Let us now report the main results of our measurements. In Fig. 1, we display the improved analysis  (obtained with the linear combination of the residuals of the nodes of LAGEOS and LAGEOS II and perigee of LAGEOS II according to Eq. (3)) using the static and tidal Earth gravitational model EGM-96, we also refined the non-gravitational perturbations model, the total period of observations was of 4 years, longer of about 1 year with respect to the previous observational period corresponding to ref. . We only removed three small periodic residual signals and the small observed inclination residuals. The removal of the periodic terms was achieved by a least squares fit of the residuals using a secular trend and three periodic signals with 1044-, 820-, and 569-day periods, corresponding to, respectively, the nodal period of LAGEOS, and the perigee and nodal periods of LAGEOS II. The 820-day period is the period of the main odd zonal harmonics perturbations of the LAGEOS II perigee; the 1044- and 569-day periods are the periods of the main tidal orbital perturbations, with $`l`$=2 and $`m`$=1, which were not eliminated using Eq. (3). Some combinations of these frequencies correspond to the main non-gravitational perturbations of the LAGEOS II perigee. We notice that this analysis, using EGM-96 and its accurate tidal model, is substantially independent of the removed signals, whereas the previous analysis , was in part sensitive on the periodic terms included in the fit. In other words, our value (Fig. 1) of the secular trend does not significantly change by fitting additional periodic perturbations, and indeed, even the fit of the residuals with a secular trend only, with no periodic terms, just changes the slope by about 10 $`\%`$. Nevertheless, in this case, the root mean square of the post-fit residuals increases by about four times with respect to Fig. 1.

Our best-fit straight lines of Fig. 1, through the combined residuals of nodes and perigee, has a slope $`\mu ^{Measured}1.1\pm 0.03,`$ where 0.03 is the standard deviation of the fit. This combined, measured, gravitomagnetic perturbation of the satellites’ orbits corresponds, in a 4-year period, to about 16 m at the LAGEOS altitude, that is, about 265 mas.

The root mean square of the post-fit combined residuals corresponding to Fig. 1 is about 9 mas. Our total systematic error is estimated to be of the order of 30$`\%`$-40$`\%`$ of $`\mu _{GR}`$ corresponding to the previous analyses of ref. , and of the order of 20$`\%`$-25$`\%`$ of $`\mu _{GR}`$ corresponding to Figure 1 .

Using the JGM-3 covariance matrix, we found the errors due to the uncertainties in the even zonal harmonics $`J_{2n}`$, with $`2n6`$, to be: $`\delta \mu ^{evenzonals:J_{2n}J_6}`$ $`\stackrel{<}{}17\%`$ of $`\mu _{GR}`$, and using the EGM-96 covariance matrix: $`\delta \mu ^{evenzonals:J_{2n}J_6}\stackrel{<}{}13\%`$ of $`\mu _{GR}`$. The errors in the modeling of the perigee rate of LAGEOS II due to the uncertainties in the odd zonal harmonics $`J_{2n+1}`$ are, with EGM-96: $`\delta \mu ^{oddzonals}\stackrel{<}{}2\%`$ of $`\mu _{GR}`$. Using the EGM-96 tidal model, we estimated the effect of tidal perturbations and other variations of Earth gravitational field to be $`\delta \mu ^{tides+othervariations}\stackrel{<}{}4\%`$ of $`\mu _{GR}`$. On the basis of analyses  of the non-gravitational perturbations — in particular, those on the perigee of LAGEOS II — we found $`\delta \mu ^{nongravitational}\stackrel{<}{}13\%`$-20$`\%`$ of $`\mu _{GR}`$, including uncertainties in the modeling of the satellites’ reflectivities, and the error due to uncertainties in the orbital inclinations of LAGEOS and LAGEOS II was estimated to be $`\delta \mu ^{inclination}\stackrel{<}{}5\%`$ of $`\mu _{GR}`$.

Taking into account all of these error sources, we arrived at a total root-sum-square error $`\stackrel{<}{}20\%`$-25$`\%`$ of $`\mu _{GR}`$. Therefore, over an observational period of 4 years and using EGM-96, we determined $`\mu ^{Measured}=1.1\pm 0.25`$ , where 0.25 is the estimated total uncertainty due to all the error sources — see Fig. 1.

We finally briefly report on our latest, 2002, measurement of the Lense-Thirring effect over 7.3 years of data of LAGEOS and LAGEOS II, i.e. over an observational time nearly double than the longer period of our previous analyses, obtained by only modeling the radiation pressure coefficient of LAGEOS II (see Fig. 2) .

This recent measurement fully confirms and improves our previous results: the Lense-Thirring effect exists and its experimental value, $`\mu (1\pm 0.02)\pm 0.2`$ (where $`\pm 0.02`$ is the standard deviation of the fit and $`\pm 0.2`$ is the estimated total systematic error), fully agrees with the general relativity prediction of frame-dragging. It is important to notice that: (1) in the analysis corresponding to Fig. 2 we only modeled on LAGEOS II the radiation pressure coefficient of the satellite, i.e. the reflectivity coefficient, $`C_R`$, and no other parameters such as the along-track accelerations as in our previous analyses corresponding to Fig. 1 and ref. ; (2) the RMS of the residuals corresponding to Fig. 2 is about 10 milliarcsec whereas the total measured signal is about 440 milliarcsec, and finally (3) the quality of the fit and of the corresponding measurement can be further improved by further reducing the RMS of the 15-day fits (corresponding to each point of Fig. 2) with further processing of the data using GEODYN/SOLVE, thus further reducing the RMS of the final fit of Fig. 2.

In conclusion, the Lense-Thirring effect exists and its experimental value, $`\mu 1\pm 0.2`$, fully agrees with the prediction of general relativity .


text:
# The radial distributions of a heavy-light meson on a lattice

## 1 INTRODUCTION

A knowledge of the radial distributions of quarks inside hadrons could lead to a better understanding of the QCD description of these hadrons and possibly suggest forms for phenomenological models. As a step in this direction, we have measured the charge (vector) and matter (scalar) densities of a heavy-light meson.

More explicitly, the heavy-light meson is simplified to being an infinitely heavy quark ($`Q`$) and an antiquark ($`\overline{q}`$) with a mass approximately equal to that of the strange quark. The physical meson nearest to this idealised meson is the $`B_\mathrm{s}`$(5.37 GeV). We consider two flavours of dynamical quarks ($`u`$ and $`d`$).

## 2 MEASUREMENTS

### 2.1 Energies

The basic quantity for evaluating the energies of a heavy-light meson is the 2-point correlation function – see Fig. 1 a) –

$$C_2(T)=U^Q(𝐱,t,T)P(𝐱,t+T,t),$$
(1)
where $`U^Q(𝐱,t,T)`$ is the heavy (infinite mass)-quark propagator and $`P(𝐱,t+T,t)`$ the light anti-quark propagator. The $`\mathrm{}`$ means the average over the whole lattice.

### 2.2 Heavy-light radial correlations

When measuring the radial correlations, the basic quantity is the 3-point correlation function – see Fig. 1 b) –

$$C_{3,c(m)}(T,𝐫)=U^QP_1\mathrm{\Theta }P_2,$$
(2)
where the $`P_{1,2}`$ are the light anti-quark propagators that go from the $`Q`$ at time $`t_1`$ to the point $`𝐫`$ at $`t=0`$ and then return to $`Q`$ at time $`t_2`$. Here $`T=t_1+t_2`$. The probe $`\mathrm{\Theta }`$ at $`𝐫`$ is here considered to have two forms: i) $`\mathrm{\Theta }=\gamma _4`$ for measuring the charge distribution (c) of the $`\overline{q}`$ and ii) $`\mathrm{\Theta }=1`$ for measuring the matter distribution (m).

The strategy is to first fit the $`C_2(T)`$ by the approximate expression

$$\stackrel{~}{C}_2(T)=\underset{\alpha }{}v^\alpha \mathrm{e}^{E_\alpha T}v^\alpha .$$
(3)
This results in the energy eigenvalues $`E_\alpha `$ and their eigenvectors $`v^\alpha `$ (Ref. ).

Given the $`E_\alpha `$ and $`v^\alpha `$, the $`C_{3,c(m)}(T,𝐫)`$ are then fitted by

$$\stackrel{~}{C}_{3,c(m)}(T,𝐫)=\underset{\alpha \beta }{}v^\alpha \mathrm{e}^{E_\alpha t_1}x^{\alpha \beta }(r)\mathrm{e}^{E_\beta t_2}v^\beta ,$$
(4)
where the $`x^{\alpha \beta }(r)`$ are varied giving directly the desired charge or matter density. In this work two levels of fuzzing with 2 and 8 iterations were used in addition to the unfuzzed lattice. This enabled us to extract information of the excited states as well.

## 3 RESULTS

### 3.1 Fits to the charge and matter densities

It is convenient to parametrize the radial correlations in some simple way. We have considered pure exponential and Yukawa fits to the data as well as their lattice forms

$$\frac{\pi a}{2r^\mathrm{E}L^3}\underset{𝐪}{}\frac{\mathrm{cos}(𝐫𝐪)}{\left[_{\begin{array}{c}i=1\end{array}}^3\mathrm{sin}^2(\frac{aq_i}{2})+(\frac{a}{2r^\mathrm{E}})^2\right]^2}$$
(5)
(lattice exponential) and

$$\frac{\pi }{aL^3}\underset{𝐪}{}\frac{\mathrm{cos}(𝐫𝐪)}{_{i=1}^3\mathrm{sin}^2(\frac{aq_i}{2})+(\frac{a}{2r^\mathrm{Y}})^2}$$
(6)
(lattice Yukawa).

They all fit the data well, and the lattice forms are able to reproduce some of the structure (lattice artefacts) — see Fig. 2. The exponential decay of the charge and matter distributions can be described by mesons of mass $`0.9\pm 0.1`$ and $`1.5\pm 0.1`$ GeV respectively. The matrix elements with excited states are shown in Figs. 34.

### 3.2 Sum rules

It is also of interest to consider the sum rules that sum over all values of $`r`$ i.e.

$$I_{\text{c(m)}}=\frac{_𝐫C_{3,c(m)}(T,𝐫)}{C_2(T)}.$$
(7)
For the charge and matter densities this yields $`I_\text{c}=1.4(1)`$ and $`I_\text{m}=0.9(1)`$ respectively. The corresponding values for the quenched approximation are $`I_\text{c}=1.30(5)`$ and $`I_\text{m}=0.4(1)`$. By charge conservation we should get $`I_\text{c}=1`$ in the continuum limit. Renormalisation effects of $`0.7`$ enter due to the finite lattice spacing.

## 4 CONCLUSIONS AND FUTURE

The $`S`$-wave charge and matter densities can be measured quite reliably out to $`0.8`$ fm and $`0.5`$ fm respectively. No difference can be seen when comparing dynamical and quenched fermions. This is probably due to the rather heavy mass of the light quarks.

This work is by no means the last word on the subject, and a lot remains to be done. The next step could be, for example, measuring the $`P`$-, $`D`$\- and $`F`$-wave densities. The correlations in the baryonic and $`(Q^2\overline{q}^2)`$ systems should be studied as well. We would also like to understand the densities phenomenologically using the Dirac equation. To check the continuum limit we should, of course, repeat the measurements using larger $`\beta `$ and bigger lattices.

The authors wish to thank the Center for Scientific Computing in Espoo, Finland for making available resources without which this project could not have been carried out. One of the authors (J.K.) wishes to thank the Magnus Ehrnrooth Foundation for financial support.