text:
# 1 Introduction

## 1 Introduction

Mathematical models of athletic sprinting/running performances, in some form or another, date back to at least the early 1900s. A. V. Hill  was one of the pioneers to study a sprinter’s “velocity curve”, or i.e. the runner’s speed as a function of time. Mathematician Joseph Keller  formulated a simplistic model to predict potential race times over a variety of distances. His method involved optimizing the set of coupled differential equations

$`\dot{d}(t)`$ $`=`$ $`v(t)`$
$`\dot{v}(t)`$ $`=`$ $`f(t)\alpha v(t),`$ (1)
with the race distance calculated as the integral

$$d=_0^Tv(t)𝑑t.$$
(2)
In Equation (1), $`\alpha `$ is a decay constant which places upper limitations on maximum velocities and velocity profiles. Solutions are obtained by solving Equation 2 subject to the initial conditions $`v(0)=0`$, $`d(0)=0`$. Keller’s original model for sprint races proposed a constant propulsive force, $`f(t)=F`$, arguing that the athlete must use his/her entire strength to maximize their performance. In the works of Tibshirani  and Mureika , the propulsive term is adjusted to vary with time, since such intense physical exertion will inevitably introduce muscular fatigue. Explicitly, a linear decrease was chosen, $`f(t)=Fct`$, where $`c>0`$.

## 2 Quasi-physical model

While Tibshirani’s model yields an approximate match of final times, it does not completely represent a simulation of an actual race. A realistic model of a sprint race should be able to accurately reproduce the critical 10 m split data which is available111A “split” is defined as the elapsed time at each 10 m division of the race, or the time-interval in which a 10 m stretch is covered. The athlete’s instantaneous velocities may be obtained for each 10 m increment, but due to technological and financial restrictions, this type of data is currently rare.. In order to better realize this approximation, the equation of motion (1) is modified as follows:

$$\dot{v}(t)=f_s+f_mf_vf_d.$$
(3)
The term “quasi-physical” has been adopted to highlight the fact that Equation (3) is actually a mix of both mathematical and physical components. It is not a fully physical representation, nor is is a purely mathematical model, although it can be used to effectively study and estimate key physical quantities (such as the elusive drag area of sprinters).

The components of Equation (3) are defined in the following subsections.

### 2.1 Drive term $`f_s`$

A sprint race can be broken into roughly three different phases: the drive, transition, and maintenance phases. Previous models of sprint performances (e.g.  and references therein) assign a singular propulsive term in the equation of motion. These do not explicitly account for the drive phase, in which a sprinter begins a race from a crouched start. From this position, the athlete is able to achieve greater speed due to increased application of forces (efficient drive posture). This phase of the race lasts for only about the first 25-30 metres, at which time the sprinter has transitioned to an upright running stance.

In the current model the following form of the drive term is adopted:

$$f_s=f_0\mathrm{exp}(\sigma t^2),$$
(4)
where $`f_0`$ is the magnitude of the drive, and $`\sigma `$ a constant to be determined. The $`t^2`$ dependence ensures that the term drops off rapidly as the race progresses. In fact, after roughly $`t=3`$ seconds, the magnitude of the term drops to less than $`0.1\%`$ of its original value, which would correspond roughly to the appropriate distance described above.

### 2.2 Maintenance term $`f_m`$

This term is left relatively intact from the previous models, but assumes a slightly different interpretation. Once the drive term has dissipated, this term is the only remaining propulsive component of Equation 3, and represents the maintenance phase of the race. Explicitly,

$$f_m=f_1\mathrm{exp}(ct).$$
(5)
Note that, unlike the linear maintenance term assumed in , the time-dependence assumed herein is non-negative $`t`$. Certainly, this is a much more realistic assumption.

Due to the adoption of the drive term (4), the value of $`f_1`$ tends to be lower than predicted in previous models. While this does not immediately effect predictions for 100 metre performances, it does have significant implications for 200 metre dash simulations (for which previous models have made rather generous performance predictions; these are discussed in detail in a forthcoming article ). This difference is further discussed in Section 3.3.

### 2.3 Velocity term $`f_v`$

The first of two counter-propulsive terms, the velocity term is another relic of the original Keller model, which accounts for a predicted $`v(t)`$ dependence in the equation of motion. Such a term has some conceivable physical interpretations: there must exist a physical barrier which limits the maximum speed of a human, based more than just on pure muscular strength. It seems logical to assume that a sprinter’s acceleration is curtailed with increasing speed (e.g. leg-turnover or stride rate are physically and physiologically-constrained quantities). The term is written

$$f_v=\alpha v(t),$$
(6)
with $`\alpha >0`$ a positive constant. A reasonable assumption might be that the value of this parameter (based on the interpretation above) is a “physiological constant” for most sprinters, and should not vary significantly from a prescribed value. Most world-class sprinters show stride-rates of about 4-5 strides/second (see e.g. the analysis in ).

### 2.4 Drag term $`f_d`$

Along with the drive term, this is the most significant adjustment to the model. Explicitly,

$$f_d=\left(1\frac{1}{4}\mathrm{exp}\{\sigma t^2\}\right)\rho (H)A_d(v(t)w)^2.$$
(7)
The athlete is thus modeled as a thin slab of frontal cross-sectional area $`A`$, which is a sufficient approximation for the effects considered herein. A critical factor in (7) is the “modified drag area” $`A_d=C_dA/M`$. Here, $`C_d`$ the drag coefficient, and $`M`$ the sprinter’s mass (since these models are expressed in force per unit mass). The exact values of these parameters are unknown, since they can only be measured experimentally. Previous works  have suggested that $`A`$ falls around $`0.45`$ m2. Additionally, these authors suggest that the drag coefficient $`C_d`$ assumed a value on the order of $`0.91.0`$, but a recent suggestion by Linthorne  pegs the value as closer to $`0.6`$. The data obtained from this model is consistent with this statement, and in fact indicates that $`C_d`$ may assume an even lower value (see Section 3).

The expression $`(11/4\mathrm{exp}(\sigma t^2))`$ in (7) is designed as an initial correction to the cross-sectional area $`A`$, representative of the transition from a crouched to upright position. Although the value $`1/4`$ is purely subjective, the actual magnitude of the correction does not significantly affect the results. It does, however, allow for a slightly greater acceleration over the first 10 m, and is a more realistic assumption than having a constant $`A`$.

The important of an ambient wind in a sprint race is exemplified by this term, since its magnitude is a function of $`(v(t)w)^2`$. Any non-zero tail-wind reduces the effective drag experience by the sprinter, allowing him/her greater acceleration, a greater top speed, and hence a faster overall time. A positive wind-speed corresponds to a tail-wind (i.e. wind in the direction of motion of the sprinter), while a negative wind-speed denotes a head-wind. The IAAF (International Amateur Athletic Federation) has adopted a limit of +2.0 ms-1 for a supporting wind, above which a race is deemed “wind-assisted” . Since such times are not recognized as legal performances, they cannot be ratified as World Records.

Wind-speed has the dominant effect on drag, but an additional factor which affects this term is the atmospheric density, $`\rho (H)`$, where

$$\rho (H)=\rho _0\mathrm{exp}(0.000125H).$$
(8)
Here, $`\rho _0=1.184`$ gm-3 is the sea-level density at 25 degrees Celsius, and $`H`$ is the measured elevation/altitude (in metres) (Dapena and Feltner  propose a second-order-in-$`H`$ correction to the altitude, which is not included here).

Races which are run at altitudes above $`H=1000`$ m are deemed altitude-assisted, but unlike wind-assisted marks, these can be ratified as records. For example, at the 1968 Olympics in Mexico City ($`H2250`$ m), World Records were set in both the 100 m and 200 m dashes, thanks to the considerable altitude (the density of air in Mexico City is roughly 76$`\%`$ of $`\rho _0`$). Pietro Mennea’s former 200 m record of 19.72 s was also set at altitude. Furthermore, the fastest-ever recorded 100 m clocking is 9.69 seconds by Obadele Thompson of Barbados, in April 1996 (at which time the official World Record was 9.85 s). This race was run in El Paso, Texas ($`H1300`$ m) with a tail-wind of $`+5.7`$ ms-1. Thompson’s previous 100 m performances were all slower than 10.00 s, which serves to demonstrate the extremal benefits of drag reduction.

## 3 Model parameters and simulation results

The coupled equations (3) and $`\dot{d}(t)=v(t)`$ cannot be solved analytically, so one must resort to numerical methods. Thanks to the significant increase in modern processor speeds, the values of the parameters $`(f_0,\sigma )`$; $`(f_1,c)`$; $`\alpha `$; $`A_d`$ can be rapidly isolated. This was done using a fourth-fifth order Runge-Kutta integration scheme (written in C) run on a 500 MHz Pentium III processor supporting Linux 6.0. An iterative time-step of 0.001 s was chosen for the integration.

In order to determine extremal parameters for the model, key 10 m split data for several world class 100 m races are matched. Such data has been obtained at various world-class track meets, including the 1997 and 1999 World Championships in Athletics , as well as the 1988 Olympic Games . The instantaneous velocity splits obtained in  were measured with a laser-based device known as a LAVEG, which collects its data by sensing the reflectivity of 20 ns, 904 nm pulses directed at a (linearly) receding target. The cited reliability of the distance measurements is $`\pm 1`$cm, at a sampling rate of up to 100 measurements per second .

After repeated numerical runs and adjustments to the parameters, the following set was obtained:

$`(f_0,\sigma )=(6.10,2.22)`$;

$`(f_1,c,\alpha )=(5.15,0.0385,0.3225)`$;

$`A_d=0.002875`$;

Also, $`w=0.0`$ ms-1; $`H=0`$ m. The simulation results are displayed in Table 1, and yield a raw (i.e. excluding reaction222Although they depend on the individual athlete and the conditions at the time of the race, reaction times tend to range between $`0.120.17`$ s. Reaction times below $`0.100`$ s are not allowed, as it is believed that it is not physiologically possible to surpass this limit.) time of 9.70 s, with a maximum velocity of $`11.845`$ ms-1 at 59.18 m. Such figures are in agreement with actual world class sprinters running sub-9.90 second races to better than $`\pm 0.01s`$, with the instantaneous velocity splits also matching to within 1 $`\%`$ of the actual value (see e.g. the data in ).

Of course, each individual sprinter would most likely be described by a unique set of parameters, and such a task is not within the immediate scope of this paper. Also, due to the nature of the drag term (7), the effects of cross-winds are not considered in the current form of the model.

The drive-phase correction (4) is necessary to accurately reproduce the observed velocity profile over the first 30-40 metres. For example, the velocity profiles presented in  for sprinter Donovan Bailey of Canada give 9.32 ms-1 (10 m), 10.95 ms-1 (20 m), 11.67 ms-1 (30 m), and 11.99 ms-1 (40 m). These velocity figures are doubtful, in light of the profiles presented by actual split data . It is unlikely that Bailey achieved a speed of 11.99 ms-1 as early as 40 m, and furthermore sustained speeds in excess of 12 ms-1 for an additional 40 m .

### 3.1 Wind assistance: determination of the drag coefficient and frontal cross-sectional area

The value of $`A_d`$ is isolated in conjunction with the measured effect of wind assistance by past authors , most of whom agree that a tail-wind of $`+2.0`$ ms-1 will boost a 10 s 100 m sprint by about $`0.1`$ s. The results of an earlier study  cite corrections of $`0.07`$s for a $`+2`$ ms-1 wind, although one of the authors has since produced updated results which are more commensurate with the literature . Table 1 also shows the effects of such a tail-wind on the predicted 100 m times, and accordingly predicts a boost of $`0.104s`$. Conversely, a head-wind of equal magnitude ($`w=2.0`$ ms-1) will increase the time by $`0.130`$ s to 9.830 s. Clearly, the non-linear nature of the drag term implies that equal but opposite wind speeds will not provide equal boosts. It is reasonable to assume that there could be mild variations in these corrections, depending on variations in $`A_d`$, but a full study of this effect is not the aim of this paper.

For a sprinter of mass 80 kg, the value $`A_d=0.002875`$ indicates that the effective drag area is $`AC_d=0.23`$ m2. Note that this is slightly less than the 0.3 m2 estimate of Linthorne . For a cross-sectional area between $`0.400.50`$ m2, this suggests that the drag coefficient is between 0.46 and 0.57. Note that it may actually be erroneous to assume that the drag coefficient is constant, since the could potentially be frequent transitions to laminar flow occurring throughout the race. Also, the drag area may vary slightly from athlete to athlete.

It should be noted that while these simulation times are recorded to 0.001 seconds, official race times are reported to only 0.01 s, and official wind gauge readings to 0.1 ms-1. In actual fact, electronic/photo-finish timing is recorded to 0.001 s, and the following “rounding-up” algorithm is applied to the times: unless the third decimal place is 0, the hundredth position is rounded up. For example, Canadian sprinters Donovan Bailey and Bruny Surin have both recorded 100 m times of 9.84 s, but their pre-rounded times were actually 9.835 s and 9.833 s respectively . A time of 9.831 s would be reported as 9.84 s, while 9.830 s would earn a 9.83 s rounding.

The reported value is measured from a gauge of height no more than 1.22 m, placed within 2 m from lane 1 on the in-field at a distance of 50 m from the finish line. The wind-speed is sampled for a duration of 10 s from the start of the race . Wind speeds are measured to 0.01 ms-1, but are rounded in a similar fashion. However, Linthorne  has recently noted that measuring the wind speed in such a fashion gives an error of $`\pm `$ 0.7 ms-1 on the actual value. Such a discrepancy would certainly impact the validity of the reported times.

### 3.2 Altitude effects

Although altitude effects are secondary insofar as drag modification is concerned, their influence is not negligible. Table 2 demonstrates the modifications to the times of Table 1, subject to an increasing elevation. Non-zero wind effects are not included, but will be discussed in detail later on. The model predicts an advantage of $`0.069`$ seconds from sea level to a 2000 metre elevation. For a “Mexico City” altitude (approximately 2250 m), the correction to this 100 m time (with no wind) would be $`0.075`$ s, in rough agreement with the predictions of Linthorne  and Dapena . It should be noted that the correction figures given in Reference  are considered to be overestimated, due to an inaccurate value of the drag area.

### 3.3 Graphical analysis of propulsive forces

While hypothesized from a mathematical basis, the quasi-physical model herein provides accurate matches to real data. A study of the propulsive forces reveals a striking agreement to that presented by Dapena and Feltner . Although their analysis is for 10.90 s-caliber sprinters, the figures presented herein for sub-10 s sprinters show the same overall structure. The authors suggest that a plot of mass-normalized propulsive force versus the athlete’s velocity can be modeled by two straight lines, adjoined at a “boundary velocity” roughly 3-5 ms-1 less than maximum.

Figure 5 demonstrates the identical analysis for the model parameters listed in Section 3. Unlike Dapena and Feltner’s model, the acceleration is highly non-linear before 4 ms-1, at which point it assumes a roughly linear form to about 7 ms-1. The “boundary velocity” in this case is representative of the point where the drive term $`f_d`$ has dropped to roughly 10 % of its original value (about 1 second into the race). At this point, the curve assumes another linear form of differing slope, before approaching the maximum velocity (and again diverging from linearity).

Such an analysis can in fact be used to rule out the Keller model and its variants. Although Keller’s value of $`F`$ was in the $`12.1`$ ms-2 range , the overall form of the function does not provide an adequate global match.

Similarly, using Tibshirani’s modification , the statistical fit of race data from Donovan Bailey’s 1996 Olympic 100 m race  predicted an initial $`f(0)=F`$ value of 7.96 (Tibshirani’s own analysis predicted a lower value of 6.41, almost half that of Keller). These models are also presented in Figure 5, and are representative of roughly a 9.8-10 second race. On their own these two approximations do not support the predictions of the quasi-physical model, however it is interesting to note that Keller’s model before, and the Tibshirani-Mureika model after their intersection are a rough first-order approximation to the data. In fact, this is similar to linear model proposed in .

### 3.4 Assistance for women sprinters

Women sprinters are certainly lighter than an average 80 kg, yet their cross-sectional area is correspondingly smaller. Similarly, their stride rates tend to be lower. Thus, a first-pass approximation for women can be made by a mild adjustment to the parameters of Equation (3). Following the interpretation of Section 2.3, a lower stride rate can be taken to correspond to higher $`\alpha `$. Similarly, adjustments to the propulsive terms are reasonably made.

By setting $`f_0=5.60;f_1=4.57;c=0.037;\alpha =0.32`$, and $`A_d=0.00275`$m2 (which corresponds to a drag area of 0.18 m2 for a mass of 65 kg), one can reproduce splits which roughly replicate those observed by world class female sprinter Marion Jones (see ), and match the predicted wind and altitude advantages of  and . For the cited drag area, a drag coefficient between 0.5-0.6 would require a cross-sectional area of 0.30-0.36 m2, roughly 80% that of men (a reasonable approximation).

Using these values gives a raw time of 10.653 s. At sea-level ($`H=0`$ m), an advantage of $`0.117`$ s is predicted for a tail-wind of $`+2.0`$ ms-1 (10.536 s), and $`0.086`$ s for $`H=2500`$m. These figures are again in close agreement with those of  and Linthorne , who predicts a correction of $`0.12\pm 0.02`$ s. A wind of $`2.0`$ ms-1 at $`H=0`$ m will increase the time to 10.802 s, a difference of $`0.149`$ s.

Figures 1,2,3, and 4 give performance corrections for men and women, for altitudes ranging between $`H=02500`$ m. The sign of the correction indicates the impact of the associated wind and altitude on the base time $`t_0`$ at 0 m altitude and 0-wind, i.e. $`\mathrm{\Delta }t=t_wt_0`$. A negative correction indicates that the time $`t_w`$ is faster than the base (e.g. for tail winds), and vice versa (note the change in sign relative to the Tables). The predicted corrections for women cannot be applied to slower men’s races, since the assumed drag area is much lower. If the modified drag area is kept as 0.002875 m2kg-1, the resulting boost is $`0.122`$ s for a $`+2.0`$ms-1 wind and a lapse of $`0.154`$ s for a $`2.0`$ms-1 wind.

## 4 Sensitivity to parameter variation

### 4.1 Correction for late-race velocity drop-off

The measured velocity data presented in Reference  indicate that the model predictions drop off faster than the measured data for the latter half of the race. It should be noted that this discrepancy can be accounted for by lowering the power of $`t`$ in maintenance term (5), and hence the overall shape of the associated “energy envelope”. In fact, the author notes that a substitution $`tt^\beta `$, with $`\beta 0.950.98`$ (and mild adjustments to $`f_1`$) provides a closer match to the data.

While such modifications are of interest to providing a fully realistic simulation, variation of this parameter is an over-complication of an already-intricate model, and is not essential for the study conducted herein. A full study of such time dependence is postponed for future work, and in particular its implications for the 200 m dash are discussed in .

Furthermore, it was noted that variation of the $`(f_1,\alpha ,c)`$ parameters did not significantly affect the velocity profile in the early part of the race, but could be used to adjust the late-race profile. Individual adjustments of these parameters to produce a 0.1 s advantage in the final time presented minimal effects on the profile for $`d<50`$ m. Conceivably, individual athletes could be represented by their own set of parameters, which could help account for variability in observed race times. Despite these minor velocity adjustments over 60-100 m, a $`+2.0`$ ms-1 assisting wind again provided boosts between $`0.1010.103`$ s.

### 4.2 Slow starts and the time dependence of $`f_s`$

It is possible to use the model to simulate “slow starts” by suitable readjustment of the drive term parameters. In fact, by modifying the time-dependence of Equation 4, one may reproduce starts in which athletes did not accelerate as quickly as others, but achieved higher overall velocities at each split in the earlier portions of the race.

Figure 6 demonstrates such a simulation, with adjusted time-dependence $`f_d=f_1\mathrm{exp}(\sigma t)`$, and $`\sigma =1.68`$, as compared to the simulation of Section 3. The simulated splits are shown in Table 4; such data is consistent in profile with those of Canadian sprinter Donovan Bailey from the 1997 World Championships in Athens . One can note this simulation lags behind the other until about 70 metres, although achieving higher maximal velocities. In this case, the finish is within 0.005 seconds. This type of behavior is frequently observable at any competition, and is characteristic of many famous world class sprinters (including the aforementioned Bailey, as well as American track legend Carl Lewis).

## 5 Variable wind speeds

As mentioned previously, the wind speeds are sampled for a period of 10 seconds following the start of the race, after which the average value of the wind-speed is taken as the official gauge reading. The examples cited previously have assumed that the wind velocity is constant. However, this is probably more of ideal situation than not. A more realistic scenario is one in which the wind speed is variable, but averages out to a value which may be unrepresentative of the true conditions.

Following the example of Dapena and Feltner , the results of a variable-wind race are given in Table 3. Four separate variability conditions are simulated, with the constraint that the gauge reading (i.e. the mean wind velocity over 10 seconds) be the legal limit of $`+2.0`$ ms-1. These conditions include:

1. Step function: $`w(t)=4\mathrm{\Theta }(t5)`$
2. Step function: $`w(t)=4[1\mathrm{\Theta }(t5)]`$
3. Linear: $`w(t)=2/5t`$
4. Linear: $`w(t)=2/5(10t)`$

Here, $`\mathrm{\Theta }(t)`$ is the Heavyside function, with $`\mathrm{\Theta }(tt_0)=0t<t_0`$, and $`\mathrm{\Theta }(tt_0)=1`$ otherwise. The time-averaged wind-speed is calculated in the usual fashion, $`\overline{w}=T^1_0^Tw(t)𝑑t`$, with $`T=10`$seconds.

While the constant tail-wind speed of $`+2.0`$ ms-1 predicts a boost of 0.104 s for the present model, the above conditions predict a rather large range of variation. Cases 1 and 3 show the smallest advantages, due to minimal (or zero) wind conditions in the drive phase (hence lower overall accelerations in the first half of the race). These cases also show the lowest peak velocities in this boost scenario, although still lower than the constant wind case. Note that the peak velocities occur much later than in the base case and the constant-wind case.

Case 4 shows a peak velocity at roughly the same location as in the base 9.700 s run (a difference of only 19 cm), although with a much higher magnitude ($`+0.182`$ ms-1). Case 2, on the other hand, shows almost the same velocity, but at a much earlier mark (54.95 m). This decidedly premature maximum is no doubt due to the overpowering contributions of the velocity term $`f_v`$ once the assisting tail-wind has subsided. The sprinter is physically unable to achieve a higher velocity.

Dapena and Feltner only consider Cases 1 and 2 in their analysis (for the time-averaged +2 ms-1 wind). However, it is interesting that the predictions of this model in fact are opposite of their assertions. For Case 1, they predict a boost of 0.066 s, while their Case 2 shows a boost of 0.060 s, both lower than their constant-wind case (0.070 s). The conditions of Case 1 and Case 3 allow the maximum velocity to be achieved later in the race, which may not have an issue addressed by the authors in question. Their underestimation of the initial boost may also be a contributing factor, and adjustment of their parameters to suit the current estimates of Dapena  may yield differing results. Whether this is a false estimate in the figures of Dapena and Feltner or a short-coming of the current model is an open question. Observation of this velocity-limitation could yield credibility to form of $`f_v`$.

Certainly, such wind speed variations considered in this section can further contribute to small velocity discrepancies between the model with constant wind speed, as well as variations in the actual race splits.

### 5.1 Variable 0-wind average

The analysis above begs the question: what kind of performance variations could be expected for a time-averaged wind speed of 0 ms-1? Table 5 demonstrates such potential discrepancies, with linear and sinusoidally-varying wind speed functions (the unrealistic step functions are omitted here). The linear winds simply range from equal-but-opposite maximum/minimum values (with 0 ms-1 occurring at the 5 s midpoint), the cosine functions begin and end at the indicated maximum and minimum, and the sine functions achieve a the maximum indicated with 0 ms-1-endpoints. Absolute maxima of 2 and 4 ms-1 are chosen.

The simulations for the 2 ms-1 absolute wind show variations between 0.034 s slower and +0.023 s faster than the actual 0 ms-1 wind-speed values of Table 1. Interestingly enough, the greatest advantage comes from a cosine wind which begins and ends as a head-wind, the form of which allows for greater acceleration through middle of the race. Conversely, the greatest disadvantage is from the equivalent negative-cosine wind. As the overall extremum of the wind increases, the differences become more acute (although note a much higher disadvantage when the wind is negative from about 30 m onward). Figure 7 presents the velocity profiles of several of these variable conditions, compared to the base case of Section 3.

## 6 Comparison to previous wind-correction estimates

Recently, a “back-of-the-envelope” expression for potential wind and altitude corrections was obtained in :

$$t_{0,0}\left[1.0270.027\mathrm{exp}(0.00125H)(1wt_{w,H}/100)^2\right]t_w,$$
(9)
Here, $`t_{w,H}`$ is the official race time (run with wind $`w`$ at altitude $`H`$), and $`t_{0,0}`$ the time for $`w=0`$ ms-1 at sea level. This is derived in part from Equation 7 by selecting a constant propulsive force, and an average velocity $`v=100/t_{w,H}`$. The assumption is made that a sprinter expends roughly $`\delta =2.7\%`$ of his/her energy fighting drag at sea level. For non-zero $`w,H`$, $`\delta \rho (H)A_d(v_{w,H}w)^2`$. The effects of wind are thus calculated by assuming the adjusted drag impacts the sprinter’s velocity by the ratio of forces $`F(v_{w,H}w)/F(v_0)`$, where $`v_{0,0}100/t_{0,0}`$ is the average velocity at sea level with no wind for a race-time $`t_{0,0}`$, with $`v_{w,H}=100/t_{w,H}`$ the wind/altitude-influenced velocity/time.

This ready-to-use expression provides a good match to the predictions of this model, as well as those of Dapena  and Linthorne . Table 6 demonstrates the corrective potential of Equation (9) subject to the model parameters of Section 3. Note that the approximation becomes slightly worse for increasing absolute wind speeds, however the magnitude of these variations is less than 0.5% for $`w[5,+5]`$ ms-1. Such errors can certainly be accounted for by mis-estimations of the input variables.

The mismatch may also result from the constant velocity assumption in the derivation. The interested reader is referred to  for the complete derivation.

## 7 What would have been the 100 m world record in 1988?

On September 24th at the 1988 Seoul Olympics, Ben Johnson of Canada clocked an astounding 100 metre World Record time of 9.79 s, bettering his previous record of 9.83 s. Although this mark (and most of Johnson’s other records) were stricken from the books due to the infamous steroid scandal which followed, the Seoul mark did not reflect the true potential of this remarkable athlete. At about 15 m from the finish line, Johnson looked over his shoulder to gauge his lead over American Carl Lewis, and then raised his arm in victory as he coasted through the remained of the race. The question lingers: what would have been his World Record time, had he not “stopped” at about 85 m?

The model discussed herein can be used to obtain projections of his potential performances. Since no instantaneous velocity data is available for this race, the parameters in Equation 3 are selected to match the 10 m splits. Table 7 shows these values (obtained from ), along with two sample simulations. The measured wind-speed of the race was $`w=+1.1`$ms-1, and the modified drag area is kept as in Section 3. The elevation of Seoul is roughly that of sea level, so altitude corrections will be minimal in this case.

Effectively, Johnson could have lowered his previous record by almost 0.1s, a quantum leap in the event. Depending on the chosen parameters, the simulations predict times in the range of 9.60-9.62s, which round to about 9.73-9.75s once reaction time is included. In this case, two sets are selected to match the recorded splits, with the parameter values the same as in Section 3, except variation in the values of $`(f_0,f_1)`$. The modified values are $`(f_0=6.20;f_1=5.16)`$ (Simulation 1), and $`(f_0=6.18;f_1=5.17)`$ (Simulation 2). Note that the maximum velocity of the simulations do not exceed 12 ms-1. While Johnson may have slightly topped this value in the real race, he most certainly did not achieve 13.1 ms-1, a feat periodically attributed to him .

In 1999, American Maurice Greene reset the World Record to 9.79 s during a race in Athens, Greece. Up to slight variations in altitude, the city is also at sea level, where the wind gauge read a calm $`+0.1`$ ms-1. Adjusting this variable in the Johnson simulations, the time is scaled down to about 9.665 s (9.797 s after reaction) for Simulation 1, and 9.653 s for Simulation 2. (9.785 s after reaction). Thus, Johnson’s former World Record would have effectively been equal to the 9.79 s set in 1999 by Greene.

## 8 General conclusions and future considerations

The quasi-physical model considered herein provides a good match to measured split data, particularly in the drive phase of the race. The wind and altitude corrections for both men and women are consistent with those proposed by previous authors, although this models yields stronger deviations from the norm for variable wind conditions. If this type of scenario occurs in a real competition, then the times may appear much faster (or slower) than expect from the measured (time-averaged) wind.

Upon comparison with the previous theoretical estimates and complementary field studies, it is clear that this model is both a realistic representation of short-sprint races, as well as a useful tool which can assist in the physiological and biomechanical study of the sport.

A forthcoming manuscript  will discuss the effects of wind and altitude assistance in the 200 metre dash.


format: text

text:
# EFFECTS OF THE AUGUST 11, 1999 TOTAL SOLAR ECLIPSE AS DEDUCED FROM TOTAL ELECTRON CONTENT MEASUREMENTS AT THE GPS NETWORK

## 1 Introduction

Experimental observations of the ionosphere at the time of solar eclipses provide the source of information about the character of behavior of the various ionospheric parameters. Regular ionospheric effects of solar eclipses are fairly well understood. They imply an increase of effective reflection heights, a reduction in concentration in the F-layer maximum, and a decrease in total electron content (TEC) in the ionosphere, which is typical of the transition to the nightside ionosphere (Cohen, 1984). The behavior of the above parameters can be modeled using appropriate ionospheric models (Stubbe, 1970; Boitman et al., 2000).

The main parameters of the ionospheric response are the value of the delay $`\tau `$ with respect to the eclipse totality phase, as well as its amplitude $`A`$ and duration $`\mathrm{\Delta }T`$. Almost all publications devoted to the study of the ionospheric response to solar eclipses present estimates of these parameters. A knowledge of these values makes it possible to refine, in terms of the respective aeronomic ionospheric models, the time constants of ionization and recombination processes at different heights in the ionopshere.

The measurements of $`\tau `$ were made by analyzing the characteristics of the ionosphere-reflected radio signal at vertical-incidence soundings at a network of ionospheric stations (Marriott et al., 1969; Goncharov et al., 1982; Boitman et al., 2000) by measuring the frequency Doppler shift at vertical- and oblique-incidence soundings (Boitman et al., 2000). In the cited references, the value of $`\tau `$ was found to vary from 0 (Marriott et al., 1969) to 20 min, with the amplitude $`A`$ of a decrease in local electron density of order $`9\times 10^4`$ cm-3, and the response duration of about 1 hour (Boitman et al., 2000).

Similar measurements at an ionospheric station were made during the solar eclipse of September 23, 1987 in the south-eastern Asia (Cheng, 1992). The amplitude, delay and duration of the ionospheric response were found to be $`0.3\times 10^6`$ cm-3, 18 min and 1 hour 30 min, respectively. At the same time, measurements during the eclipse of October 24, 1995 in the same region provided an estimate of $`\tau `$ of about 80 min (Huang, 1999), $`\mathrm{\Delta }T`$ = 1 hour 20 min, and $`A`$ in excess of $`1\times 10^6`$ cm-3. The measurements of $`f_0F2`$ over the Scaramanga station during a partial solar eclipse of May 20, 1996 (Anastassiadis, 1970) give the value of $`\tau `$ of about 38 min. In paper Zherebtsov et al. (1998a), the delay time of the ionospheric response to the eclipse was between 25 and 30 min (at 300 km altitude).

Interesting results were obtained in observations of an annular solar eclipse of May 30, 1984 at the Millstone Hill incoherent scatter radar (Salah, 1986). A decrease in electron density with respect to the eclipse totality phase occurred in this case 20–30 min later, with $`A`$ = $`4\times 10^6`$cm-3, and $`\mathrm{\Delta }T`$ = 2 hours 15 min. However, facilities of this kind are too few to be extensively used in experiments during solar eclipses.

The ionospheric response to the March 9, 1997 total solar eclipse was investigated using a network of ionosondes (Zherebtsov et al., 1998b). For this event, the position of the minimum of the trough in the time dependence of electron density was shifted with respect to the time of maximum occultation of the solar disk by 10–15 min, with $`A`$ = $`1\times 10^5`$ cm-3, $`\mathrm{\Delta }T`$ = 1 hour 15 min.

A large amount of data was obtained by measuring the Faraday rotation of the plane of polarization of VHF signals of geostationary satellites (Klobuchar et al., 1970; Hunter et al., 1974; Davies,1980; Rama Rao, 1982; Essex, 1982; Deshpande, 1982; Singh, 1989). These measurements revealed the eclipse-induced effect of a deep depression (decrease) in TEC with the amplitude $`A`$ varying from 2 to 14 TECU, and with a typical time of TEC depression and recovery of about several hours (the unit of measurement of TEC that is adopted in the literature corresponds to $`10^{16}`$ -2). The spread in the values of $`\tau `$ was also found to vary over a wide range, from 5 to 40 min.

According to Klobuchar et al. (1970), for a total solar eclipse that occurred in spring (March 7, 1970), the amplitude $`A`$ was 11 TECU, $`\mathrm{\Delta }T`$ = 2 hours 44 min, and $`\tau `$ of order 33 min. For the total solar eclipse of February 16, 1980, the observations at Waltair (Singh, 1989) gave the following estimates: $`A`$ = 0.11 TECU, $`\mathrm{\Delta }T`$ = 1 hour, and $`\tau `$ = 20 min. Essex et al. (1982) inferred that for the total solar eclipse of February 20, 1979, the values of $`A`$ and $`\tau `$ are 14 TECU and 33 min, respectively.

Hence a large body of experimental data do not permit us to make any reliable estimates of the basic parameters of the ionospheric response. One of the reason for such a great difference is attributable to the use of different methods of measurements which differ greatly by spatial and temporal resolution. However, the main reason is caused by dissimilar characteristics of the eclipse itself, by geophysical conditions of individual measurements, and by a large difference of the latitude, longitude and local time when experiments are conducted.

To obtain more reliable information about the behaviour of the ionosphere during an eclipse, it is necessary to carry out simultaneous measurements over a large area covering regions with a different local time. Furthermore, high spatial (of some tens of kilometres at least) and temporal (at least 1 min) resolution is needed. However, none of the above familiar methods meets such requirements.

The development of the global navigation system GPS and the creation, on its basis, of extensive networks of GPS stations (which at the end of 1999 consisted of no less than 600 sites), the data from which are placed on the INTERNET (Klobuchar, 1997), open up a new era in remote sensing of the ionosphere. At almost any point of the globe and at any time at two coherently-coupled frequencies $`f_1=1575.42`$ MHz and $`f_2=1227.60`$ MHz, two-frequency multichannel receivers of the GPS system are used to carry out high-precision measurements of the group and phase delay along the line o sight between the ground-based receiver and satellite-borne transmitters in the zone of reception. The sensitivity afforded by phase measurements in the GPS system permits irregularities to be detected with an amplitude of up to $`10^3`$$`10^4`$ of the diurnal variation of TEC.

Afraimovich et al. (1998) were the first to use in a detailed analysis of regular ionopsheric effects measurements of TEC and its gradients for the total solar eclipse of March 9, 1997, using for this purpose the GPS interferometer at Irkutsk. Their results bear witness to significant changes occurring in the process of ion production in the ionosphere during the solar eclipse, simultaneously in a large volume of space with a radius of at least 300 km at 300 km altitude. The delay of a minimum TEC value with respect to the phase of totality is about 10 min, the depth of TEC depression varies from 1 to 3 TECU, and $`\mathrm{\Delta }T`$ = 1 hour 15 min.

Boitman et al. (2000) obtained a good agreement of these data, as well of measurements for several HF Doppler sounding paths with results of a numerical simulations using the ionosphere-plasmasphere coupling model (Krinberg, 1984). Unfortunately, during this eclipse the total shadow’s path in 1997 was monitored by only few GPS stations; therefore, it was not possible to determine the spatial-temporal characteristics of the ionopsheric response.

According to the data from five GPS stations reported by Tsai and Liu (1999) for the solar eclipses of October 24, 1995 and March 9, 1997, the delay $`\tau `$ varied over a broad range, from 0 to 30 min, the amplitude $`A`$ was close to 7 TECU, and $`\mathrm{\Delta }T`$ varied from 40 min to 1 hour.

A unique opportunity to exploit the potential of the GPS network was provided by the total solar eclipse of August 11, 1999. For this period of time, the INTERNET made available the data from at least 100 GPS stations located in Western and Central Europe within and near the totality path. A massive flow of publications on the various ionospheric effects of this eclipse should be forthcoming, based on the GPS data in conjunction with the data from other observing facilities. The objective of this paper is to determine the basic parameters of the ionospheric response to the August 11, 1999 total solar eclipse using these data.

General information about the eclipse, and a description of the experimental geometry are given in Section 2. The ionospheric response to the eclipse is discussed in Section 3 using the data from reference ionospheric station Chilton. The processing technique for the GPS network data, and results derived from analysing the ionospheric response of the August 11, 1999 solar eclipse in Europe are outlined in Section 4. Results obtained in this study are discussed in Section 5.

## 2 The geometry and general information of total solar eclipse AUGUST 11, 1999

The last solar eclipse in the 20–th century began in North Atlantic, and path of the Moon’s shadow made first landfall in south-western England at 10:10 UT; the Sun at that time was at an angle of $`45^{}`$ over the eastern horizon. The centre line duration of total eclipse averaged 2 min, and the total eclipse was confined to a narrow corridor 103 km wide.

Fig. 1 shows a schematic map of the path of the Moon’s shadow crossing parts of Western and Central Europe (the data from Espenak, 1999 were used in constructing this map). The centre line of eclipse at ground level is shown as a thick line, and thin lines correspond to its southern and northern boundaries. The location of reference ionospheric station Chilton (RAL) is marked by the symbol $``$. Dots $``$ and symbols $`+`$ show the locations of the GPS stations used in the analysis; their geographic coordinates are presented in Table 1. Dots $``$, together with symbols $`+`$ are used to represent a total set of GPS stations. Symbols $`+`$ form a network of GPS stations located near the eclipse path; therefore, we designated this group as the near zone. Numbers for the longitudes $`10^{}W`$, $`0^{}`$, $`10^{}E`$, $`20^{}E`$, $`30^{}E`$, $`40^{}E`$ correspond to the local time for these longitudes.

In this paper we confine ourselves to analysing only to a region of Western and Central Europe from the coast of southern England to the point with coordinates $`56.03^{}N`$, $`37.2^{}E`$, where the totality phase was observed at 11:20 UT (14:20 LT). Thus the solar eclipse effect occurred for the conditions of the daytime summer ionosphere.

The distance along the great-circle arc between the above-mentioned extreme points is about 2900 km, with the time difference of 67 min only. Hence a distinguishing characteristic of this eclipse was the supersonic speed of the Moon’s shadow sweeping through the terrestrial surface, exceeding 720 m/s.

At ionospheric heights the totality path was travelling somewhat further south. The onset time of the totality phase for the height h=300 km over Budapest is 1.3 min ahead of that at ground level. The difference in the values of the totality phases and their onset times is caused by the Sun’s altitude over the horizon. At the time of the totality phase in Budapest (11:05 UT, or 14:05 LT), it was as small as $`59^{}`$.

The period under consideration was characterized by a low level of geomagnetic disturbance ($`Dst`$-variation from -10 to -20 nT), which simplified greatly the problem of detecting the ionospheric response of eclipse.

## 3 The ionospheric response by eclipse from date of ionospheric station Chilton

First we consider the variations of critical frequencies $`f_0F2`$ over the time interval 00:00-24:00 UT on August 11, 1999 according to the data from station Chilton (RAL) - Fig. 2a. The onset time of the totality phase of eclipse in the area of station Chilton (RAL) at 300 km altitude is shown by a vertical solid line. As might be expected under the conditions of the summer ionosphere, the mean level of $`f_0F2`$ differs only slightly for the night-time and daytime. Nevertheless, a decrease of $`f_0F2`$ during the totality phase of eclipse is sufficiently clearly distinguished.

Consider the variations of ionospheric parameters for the time interval 06:00–15:00 UT on August 11, 1999, and on the background days of August 10 and 12, using the data from station Chilton (RAL) - Fig. 3. Dots correspond to variations of critical frequencies $`f_0F2`$ (panel a), $`f_0F1`$ (panel c), apparent heights $`h^{}F`$ (panel b), and $`h^{}F2`$ (panel d) for August 11, 1999. Solid curves plot the same values that are smoothed with a time window of 60 min. For August 10 and 12, only smoothed curves are given, with the same time window of 60 min. The onset time of the totality phase of eclipse (10:16 UT) at 300 km altitude over the station is shown by a thin vertical line.

The eclipse effect is most conspicuous in the variations of critical frequencies $`f_0F2`$, whose maximum difference from the background values on August 10 and 12 at the time of reaching a minimum (10:20 UT) was up to 2 MHz. On the other hand, the amplitude of a decrease of $`f_0F2`$ for the August 11 event (after the totality phase of eclipse) does not exceed 1 MHz. The eclipse effect on other parameters is not as clearly distinguished, and becomes evident only when comparing the smoothed curves. Similar results of measurements at an ionospheric station were obtained during the total solar eclipse of September 23, 1987 in South-Eastern Asia (K. Cheng, 1992).

The $`f_0F2`$ – variations, measured at station Chilton (RAL) at time intervals of 4 min, are presented in greater detail in Fig. 4a (heavy dots). The solid curve connecting these dots is an approximating one for these values. This panel plots also the geometrical function of eclipse $`S(t)`$ at 300 km altitude, calculated for station Chilton. Minimum values of $`f_0F2`$ and $`S(t)`$ correspond to the points A and B in this figure. The respective delay $`\tau `$ between the time of a minimum of $`f_0F2`$ and of the function $`S(t)`$ is close to 4 min in this case.

The data from the station Chilton (RAL) and from GPS station HERS nearest to it are compared in the next Section.

## 4 The process of GPS–network data and results of analysis of ionospheric effect by total solar eclipse of August 11, 1999

We now give a brief account of the sequence of procedures used in processing the GPS data. Input data are represented by series of ”oblique” values of TEC $`I(t)`$, as well as by the corresponding series of elevations $`\theta (t)`$ measured from the ground, and azimuths $`\alpha (t)`$ of the line of sight to the satellite measured clockwise from the northward direction. These parameters are calculated by our developed CONVTEC program by converting the GPS-standard RINEX-files from the INTERNET. Series of elevations $`\theta (t)`$ and azimuths $`\alpha (t)`$ of the line of sight to the satellite are used to determine the location of subionospheric points. In the case under consideration, all results were obtained for larger than $`45^{}`$ elevations $`\theta (t)`$.

Fig. 4c presents the experimental geometry in the area of ionospheric station Chilton (RAL) - $``$, and GPS station HERS ($``$). Heavy dots $``$ show the centre line of eclipse at ground level, and smaller dots correspond to its southern and northern boundaries. The symbol $``$ shows the position of the subionospheric point at the time of a maximum response of TEC (see below).

Various methods for reconstructing the absolute value of TEC are currently under development using measurements of both phase and group delays; however, effective algorithms for an accurate solution of this problem are still unavailable for the different types of two-frequency receivers and operating modes of the GPS system. In this connection, for purposes of this paper, we limit our attention to considering only those TEC variations which were obtained from phase delay measurements by formula (Afraimovich et al., 1998):

$$I_p=\frac{1}{40.308}\frac{f_1^2f_2^2}{f_1^2f_2^2}[(L_1\lambda _1L_2\lambda _2)+const+nL]$$
(1)
where $`L_1\lambda _1`$ and $`L_2\lambda _2`$ are additional paths of the radio signal caused by the phase delay in the ionosphere, (m); $`L_1`$ and $`L_2`$ represent the number of phase rotations at the frequencies $`f_1`$ and $`f_2`$; $`\lambda _1`$ and $`\lambda _2`$ stand for the corresponding wavelengths, (m); $`const`$ is the unknown initial phase ambiguity, (m); and $`nL`$ are errors in determining the phase path, (m).

For an approximate specification of the TEC constant component $`I_0`$ at time intervals of 2 hours, we made use of the INTERNET data on corresponding global maps of the absolute vertical value of TEC in the IONEX format (Mannucci, 1998) - see also Fig. 2b. To normalize the response amplitude we converted the ”oblique” TEC to an equivalent ”vertical” value:

$$I(t)=I_p(t)sin(\theta (t))$$
(2)
Although, because of a strong horizontal TEC gradient and without regard for the spericity of the problem, this procedure gives a very rough result, but this result is quite acceptable because all our results were obtained for larger than $`45^{}`$ elevations $`\theta (t)`$.

With the purpose of eliminating variations of the regular ionosphere, as well as trends introduced by the satellite’s motion, we employ the procedure of eliminating the trend by preliminarily smoothing the initial series with the time window in the range from 40 to 100 min which is fitted for each TEC sampling. Such a procedure is also required for a clearer identification of the ionospheric response of eclipse which is characterized by a relatively small amplitude (see below) under the presence of space-time TEC variations that are not associated with the eclipse.

For the purposes of illustration of the data processing procedure Fig. 5 presents the filtered TEC variations $`dI(t)`$ for August 11 (thick line), and for the background days of August 10 and 12, 1999 (thin lines) for station HERS for satellite N14 (PRN 14) - c), as well as for stations ZIMM (PRN 1) - b) and WTZT (PRN 14) - a), which are separated from station HERS, respectively, by $`7^{}`$ and $`13^{}`$ in longitude eastward. Figures correspond to GPS satellite numbers. The onset time of the totality phase of eclipse at 300 km altitude for the corresponding subionospheric points is shown by vertical solid lines. As is evident from the figure, responses to the eclipse at these stations are very similar in both form and amplitude, but the response delay increases with the longitude. $`dI(t)`$ – variations for the background days of August 10 and 12, 1999 for different stations differ substantially not only from the TEC response to eclipse but also from one another.

Fig. 4b presents the filtered variations of TEC $`dI(t)`$ for station HERS for PRN 14 during August 11, 1999 (thick line). This panel plots also the geometrical function of eclipse at 300 km altitude, $`S(t)`$, that is calculated for the subionospheric point of PRN 14. A minimum value of $`S(t)`$ corresponds to the point A in the figure.

As is apparent from this figure, the form of the filtered variations is similar to a triangle whose vertex (point B) corresponds to the time at which the TEC attains its minimum value. The value of $`dI_{min}`$ itself can serve as an estimate of the amplitude of TEC response to eclipse, and the time interval between the times of intersection of the line $`dI`$=0 (points C and D) can serve as an estimate of the duration of the response $`\mathrm{\Delta }T`$. The corresponding delay $`\tau `$ between the times of the minimum of $`dI_{min}`$ and of the function $`S(t)`$ in this case was found to be 4 min, which coincides with the estimate of $`\tau `$ for ionospheric station Chilton (see Section 3). The response amplitude in this case was close to 0.18 TECU, and the response duration was 46 min.

Such $`dI(t)`$ – variations are characteristic for all GPS stations and satellite numbers 01 and 14 listed in Table 1. The choice of the same satellites, PRN 01 and PRN 14, for the entire selected set of GPS stations was dictated by the fact for these satellites a maximum value of the elevation $`\theta `$ of the line of sight to the satellite exceeded $`45^{}`$ for the time interval 10:00-12:30 UT, which reduced to a minimum the possible error of conversion to the ”vertical” value of TEC as a consequence of the sphericity.

The first line of Table 2 presents the results of a statistical processing for the entire set of GPS stations listed in Table 1. The second line includes only those stations whichlie in the immediate vicinity of the eclipse path, within $`\pm `$ $`5^{}`$ with respect to the centre line (near zone). In Table 2, the values before the bar are mean values, and those after the bar correspond to the standard deviation. The mean value of $`\tau `$ for the entire set of stations is 16 min, while for the stations in the near zone this value is 7 min. The mean value of the amplitude A=0.3 TECU for all stations and A=0.1 TECU for the near zone. The width of the TEC trough for the far and near zones $`\mathrm{\Delta }T`$ = 60 min, on average.

Fig. 6 presents a longitudinal dependence of the time position $`t_{min}`$ of minima of the curves $`dI(t)`$ for subionospheric points lying in the immediate vicinity of the eclipse band (satellites 1 and 14) - heavy dots. Dark symbols $`\mathrm{}`$ designate the times of totality phases of the geometrical function of eclipse versus longitude which are calculated for the subionospheric points. Shaded symbols $``$ correspond to delay times between maximum of the geometrical function of eclipse and a minimum in TEC.

It was found that the delay $`\tau `$ increases gradually from 4 min at the Greenwich longitude (10:23 UT, LT) to 16 min at the longitude of $`16^{}`$ (12:09 LT).

## 5 Conclusions

Our results are in good agreement with earlier measurements and theoretical estimates (see a review of the data in the Introduction). The key feature of our data is a higher reliability of determining the main parameters of the response to eclipse which is due to high space-time resolution and to the increased sensitivity of detection of ionospheric disturbances inherent in the GPS-array method which we are using.

Also, due regard must be had to the fact that the distinctive property of the eclipse under consideration was a relatively small response amplitude, which required special filtering of the TEC series (see preceding Section). The reason is that, unlike a number of eclipses for which more-or-less reliable data were obtained, this eclipse occurred in the summer season characterized by only moderate differences of the daytime and night-time ionization. Furthermore, in this situation the time variation of $`f_0F2`$ and vertical TEC near noon usually shows a minimum which, in essence, masks the eclipse effect (see also Fig. 2b). It is also vital to note that the time constant of a decrease in ionization in the $`F_2`$ maximum exceeds substantially the duration of the totality phase of eclipse, which leads also to a decrease in response amplitude.

The local time-dependence of $`\tau `$ that is revealed in this paper is in agreement with theoretical estimates reported in (Stubbe, 1970). The value of $`\tau `$ for foF2, approaching 6 min, corresponded to 13:40 LT. Using modeling methods (Ivelskaya et al., 1977) showed that the variations of the delay time $`\tau `$ of minimum local electron density $`N_e(t)`$ with respect to a minimum of the ion production function are as follows: $`\tau `$ = 1-2 min at 150 km altitude, $`\tau `$ = 3 min at 200 km, $`\tau `$ = 20 min at 300 km, and $`\tau `$ = 45 min above 600 km. In this paper $`\tau `$ is estimated, respectively, at about 3 min for 200 km altitude and 20 min for 300 km altitude for 12 LT.

## 6 Acknowledgements

We are grateful to K. S. Palamartchouk, A. V. Tashchilin and A. D. Kalikhman for their interest in this study, helpful advice and active participation in discussions. Thanks are also due to V. G. Mikhalkosky for his assistance in preparing the English version of the manuscript. This work was done with support from the Russian foundation for Basic Research (grant 99-05-64753) and RFBR grant of leading scientific schools of the Russian Federation No. 00-15-98509.


text:
# Kinks, energy conditions and closed timelike curves

## 1 Introduction

Chamblin has recently stressed1 that the possibility of extending a non-spacelike geodesically incomplete kinked spacetime to a spacetime that is geodesically complete hinges on the energy conditions being satisfied. On the other hand, Gibbons and Hawking proposed2 that kinked space-times should be related to the existence of closed timelike curves (CTCs). This proposal was subsequently critisized by Chamblin and Penrose3, who showed that this relation could not be established, at least for classical space-times and matter.

The above alluded connections have been already discussed in the case of de Sitter space kink4 and other spherically-symmetric gravitational topological defects5,6. It turned out that in these examples the energy conditions were not satisfied7 and, however, they all can be maximally extended to a geodesically complete spacetime. Moreover, none of these space-times may contain CTCs when their matter contents are considered classically4,8, in accordance with the Chamblin-Penrose claim3.

As a first step to investigate the influence of space-time symmetry on these issues, they will be discussed in the present paper for the case that the kink occurs in a cylindrically-symmetric space-time.

## 2 Cylindrically-symmetric kink

The typical example of a cylindrically-symmetric kink is the kink in the interior of an extreme cosmic string9,10. This is characterized by a gravitational coupling $`G\mu =\frac{1}{2}`$, where $`\mu `$ is the string linear density. In what follows, let us briefly first review the topological properties of the kinked extreme string, and then comment on some aspects of its geometry.

The general concept of a gravitational kink can be introduced by starting with the Lorentz metric $`g_{ab}`$ of a four-dimensional spacetime as given by a map, $`P`$, from any connected three-manifold, $``$M, of the spacetime four-manifold, M, into the set of timelike directions in M 11. Metric homotopy can then be classified by the degree of this map, and the kink number (or topological charge) of the Lorentz metric, with respect to a hypersurface $`\mathrm{\Sigma }`$, can be defined by11

$$\mathrm{Kink}(\mathrm{\Sigma };g_{ab})=\mathrm{deg}\left(P\right),$$

so that the gravitational kink can be viewed as a measure of how many times the light cones rotate around as one moves along hypersurface $`\mathrm{\Sigma }`$.

In the case of the spacetime of an extreme cosmic string, whose interior geometry can be visualized as that of a sphere when the corresponding two-metric is embedded in an Euclidean three-sphere12, the pair ($`\mathrm{\Sigma };g`$) will describe a gravitational kink with topological charge $`\kappa =+1`$ if $`\mathrm{Kink}(\mathrm{\Sigma };g)=1`$. From the above discussion, one may also visualize the internal geometry of the extreme string by enforcing the constant-time sections, $`\tau =\tau _0`$, of the interior metric of a string9,12 with uniform density $`ϵ`$, out to some cylindrical radius,

$$ds^2=\frac{dr^2}{1\frac{r^2}{r_{}^2}}+dz^2+r^2d\varphi ^2,$$
(2.1)
where

$$r=r_{}\mathrm{sin}\frac{\rho }{r_{}},$$
(2.2)
with $`r_{}=\left(8\pi Gϵ\right)^{\frac{1}{2}}`$, and $`\mathrm{}<z<\mathrm{},0\varphi 2\pi ,0\rho r_{}\mathrm{arccos}\left(14G\mu \right)`$, to be isometrically embedded in the kinked spacetime. The corresponding cylindrically-symmetric standard, kinked metric is given by9,13

$$ds^2=\mathrm{cos}2\alpha d\widehat{t}^22kd\widehat{t}dr+dz^2+r^2d\varphi ^2,$$
(2.3)
where the upper/lower sign of the second term corresponds to a positive/negative topological charge, $`k=\pm 1`$, depending on which of the two coordinate patches required for a complete description of the kink is being considered9, and $`\alpha `$ is the tilt angle of the light cones in the kink, $`0\alpha \pi `$. The isometric embedding will hold if in metric (2.3) we have furthermore

$$\mathrm{cos}2\alpha =1\frac{r^2}{r_{}^2}$$
(2.4)
and

$$\widehat{t}=\tau _0k\frac{dr}{\mathrm{cos}2\alpha }.$$
(2.5)
Actually, a gravitational kink depends only on D-1 of the D spacetime coordinates, and is spherically symmetric on them. However, the cylindric coordinate $`z`$ in metric (2.1) and (2.3) is not going to play any role in the analysis to follow and, therefore, one could reduce these metrics just to their hemispherical $`z`$=const. sections. On the other hand, one can also embed the $`z`$=const. sections of metric (2.3) in an Euclidean space and, hence re-express that metric in an explicit spherically-symmetric form:

$$ds^2=\mathrm{cos}2\alpha d\widehat{t}^2kd\widehat{t}dr+r_{}^2d\mathrm{\Omega }_2^2,$$

where we have specialized to the case of a gravitational kink with positive topological charge, $`\kappa =+1`$, $`d\mathrm{\Omega }_2^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$ is the metric on the unit two-sphere, and we have used Eq. (2.5).

The metric (2.3) can be then regarded as the metric for the embedding of metric (2.1), and the kinked time $`\widehat{t}`$, as the corresponding embedding function. Hence, one can obtain an embedding “rate”

$$\frac{d^2r}{d\widehat{t}^2}=\frac{2r}{r_{}}\left(\frac{r^2}{r_{}^2}1\right),$$
(2.6)
which tells us that the embedding surface would flare either outward if $`r<r_{}`$, or inward if $`r>r_{}`$. The string metric (2.1) should now be interpreted as a kinked boundary in the space with kinked spacetime (2.3).

If the isometric embedding of metric (2.1) in metric (2.3) holds, from (2.2) and (2.4) we have $`\mathrm{cos}^2\theta =\mathrm{cos}2\alpha `$, with $`\theta =\frac{\rho }{r_{}}`$, and if the one-kink is conserved, then $`G\mu `$ is enforced to be $`\frac{1}{2}`$ and $`r`$ should be analytically continued beyond $`r_{}`$, up to $`\sqrt{2}r_{}`$9. This extension creates a spherical shell filled with broken phase at each $`z`$-const. section, preventing the extreme string with $`G\mu =\frac{1}{2}`$ from disappearing, and converts the conical singularity at $`r=r_{}`$ into a de Sitter-like cosmological singularity (horizon)9. All of the topological charge of the kink would then be confined within the shell, that is within a finite compact region beyond the cosmological horizon that extends up to $`r=\sqrt{2}r_{}`$. Inside the horizon all hypersurfaces $`\mathrm{\Sigma }`$ are everywhere spacelike. Thus, as a consequence from the back reaction of the gravitational field of the one-kink, the lost picture of a cosmic string with a core region of trapped energy would be recovered for the extreme string with $`G\mu =\frac{1}{2}`$.

We have established a consistent and regular embedding of the extreme string metric in each of the two patches of the kinked spacetime whose surfaces would, according to expression (2.6), flare outward at $`\sqrt{2}r_{}`$, with a maximum “rate”

$$\frac{d^2r}{d\widehat{t}^2}|_{r=\sqrt{2}r_{}}=\frac{2\sqrt{2}}{r_{}}.$$

To stationary observers at the center of the sphere corresponding to each surface $`z`$=const., $`\tau `$=const., the compact shell containing all the topological charge of the kink14 locally coincides with a finite region of the exterior of either a de Sitter space when the light cones rotate away from the observers (positive topological charge), or the time-reverse to de Sitter space if the observers see light cones ratating in the opposite direction (negative topological charge). In the latter case, only the region outside the cosmological horizon would be accessible to stationary observers.

The kink metric (2.3) is geodesically incomplete as it shows an apparent horizon at $`r=r_{}`$ on the two coordinate patches. However, it can be converted into a geodesically complete metric by introducing Kruskal coordinates. It was obtained in Ref. 9 that the maximally-extended metric describing the spacetime of an extreme string kink can be written as

$$ds^2=\frac{4kr_{}^2}{\left(kUV\right)^2}dUdV+dz^2+r^2d\varphi ^2,$$
(2.7)
where again $`k\left(=\pm 1\right)`$ labels the two coordinate patches required to describe a complete one-kink, and $`U`$ and $`V`$ are the Kruskal coordinates9

$$U=e^{\frac{k\widehat{t}}{r_{}}}\left(\frac{r_{}r}{r+r_{}}\right),V=\pm e^{\frac{k\widehat{t}}{r_{}}},$$
(2.8)
in terms of which the radial coordinate can be defined as

$$r=r_{}\left(\frac{k+UV}{kUV}\right),$$
(2.9)
with the time $`\widehat{t}`$ given by

$$\widehat{t}=tkr_{}\sqrt{2\left(1\frac{r^2}{2r_{}^2}\right)}$$

$$+\frac{1}{2}kr_{}\mathrm{ln}\left[\frac{\left(1+\sqrt{4\left(1\frac{r^2}{2r_{}^2}\right)}\right)\left(rr_{}\right)}{\left(1\sqrt{4\left(1\frac{r^2}{2r_{}^2}\right)}\right)\left(r+r_{}\right)}\right],$$
(2.10)
where $`t`$ is the metrical kinked time which is related to the time entering the metric of the kinkless cosmic string9.

The interesting feature of metric (2.7) is that its $`z`$-const. sections are exactly the metric which describes a hemispherical section of the de Sitter spacetime kink4 for a positive cosmological constant $`\mathrm{\Lambda }=\frac{3}{r_{}^2}`$.

## 3 Energy conditions in cylindrically-symmetric kinks

The simplest, general metric describing the spacetime of a cylindrically-symmetric kink can be written as

$$ds^2=\mathrm{cos}2\alpha \left(dt^2dr^2\right)2\mathrm{sin}2\alpha dtdr+dz^2+r^2d\varphi ^2.$$
(3.1)
In order to investigate the possibility that geodesically incomplete kinks can be extended to geodesically complete ones when the energy conditions are satisfied in a cylindrically-symmetric space-time, let us review these conditions and also physical conditions of the kind considered by Finkelstein and McCollum13, by using the Hawking-Ellis procedure15. Thus, we write the eigenvalue equation

$$\left(G_{\mu \nu }\lambda g_{\mu \nu }\right)\xi ^\nu =0.$$
(3.2)
For this equation to be implemented in the case of a cylindrically-symmetric space-time kink, we obtain first the nonzero Chrystoffel symbols for metrics (3.1). These are:

$$\mathrm{\Gamma }_{tt}^t=\mathrm{sin}^22\alpha _r\alpha ,\mathrm{\Gamma }_{tr}^t=\mathrm{\Gamma }_{tt}^r=\mathrm{sin}2\alpha \mathrm{cos}2\alpha _r\alpha ,$$

$$\mathrm{\Gamma }_{rr}^t=\left(1+\mathrm{cos}^22\alpha \right)_r\alpha ,\mathrm{\Gamma }_{tr}^r=\mathrm{sin}^22\alpha _r\alpha ,$$

$$\mathrm{\Gamma }_{\varphi \varphi }^t=r\mathrm{sin}2\alpha ,$$
(3.3)
$$\mathrm{\Gamma }_{rr}^r=\mathrm{sin}2\alpha \mathrm{cos}2\alpha _r\alpha ,\mathrm{\Gamma }_{r\varphi }^\varphi =r^1,\mathrm{\Gamma }_{\varphi \varphi }^r=r\mathrm{cos}2\alpha .$$

Hence we can obtain the nonzero components of the Ricci tensor which are:

$$R_{tt}=R_{rr}=\frac{1}{2r}\mathrm{cos}2\alpha _r\left(r^2\mathrm{}\right),$$

$$R_{tr}=\frac{1}{2r}\mathrm{sin}2\alpha _r\left(r^2\mathrm{}\right),$$
(3.4)
$$R_{\varphi \varphi }=r^2\mathrm{},$$

where

$$\mathrm{}=\frac{2}{r^2}_r\left(r\mathrm{sin}^2\alpha \right)=\frac{2}{r^2}_r\mu ,$$
(3.5)
with $`\mu `$ ($`=r\mathrm{sin}^2\alpha `$) the function introduced by Finkelstein and McCollum13. For the curvature scalar we then obtain

$$R=\mathrm{}+\frac{1}{r}_r\left(r^2\mathrm{}\right),$$
(3.6)
and for the nonzero mixed components of the Einstein tensor

$$G_t^t=G_r^r=\frac{1}{2}\mathrm{}=\frac{1}{r^2}_r\mu $$

$$G_z^z=\frac{1}{r^2}_r\mu \frac{1}{r}_r^2\mu $$
(3.7)
$$G_\varphi ^\varphi =\frac{1}{r^2}_r\mu \frac{1}{r}_r^2\mu .$$

Introducing these Einstein-tensor components in the eigenvalue equation (3.2), we obtain the eigenvalues

$$\lambda _0=\lambda _1=\frac{1}{r^2}_r\mu $$

$$\lambda _2=\lambda _1=\frac{1}{r^2}_r\mu \frac{1}{r}_r^2\mu $$
(3.8)
$$\lambda _3=\lambda _1=\frac{1}{r^2}_r\mu \frac{1}{r}_r^2\mu .$$

The corresponding eigenvectors are

$$E_0=(\mathrm{cos}\alpha ,\mathrm{sin}\alpha ,0,0)$$

$$E_1=(\mathrm{sin}\alpha ,\mathrm{cos}\alpha ,0,0)$$

$$E_2=(0,0,1,0)$$
(3.9)
$$E_3=(0,0,0,\frac{1}{r}).$$

Since $`E_0`$ is timelike and the $`E_\rho `$’s ($`\rho =1,2,3`$) are all spacelike, we have a canonical form of Type I, according to the classification of Hawking and Ellis15. The spacelike eigenvectors $`\left\{E_\rho \right\}`$ form an orthonormal basis and the tetrad components of the metric tensor are

$$\overline{g}_{\rho \sigma }=g(E_\rho ,E_\sigma )=\mathrm{diag}(1,1,1,1).$$
(3.10)
Finally, we obtain for the tetrad components of the energy-momentum tensor

$$\overline{T}^{\rho \sigma }=\mathrm{diag}(ϵ,p_1,p_2,p_3),$$
(3.11)
with

$$ϵ=\frac{1}{r^2}_r\mu ,p_1=\frac{1}{r^2}_r\mu ,$$
(3.12)
$$p_2=\frac{1}{r^2}_r\mu \frac{1}{r}_r^2\mu ,$$
(3.13)
$$p_3=\frac{1}{r^2}_r\mu \frac{1}{r}_r^2\mu .$$
(3.14)
We can now discuss the energy conditions in general cylindrically-symmetric kinked space-times. The weak energy condition states15,16 that the energy density as measured by an observer must be non-negative, and this requires

$$ϵ0,ϵ+p_\rho 0,\rho =1,2,3.$$
(3.15)
Using (3.12)-(3.14), inequalities (3.15) lead to:

$$_r\mu 0$$
(3.16)
$$_r\left(\frac{1}{r^2}_r\mu \right)0$$
(3.17)
$$_r^2\mu 0.$$
(3.18)
On the other hand, for a canonical form of type I there also holds the strong energy condition15,16, provided that

$$ϵ+p_\rho 0,ϵ+\underset{\rho }{}p_\rho 0,\rho =1,2,3.$$
(3.19)
For (3.12)-(3.14) these inequalities lead to (3.18) again. We note that although the strong energy condition does not imply the weak energy condition, the vice versa is true however.

Physical conditions which function $`\mu \left(r\right)`$ must satisfy7 are: (a) $`_r\mu 0`$ for all $`r`$, (b) $`_r\left(\frac{1}{r^2}_r\mu \right)0`$ for all $`r`$, (c) $`_r^2\mu 0`$ for all $`r`$, (d) $`\mu =O\left(r^3\right)`$ as $`r0`$ in order for $`\left|G_\rho ^\sigma \right|<\mathrm{}`$ at $`r=0`$, and (e) $`0\frac{\mu }{r}1`$ in order for $`0\mathrm{sin}^2\alpha 1`$. On the other hand, to ensure the existence of an one-kink, we need to impose the boundary conditions $`\alpha \left(0\right)=0`$ and $`\alpha \left(0\right)=\pi `$, since in the present case, $`\alpha =\mathrm{arcsin}\left(\frac{r}{\sqrt{2}r_{}}\right)`$. Hence, recalling that $`\frac{\mu }{r}=\mathrm{sin}^2\alpha `$, we have the additional condition (f) $`lim_{r0}\frac{\mu }{r}=0`$. From this condition it follows $`\mu =O\left(r^{1+c}\right)`$, with $`c>0`$ for small $`r`$, and hence $`_r\mu =O\left(r^c\right)>0`$ and $`_r^2\mu =O\left(r^{c1}\right)>0`$, for small $`r>0`$. The latter inequality violates condition (c) and, therefore, not only the strong energy condition but also the weak energy condition is violated in this space-time. Thus, like for kinked de Sitter space, space-time kinks possessing cylindrical symmetry should violate the above two energy conditions, in spite of they being maximally extendible to geodesically complete space-time. Moreover, it can be readily seen that these violations also imply violation of the dominant energy condition16, $`ϵ0`$, $`p_\rho [ϵ,+ϵ]`$, and the null energy condition16 $`ϵ+p_\rho 0`$.

Let us finally consider what happens to the energy conditions when we restrict ourselves to a $`z=`$const. section of metric (3.1). It is easy to see then that the weak energy condition implies only inequalities (3.16) and (3.17) to hold, while the strong energy condition leads to (3.17) again. It follows that the former condition would only be satisfied if $`0c2`$, while the latter one holds when $`0c1`$. Since in the actual case $`c=2`$, we see that for $`z`$=const. sections of metric (3.1), the weak energy condition is satisfied but the strong energy condition remains being violated. This conclusion can be extended to the geodesically complete space-times. Thus, the Chamblin’s argument1 for cylindrically-symmetric kinks should be relaxed so that a non-spacelike geodesically incomplete kinked spacetime would be extendible to a geodesically complete one if the weak (but not the strong) energy condition is preserved on their $`z=`$const. sections.

## 4 Kinks and closed timelike curves

We shall analyse next the possible connection between kinks and CTCs, suggested by Gibbons and Hawking2, for the case of an extreme cosmic-string kink. Both the discussion and conclusion to be obtained can straightforwardly be generalized to any other spacetime kinks.

The existence of CTCs in a given spacetime can be investigated by considering the paths followed by null geodesics on the Kruskal diagrams. For the extreme string kink with Kruskal diagrams given in Fig. 1, it might seem at first glance that since the two coordinate patches are identified on the surfaces $`r=A=\sqrt{2}r_{}`$, both on the original regions I and II and the new regions III and IV (created by Kruskal extension), because of continuity of the tilt angle at $`\alpha =\frac{\pi }{2}`$, one could choose null geodesics that started at $`r=0`$ in region I+ and would somewhat loop back through the new regions to finally arrive at their starting point, after traversing both coordinate patches. However, one can easily convince oneself that such itineraries are classically disallowed, since they would require identification of the two patches also on minimal surfaces at $`r=0`$ belonging to a physical and a nonphysical region, respectively (see Fig. 1).

Nevertheless, the above conclusion is no longer valid when we consider propagation of quantum fields in the same kinked spacetime. The semiclassical regime for kinked spacetimes with event horizons can simply be achieved by considering the mathematical implications imposed by the fact that time $`\overline{t}`$ enters the Kruskal coordinates $`U,V`$ in the form of the dimensionless exponent $`k\widehat{t}/r_{}`$. The argument of the logarithm in (2.10) becomes then square rooted and, therefore, the expression for the time $`\widehat{t}`$ entering the definition of coordinates $`U,V`$ should be generalized to9:

$$\widehat{t}\widehat{t}_g=\widehat{t}+\frac{i}{2}k\kappa \left(1\kappa \right)\pi r_{},$$
(4.1)
where $`\widehat{t}`$ is given by (2.10) and $`\kappa =\pm 1`$. For $`\kappa =+1`$, $`\widehat{t}_g=\widehat{t}`$, and for $`\kappa =1`$, the points ($`\widehat{t}ik\pi r_{},r,z,\varphi `$) in each patch are actually the points in the same patch obtained by reflection in the bifurcation point $`U,V=0`$, keeping the Kruskal metric real and unchanged.

We note that one can still recover the standard kink metric (2.3) from a general Kruskal metric if we redefine the Kruskal coordinates as follows:

$$U=\pm 2b\kappa e^{\frac{k\widehat{t}_c}{r_{}}}\left(\frac{r_{}r}{r_{}+r}\right)$$
(4.2)
$$V=k\kappa e^{\frac{k\widehat{t}_c}{r_{}}},$$
(4.3)
where

$$\widehat{t}_c=\widehat{t}+i\pi r_{}\kappa k.$$
(4.4)
This choice leaves the expressions for $`UV,F,r`$ and the Kruskal metric real and unchanged. For $`\kappa =1`$, Eqns. (2.8) become the sign-reversed of (4.2) and (4.3), respectively; i.e.: the points ($`\widehat{t}ik\pi r_{},r,z,\varphi `$) on the Kruskal diagrams of the two coordinate patches are the points in the new regions IIIk on the same diagrams, obtained by reflecting in the origins of the respective $`U,V`$ planes, preserving the Kruskal metric real and unchanged. This leads to identification of hyperbolae in the new regions with hyperbolae in the original regions for the same values of $`r`$; i.e. to identification of hyperbolas IIIk with hyperbolas IIk and hyperbolas IVk with hyperbolas Ik. We note that the existence of such identifications in turn amounts to both the kind of periodicity required by Hawking thermal radiation17 in each patch, and the existence of CTCs.

Since in the semiclassical description, we can identify maximum surfaces of the physical regions with those of the nonphysical regions in each coordinate patch separately, we recover allowance for the null geodesics that start at surface $`r=0`$ in region I+ of patch $`k=+1`$ (Fig. 1) to continue propagating on patch $`k=1`$, after the maximum surface of region III-, first through region II- and then through region IV-, up to the surface at $`r=A=\sqrt{2}r_{}`$ of the latter new region. Because this surface can be identified with the similar surface in region III+ of patch $`k=+1`$, the considered null geodesics can thereafter propagate into the region IV+ and, again by quantum identification of surfaces at $`r=0`$, come back to their starting points on the surface at $`r=0`$ of region I+, in patch $`k=+1`$. Hence, null geodesics starting from original regions at $`r=0`$ can still loop back to arrive at their starting points, so completing a CTC, provided such a CTC is involved at a thermal radiation process preventing any information to flow from or to the CTC. Our general conclusion therefore is that CTCs are linked to spacetime kinks if and only if the matter traveling through these spacetimes is considered quantum-mechanically.

## 5 Conclusion

The problem of the relation between geodesic incompleteness and the holding of the energy conditions has been considered for the case of a kink in a cylindrically-symmetric spacetime. We discussed the physical conditions that such a spacetime must satisfy in relation with the weak and strong energy conditions. It was shown that, whereas none of these conditions holds for the four-dimensional case, the weak energy condition can still survive on the $`z`$ = const. (transversal) sections of the spacetime. Also considered has been the problem of the connection between spacetime kinks and the existence of closed timelike curves. We obtained that kinked spacetimes do not contain CTCs if the matter propagating on them is dealt with classically, but CTCs become a necessary ingredient in such spacetimes whenever the propagating matter shows quantum behaviour. As it was pointed out first in Ref.  and later in Ref.  there must be a close connection between CTCs and the thermal processes induced by the presence of an event horizon. Any possible observer of the quantum matter trapped in a CTC could only detect it by means of a totally incoherent radiation carrying no desciphrable information to or from the observer. Thus, the connection of CTCs with the thermal process prevents the existence of any observable causality violating process induced by the CTCs and, therefore, allows one to conjecture a censorship for causality violation, even when CTCs and time machines can exist and be operative. One might illustrate the resulting situation by re-paraphrasing Stephen Hawking : there could perfectly be hords of tourists visiting us from the future, but neither they nor we could know anything about their trip. For them, it would be a touring which costs a lot and rewards nothing; for us, the trip would simply be unnoticed.

Acknowledgements

The author acknowledges DGICYT by support under Research Project No. PB97-1218

## References

A. Chamblin, J. Geom. Phys. 13, 357 (1994).

G.W. Gibbons and S.W. Hawking, Phys. Rev. Lett. 69, 1719 (1992).

A. Chamblin and R. Penrose, Twistor Newsletter 34, 13 (1992).

K.A. Dunn, T.A. Harriott and J.G. Williams, J. Math. Phys. 35, 4145 (1994).

P.F. González-Díaz, Phys. Rev. D51, 7144 (1995).

P.F. González-Díaz, in: Theories of Fundamental Interactions; Maynooth Bicentenary Volume, ed. T. Tchrakian (World Scientific, Singapore, 1995); Grav. Cosm. 2, 621 (1996).

K.A. Dunn, T.A. Harriott and J.G. Williams, J. Math. Phys. 37, 5637 (1996).

P.F. González-Díaz, Int. J. Mod. Phys. D7, 793 (1998).

P.F. González-Díaz, Phys. Rev. D52, 5698 (1995).

A. Vilenkin and E.P.S. Shellard, Cosmic Strings and other Topological Defects (Cambridge U. P., Cambridge, UK, 1994).

G.W. Gibbons and S.W. Hawking, Commun. Math. Phys. 142, 1 (1990).

J.R. Gott, Astrophys. J. 288, 422 (1985); W.A. Hiscock, Phys. Rev. D31, 3288 (1985).

D. Finkelstein and C.W. Misner, Ann. Phys. (N.Y.) 6, 230 (1959).

A. Chamblin, Kinks and singularities, Preprint, DAMTP-R95/44 (1995).

S.W. Hawking and G.F.R. Ellis, The Large Scale Structure of Space-Time (Cambridge U. P., Cambridge, UK, 1973).

M. Visser, Lorentzian Wormholes (American I. P. P., Woodbury, USA, 1996).

J.B. Hartle and S.W. Hawking, Phys. Rev. D31, 2188 (1976).

P.F. González-Díaz, Time Walk through the quantum cosmic string, gr-qc/9412034.

S.W. Hawking, Phys. Rev. D52, 5681 (1995).

S.W. Hawking, Phys. Rev. D46, 603 (1992).


text:
# Internal Motions in Globular Clusters

## 1. Introduction

The dynamics of globular clusters is a venerable field, but new advances continue to be made. In this brief review we will deal with vistas that have opened up as a result of work in which we are currently engaged. Although we will deal with some theory, our discussion will be mainly observational.

## 2. HST astrometry

We began trying to measure proper motions with the Hubble Space Telescope (HST) in order to separate cluster stars from field stars in NGC 6397, where the cluster main sequence gets lost among the field stars a magnitude and a half above our faint limit. With WFPC2 images spanning less than 3 years we found that we were able to effect an excellent separation, and thereby detect the lowest part of the main sequence (King et al. 1998).

But the real turning point came when Georges Meylan asked us to measure his WFPC2 images of the center of 47 Tucanae, taken two years apart. His original aim was merely to detect which stars had the largest internal proper motions in the cluster, but we quickly found that our techniques could produce a motion for each individual star. After two more years of work on our methods, we have gained an improvement of nearly an order of magnitude in our accuracy, by identifying and removing the sources of systematic error that arise in the severely undersampled WFPC2 camera. The most important step, we have found, is to derive an extremely accurate point-spread function (Anderson & King 2000).

In 47 Tuc, where the dispersion of proper motions is about 0.6 milli-arcsec/yr in each coordinate, we are now in process of using 1995, 1997, and now 1999 images, to measure proper motions of individual stars at about an $`8\sigma `$ level of accuracy or better. At the time of this conference we have derived proper motions only from the 1995 and 1997 images; they are at about the $`4\sigma `$ level. Figure 1 shows examples of our measurements and their scatter about the mean positions that we derive from them.

The ability to measure proper motions of individual stars in globular clusters opens up many new avenues of research, several of which we will be discussing here. We will examine the shape of the (isotropic) velocity distribution at the cluster center, point out the value of proper motions in measuring cluster distances, and, from another data set, examine the anisotropy of the velocity distribution away from the center.

### 2.1. Velocity distribution in the cluster center

At the center of the cluster there is no anisotropy in the distribution of stellar velocities; as we shall indicate below, anisotropy manifests itself only in the outer parts of a cluster. What we can study quite well here is the form of the velocity distribution. We have motions for enough stars to make a good comparison with the distribution that should be expected. We show such a comparison in Figure 2. Most cluster models have used a “lowered Gaussian”, i.e., a Gaussian minus a constant, where the constant is chosen so as to place the cutoff at the proper escape velocity for the cluster center; we have made such a choice here. The only free parameter in the fit is the velocity dispersion of the Gaussian. Because even the 490 stars used here would give rather noisy numbers if binned into a probability-distribution histogram, we have chosen instead to exhibit a cumulative-distribution step-function, which rises to the left because we have accumulated from high velocities toward low.

The quantity that we have chosen to compare is the magnitude of the proper motion. The fit is strikingly good; the deviations near the high-velocity end are probably due to small-number statistics.

Interestingly, we have also been able to examine the “cannonball stars” (Meylan, Dubath, & Mayor 1991) whose existence was part of the original motivation for measuring proper motions in 47 Tucanae. They turn out to have nearly zero proper motions, so that their radial velocities are the whole of their motions, which do not appear to be above the escape velocity at the center of the cluster.

Another distribution that we can examine is that of the blue stragglers. They are known to be more concentrated to the cluster center than are the red giants; they should therefore have a lower velocity dispersion. Their distribution in Fig. 2 confirms this quite well.

### 2.2. Cluster distances

Although it is outside the immediate focus of this conference, we should mention that the most valuable product of such proper motions will be fundamental distances of globular clusters, found by comparing the dispersion of proper motions, an angular quantity, with that of radial velocities, a linear quantity. With 10,000–20,000 proper motions in a cluster and at least 2500 radial velocities measured by our collaborators (see, e.g., Gebhardt et al. 1997), the only completely unavoidable error will be that due to statistical sample size, which is between 1 and 2%. Our aim, then, will be to avoid all other sources of error. The errors of measurement are not a problem per se; what counts is that their size be accurately known, so that they can be correctly subtracted in quadrature in order to derive the true dispersions.

A serious concern will be the effect of orbital motion in binaries; it must be corrected for statistically, and will be different for the radial velocities and the proper motions, because of the different distribution of the observations in time. To make it worse, the binary fraction in globular clusters is poorly known.

Another problem of consequence takes us back to dynamics: the radial velocities of 47 Tuc stars show a quite prominent rotation of the cluster (Meylan & Mayor 1986). We can in fact see rotation clearly in the proper motions too; at the cluster center they show a slightly greater dispersion in the equatorial than in the polar direction. It is obviously essential that we have a good dynamical model for the cluster. (It should also go without saying that the model should have the correct mass function for this cluster, a precept that has not always been followed. Fortunately our HST observations of the cluster provide a good mass function down to quite low masses.) For the modeling itself we intend to generalize the rotating models of Wilson (1975) to the multi-mass case.

## 3. The theory of anisotropy in globular clusters

### 3.1. The origin of anisotropy

The most simple-minded dynamical models of clusters, such as the original “King models” (King 1966), have isotropic velocity distributions, but this is a serious oversimplification for real clusters. Whether a cluster is born by the collapse of a larger configuration or by expansion of a smaller one, the initial motions are largely radial; and conservation of angular momentum in each stellar orbit will then constrain the motions to remain largely radial. This initial anisotropy should be removed gradually as relaxation moves the velocity distribution toward isotropy, and this trend should spread gradually outward from the dense center to the lower-density envelope.

Intuition suggests that the velocity distribution in a cluster should have retained its anisotropy outside the radius where the local relaxation time equals the cluster age (i.e., about a Hubble time), and this is approximately what is observed.

### 3.2. The nature of anisotropy

In a spherically symmetric globular cluster, the three phase-space variables on which the velocity distribution can depend are the radial distance $`r`$, the radial component of velocity $`v_\mathrm{r}`$, and the tangential component of velocity $`v_\mathrm{t}`$. It has been customary to represent anisotropy in a globular cluster by writing the velocity distribution as

$$f(r,v_\mathrm{r},v_\mathrm{t})=F(E+cJ^2),$$
(1)
where

$$E=\frac{1}{2}\left(v_\mathrm{r}^2+v_\mathrm{t}^2\right)+U(r)$$
(2)
is the energy per unit mass ($`U`$ being, of course, the gravitational potential function), and

$$J=rv_\mathrm{t}$$
(3)
is the magnitude of the angular momentum per unit mass. The reasons for this representation are that (1) Jeans’ theorem says that in a time-independent state (which is a very good approximation in a slowly evolving system) the velocity distribution must be a function of these two integrals of the orbital motion of a star, and (2) this simple form is easily shown to lead to a velocity spheroid whose radial elongation increases with increasing $`r`$.

It is this kind of velocity distribution that Gunn and Griffin (1979) and Meylan (1987, 1988) used in their studies of anisotropy, and a modified form of this distribution that Lupton, Gunn, & Griffin (1987) used in their treatment of anisotropy in a rotating cluster. We will show below, however, that observed anisotropic velocity distributions unfortunately do not have this form.

### 3.3. Dependence of anisotropy on stellar mass

Nothing is known about how anisotropy differs from one stellar mass to another. We are not aware of any theoretical prediction, nor is there observational information, since existing studies, through radial velocities, have been confined to stars above the main-sequence turnoff, all of which have practically the same mass. We expect to derive such information from observations that will be described below, but for the time being one can only make a rough theoretical estimate.

Anisotropy is removed by relaxation, whose rate goes as $`v^3`$, where $`v`$ is a characteristic velocity for the stars in question. Since in the condition of equipartition that seems to exist in nearly all globular clusters $`v1/m`$, the time to remove anisotropy would at first appear to go as $`m^{3/2}`$. But the dependence on mass should be even stronger, since the stars of lower mass inhabit regions of lower density, on account of their higher velocities, and this makes their relaxation even slower. A quantitative study, however, promises to be rather complicated, and we will not pursue it here.

## 4. Methods of observing anisotropy

Until now, anisotropy of stellar motions in globular clusters has been studied only indirectly, by comparing the radial drop-off that is observed in the dispersion of radial velocities with the drop-off that is predicted by theoretical models with differing amounts of anisotropy. Not only is this method indirect; it is also too vulnerable to inadequacies in the dynamical modeling. By contrast, proper-motion measurements measure anisotropy directly, and are nearly independent of any modeling.

The analysis of the radial-velocity observations depends on our different geometrical perspective at the center and at the edge of the cluster. The velocities that we observe along the line of sight at the cluster center are clearly values of $`v_\mathrm{r}`$ only, whereas at points away from the center the observed dispersion has an increasing admixture of $`v_\mathrm{t}`$, so that the drop-off of $`v^2_{\mathrm{los}}`$ gives us information about the anisotropy. Unfortunately that information is only indirect. The problem is that there is a natural drop-off of $`v^2_{\mathrm{los}}`$ even in a cluster whose velocity distribution is isotropic; because the velocity distribution has a cutoff at the escape velocity, the decrease of the latter with increasing radius causes $`v^2_{\mathrm{los}}`$ to fall. Thus the observed quantity that is needed becomes in effect a higher-order one—the difference between the actual fall-off of $`v^2_{\mathrm{los}}`$ and the fall-off that would occur in a cluster model that has isotropic velocities. This weakness gives radial velocities only a poor grip on anisotropy. One should note further that such a measurement of anisotropy is very sensitive to the cluster model adopted. The results will suffer badly from any errors either in the assumed mass function or in the functional form that is used for the representation of anisotropy in the velocity distribution. In the latter respect, we will show below that the functional form that has commonly been used is contradicted by our new observations.

The proper motions play a quite different role. At a point away from the cluster center one can directly compare the dispersion in the radial direction with that in the transverse direction. There are, to be sure, projection problems along the line of sight (as there are equally in the interpretation of the radial velocities), but the proper motions give us our only really direct measure of anisotropy.

A caution to note, however, is that rotation of the cluster will distort the apparent anisotropy. Along a given line of sight, components of the rotational motion will be included in the part of the proper motions that is parallel to the cluster equator. (This is similar to the small rotation-induced anisotropy that we have already noted at the cluster center.) Interpretation of the proper motions thus has a small dependence on modeling of the cluster, but far less than in the case of the radial velocities.

## 5. Our new observations: anisotropy in 47 Tuc

In addition to our work on Meylan’s images at the center of 47 Tucanae, we have a program of our own, in which we have images of two fields on opposite sides of the center of 47 Tuc, each at about 10 core radii from the center. Meylan’s radial-velocity study (1988) suggests that anisotropy becomes strong at a radius somewhere between there and 40 core radii, so our fields should show a modest amount of anisotropy. Our time baselines are 4–5 years. As already indicated, we are able to examine anisotropy directly, by comparing proper-motion dispersions in the radial and transverse directions. We have so far made a preliminary reduction of one of these fields. Although our measurements are still only preliminary, their trend is clear.

The results are surprising. We do see anisotropy, but only among the stars of highest velocity. In other words, in our two-dimensional distribution of proper motions the lines of constant density are not ellipses of constant axial ratio. The contours near the center are circles, while only the largest contours show an elongation in the radial direction. These observations contradict the picture of similar-shaped ellipses that would follow from the sort of velocity distribution that is specified by Eq. (1).

This is a result that we had not foreseen, and as far as we know no one else had either; but it should not have been a surprise, because in clear hindsight it is exactly what we should have expected. To understand it we must keep in mind that anisotropy is removed only by relaxation, whose rate is proportional to the star density in the regions through which a star passes. In a region like the one that we have observed, the stars of low velocity have orbits that bring them into the denser inner regions, where they have seen enough relaxation to isotropize their orbits, whereas those of higher velocity have orbits in which they have spent most of their lives in the low-density outer regions, and with less exposure to the relaxing tendency they have retained some of their anisotropy.

Given the fact that anisotropic velocity distributions are not as we thought they were, how is the theory of cluster dynamics to react? Obviously we need to construct cluster models based on velocity distributions in which the roles of $`E`$ and $`J`$ differ from the simplistic ones in Eq. (1). Unfortunately it is not yet clear what sort of mathematical function will fill the bill. An important challenge to the building of realistic dynamical models of star clusters is to find a suitable function. Presumably, for the sake of providing a finite escape velocity it can then be cast in a form analogous to the

$$f(E,J^2)=e^{2j^2(E+aJ^2)}e^{2j^2aJ^2},$$
(4)
which is the form used in the now-outmoded models.

We look forward to a better understanding of anisotropy after we have fully reduced the images of these two fields in 47 Tuc, and especially after measuring an outer field in $`\omega `$ Centauri, where from Meylan’s (1987) results we expect the anisotropy to be much greater.

The interpretation of anisotropy observations will also involve modeling of the cluster. In this case we will need models that include both anisotropy and rotation. Fortunately Lupton & Gunn (1987) have shown how to build such models.

## 6. Summary

Observations of internal motions in globular clusters are obviously essential for understanding the dynamics of these systems. Radial velocities have contributed some information, but the new ability to measure internal proper motions of individual stars, using the Hubble Space Telescope, opens new vistas for observational studies of dynamics.

Proper motions of large numbers of stars show that the velocity distribution is quite close to Gaussian, as dynamical models have supposed.

Proper motions are much more suitable than radial velocities for the study of anisotropy of velocity dispersions in the outer parts of clusters. Here preliminary results indicate that the anisotropy does not take the form that is customarily assumed in dynamical models; hindsight easily shows why this should be so. New mathematical forms for the velocity distribution need to be contrived.

The most important application of internal proper motions in globular clusters will be the derivation of fundamental distances, from a comparison of the dispersion of the proper motions with that of the radial velocities.

Finally, we note that in a later paper in this volume Ken Freeman reports results of a unique ground-based study that analyzes internal proper motions in $`\omega `$ Centauri.

#### Acknowledgments.

The research supported here was supported by grant GO-7503 from the Space Telescope Science Institute.

## References

Anderson, J., & King, I. R. 2000, PASP, in press

Gebhardt, K., Pryor, C., Williams, T. B., Hesser, J. E., & Stetson, P. B. 1997, AJ, 113, 1026

Gunn, J. E., & Griffin, R. F. 1979, AJ, 84, 752

King, I. R. 1966, AJ, 71, 66

King, I. R., Anderson, J., Cool, A. M., & Piotto, G. 1998, ApJ, 492, L37

Lupton, R., & Gunn, J. 1987, AJ, 93, 1106

Lupton, R., Gunn, J., & Griffin, R. 1987, AJ, 93, 1114

Meylan, G. 1987, A&A, 184, 144

Meylan, G. 1988, A&A, 191, 215

Meylan, G., & Mayor, M. 1986, A&A, 166, 122

Meylan, G., Dubath, P., & Mayor, M. 1991, ApJ, 383, 587

Wilson, C. P. 1975, AJ, 80, 175