text:
# Phase transitions towards frequency entrainment in large oscillator lattices

## Abstract

We investigate phase transitions towards frequency entrainment in large, locally coupled networks of limit cycle oscillators. Specifically, we simulate two-dimensional lattices of pulse-coupled oscillators with random natural frequencies, resembling pacemaker cells in the heart. As coupling increases, the system seems to undergo two phase transitions in the thermodynamic limit. At the first, the largest cluster of frequency entrained oscillators becomes macroscopic. At the second, global entrainment settles. Between the two transitions, the system has features indicating self-organized criticality.

Many systems in science, engineering and social life can be described as large networks of coupled limit cycle oscillators pikovsky . Most often there is a spread in the individual, natural frequencies, and the coupling is such that it tends to even out these frequency differences. A general question is how the dynamics of such systems changes when the coupling between the oscillators increases. Is there a phase transition at which the oscillators attain a common, collective frequency, in the thermodynamic limit where the number $`N`$ of oscillators goes to infinity?

Such phase transitions are relevant in several fields. For example, the proper function of the millions of pacemaker cells in the sinus node in the heart, requires that they work at the same frequency. Cardiac arrhythmias may result if this is not the case sicksinus . The appearance of several frequencies in the sinus node may be caused by decoupling due to tissue degeneration. The brain also contains many pacemaker cells. Increasing evidence suggests that enhanced electric coupling between neighbor neurons can provoke epileptic seizures carlen . These correspond to pathologically large regions of synchronized electric discharges. Thus, in this case a phase transition to a state of collective oscillation is unfavorable, while it is vital for the coordinated function of the heart.

Most theoretical studies of such transitions assume that each oscillator is coupled to all the others equally strong. The most well-known system of this kind is the Kuramoto model kurarev . More realistic networks have local coupling. Winfree win1  and Kuramoto kura  hypothesized that in systems with nearest neighbor coupling and random natural frequencies, there should be a critical coupling strength at which the number of members, or size $`S_{max}`$, of the largest cluster of frequency entrained ent  oscillators becomes macroscopic. This can be expressed as a phase transition at which the order parameter $`r`$ becomes nonzero, where

$$r\underset{N\mathrm{}}{lim}S_{max}/N.$$
(1)
For a long time, model studies only revealed negative or inconclusive results regarding the existence of such a phase transition prev . Recently, we proved that it was present in a one-dimensional chain of pulse-coupled oscillators chain . A finite critical coupling strength $`g_c`$ was found, at which global frequency entrainment settles. At $`g_c`$, $`r`$ jumps discontinuously from zero to one. In this Letter, a two-dimensional square lattice with bidirectional nearest neighbor coupling is studied. Two phase transitions at the critical couplings $`g_{c1}`$ and $`g_{c2}`$ are found, at which the order parameter $`r`$ seems to become non-zero and one, respectively.

In our model (c.f. Ref. chain ), the state of oscillator $`k`$ is given by the phase $`\varphi _k[0,1)`$. The time evolution of the phase is given by

$$\dot{\varphi }_k=1/P_k+gh(\varphi _k)\underset{ln_k}{}\delta (\varphi _l).$$
(2)
$`P_k`$ is the natural period of oscillator $`k`$, and $`n_k`$ is the set of its nearest neighbors. An oscillator $`l`$ is said to fire when $`\varphi _l=1`$. Then $`\varphi _l0`$ and a pulse is delivered to the neighbor $`k`$, so that its phase immediately shifts according to $`\varphi _k\varphi _k+gh(\varphi _k)`$. This kind of system can model oscillators that interact with short pulses and are strongly attracted to their limit cycles. Examples include pacemaker cells in the heart, neurons, flashing fireflies, chirping crickets, and people clapping their hands in the theater. The function $`gh(\varphi _k)`$ is called the phase response curve (PRC), where $`g`$ is the coupling strength. Inspired by experiments on pacemaker cells in the heart sano , we assume the form of the PRC given in Fig. 1. This coupling tends to even out phase and frequency differences between oscillators for the following reason: If the phase $`\varphi _k`$ of an oscillator $`k`$ receiving a pulse from a neighbor $`l`$ is small, it becomes even smaller ($`h<0`$), approaching $`\varphi _l=0`$. If $`\varphi _k`$ is large, it becomes even larger ($`h>0`$), again approaching $`\varphi _l`$. We expect that it is this bipolar character of the PRC that is essential for the observed dynamics, not its exact shape.

The natural periods $`P_k`$ are taken as random numbers from a square distribution with $`P_{min}=1`$ and $`P_{max}=1.5`$ time units (t.u.). We use periodic boundary conditions. The same method of numerical integration as in Ref. chain  is used. The lattice is divided into blocks of $`10\times 10`$ oscillators, within which the integration is exact. These correspond to the segments of $`25`$ oscillators within which the integration was exact in Ref. chain .

The two phase transitions separate three phases. We shall say that states with $`r=0`$ belong to “phase 1”, states with $`0<r<1`$ to “phase 2”, and states with $`r=1`$ to “phase 3”. Figure 2 shows mean periods in a large lattice measured during $`10^4`$ t.u. after a transient of $`10^5`$ t.u. for different coupling strengths $`g`$. For $`g=0`$ the mean periods are the natural periods, which are independent random numbers. As we increase $`g`$, clusters of oscillators with nearly identical mean periods appear, as can be seen for $`g=0.40`$. The typical size of these clusters increases with $`g`$. At $`g=0.51`$, one cluster is almost percolating horizontally, indicating that this coupling is close to the critical value $`g_{c1}`$, which separates phase 1 from phase 2. At $`g=0.53`$ and $`g=0.54`$, one cluster is percolating through the lattice. This cluster can be interpreted as macroscopic, suggesting that we have entered phase 2. For $`g=0.55`$ the entire lattice attains the same frequency, suggesting that we have entered phase 3.

Around $`g=0.48`$, oscillators that do not fire during the measurement interval start to appear. They become more frequent as $`g`$ increases towards $`g_{c2}`$. However, their number decreases with increasing measurement interval, showing that they are not silent forever. In phase 3, however, there is a small number of oscillators ($`<0.05\%`$) that never fire. The others always fire with a common interval. Evidently, the silent oscillators are repeatedly perturbed by their neighbors at a phase $`\varphi `$ where the PRC is negative, so that the phase is kept in this interval, and does not reach the threshold $`\varphi =1`$. This situation can only be maintained at all times if the surrounding is stable, as it is in phase 3. To say that $`r=1`$ in phase 3, we have to exclude the silent oscillators.

To locate $`g_{c1}`$ more precisely, we study how the distribution of cluster sizes depends on $`g`$. A cluster is numerically defined to be a connected set of oscillators whose mean periods do not differ more than $`dP=0.001`$. The results are shown in Fig. 3(a). The distribution becomes critical (i.e. a power law) at $`g0.51`$, suggesting that $`g_{c1}`$ is close to $`0.51`$. At $`g=0.53`$ and $`g=0.54`$ the cluster sizes are still close to critically distributed, if we exclude the largest clusters, which are too large to fit in. For $`g=0.53`$ there are two such clusters, and for $`g=0.54`$ there is one. Thus the system might be described as critical in the entire phase 2. A cluster percolating through the lattice is not always seen in phase 2. However, the largest clusters are always larger than they would be if all cluster sizes were critically distributed. In some cases, a percolating cluster develops if the simulation is allowed to proceed further. The opposite is never seen, i.e. that a percolating cluster disappears, suggesting that a macroscopic, percolating cluster always appears sooner or later in phase 2. The critical distribution for the smaller clusters is seen in all simulations in phase 2, and it is robust with respect to the cluster discrimination parameter $`dP`$. Self-organized criticality has been observed previously in lattices of pulse-coupled oscillators with diverse natural frequencies corral , where it appeared as critically distributed avalanches of simultaneous firings. The question of frequency entrainment was not addressed in that study. In our model, such avalanches are, however, impossible, since waves of firings propagate with finite speed whenever $`g<1`$.

To confirm the existence of $`g_{c1}`$ and $`g_{c2}`$, one should ideally determine their magnitude as a function of $`N`$, and see whether they converge to finite, separate values as $`N\mathrm{}`$. We have not been able to execute this scheme fully, because of the long computation times involved. Figure 3(b) shows estimations of $`g_{c1}(N)`$ and $`g_{c2}(N)`$ with accuracy $`\mathrm{\Delta }g=0.005`$, using a single realization, meaning that a single assignment of natural periods and a single initial condition (phases at time zero) is used for a given $`N`$. More realizations indicated that the spread of the critical values is similar to $`\mathrm{\Delta }g`$ for $`N=500^2`$, i.e. considerably less than $`g_{c2}(N)g_{c1}(N)`$. The shape of the curves in Fig. 3(b) support the hypothesis that the critical couplings converge to separate finite values. Since a critical coupling $`g_c=\sqrt{2/3}0.82`$ for global frequency entrainment was shown to exist in the corresponding one-dimensional lattice chain , we are strongly inclined to believe that a finite $`g_{c2}<0.82`$ exists here, since increased connectivity in general facilitates the appearance of order. Therefore the most important observation in Fig. 3(b) is that $`g_{c2}(N)g_{c1}(N)`$ does not seem to decrease as $`N`$ increases, suggesting that $`g_{c1}`$ is indeed lower that $`g_{c2}`$. To be confident that phase 2 between $`g_{c1}`$ and $`g_{c2}`$ exists, one should also check that the lattice does not converge to a frequency entrained state in this coupling interval, albeit much slower than above the presupposed larger $`g_{c2}`$. This was done for a $`300\times 300`$ lattice at $`g=0.525`$, which was simulated during $`4.4\times 10^5`$ t.u. After a transient of about $`1.5\times 10^5`$ t.u., the standard deviation of the mean period distribution in the lattice did not show any tendency to decrease.

Mean period landscapes from this simulation (in phase 2) are shown in the bottom row of Fig. 4. The cluster configuration never seems to stabilize. However, some clusters exist for a long time. For example, the white cluster in the leftmost panel exists to the end of the simulation, i.e at least for $`1.3\times 10^5`$ t.u. These features of phase 2 are consistent with the hypothesis of criticality. Then the system is expected to be self-similar in time, without a characteristic time interval, during which mean periods could be calculated with confidence. It should also lack a characteristic life-time for the clusters, so that some clusters exist during short time intervals, and a few for very long times. The upper row of Fig. 4 shows mean periods in the same lattice for $`g=0.40`$ (in phase 1), calculated at corresponding times. It is seen that the positions of the clusters remain essentially fixed. This seems always to be the case in phase 1. A related instability in phase 2 is that the mean period landscapes at a given time from two simulations with different initial conditions look very different. This would make the landscape in phase 2 sensitive to external perturbations. In contrast, in phases 1 and 3 the system seems to approach the same mean period landscape regardless the initial condition. These instabilities made it impossible to determine if, and how, the order parameter $`r`$ changes with $`g`$ in phase 2.

The mean period of an oscillator is defined as time goes to infinity. The instability of the clusters in phase 2 indicates that these numbers, if they exist, could be different from mean periods that are measured during a finite time. We can imagine three possibilities: 1) The mean periods do not exist, 2) they exist and are the same for all oscillators, and 3) they exist, and are different. If frequency clusters appear and disappear at completely random places in the network as time passes, then alternative two could be true. If not, we could have alternative three. If this alternative is correct, one could ask how the mean period landscape would look like. It is not self-evident that there would still be a macroscopic, percolating cluster. If this is the case, or alternative 1) or 2) is true, one has to assume a finite measurement time to say that $`0<r<1`$ in phase 2.

In addition to the study of mean period landscapes, we investigate the distribution $`n(P)`$ of the mean periods $`P`$. Figure 5 shows such distributions for different values of $`g`$. It is seen that the maximum shifts towards smaller values of $`P`$ as $`g`$ increases. For a one-dimensional chain, it was shown that the entrained period was always that of the fastest oscillator in the thermodynamic limit chain . The same could very well be true here, in which case the distribution would become a delta spike $`\delta (P1)`$ at $`g_{c2}`$. For $`g0.4`$, the distribution is approximately symmetric around the maximum at $`P_0`$. Assuming a functional form $`n(P)\mathrm{exp}(\alpha |PP_0|^\beta )`$, it seems that $`\beta `$ decreases through one as $`g`$ increases. Above $`g=0.45`$, the distribution becomes progressively asymmetric, with a wider and wider tail of long periods. It is seen in Fig. 5(c) that this tail obeys a power law $`n(P)(PP_0)^\gamma `$ in phase 2. Subtracting one from the slope in this cumulative plot, we get $`\gamma 2`$.

The appearance of oscillators with long mean periods parallels the appearance of completely silent oscillators, mentioned previously. To investigate the behavior of such slow or silent oscillators, time series of firing intervals for individual oscillators are studied. It is seen that oscillators with very long mean periods most often fire with normal intervals, but also experience very long periods of silence, appearing as intermittent “bursts”. Oscillators with more normal mean periods also display intermittency, but the bursts of silence are much shorter. The bursts become more scarce as the coupling increases towards $`g_{c2}`$, as usual for an intermittent signal as we approach the bifurcation point where they disappear.

In summary, we have found strong indications for the existence of two phase transitions in a large lattice of pulse-coupled oscillators with diverse natural frequencies. At the first transition, one frequency entrained cluster becomes macroscopic. At the second, all oscillators frequency entrain. Between the two transitions, the system seems critical, with self-similar cluster size distribution for the microscopic clusters, and strong fluctuations, with possible self-similarity in time. These phenomena calls for further investigations, in particular to find out if they are generic for large locally coupled oscillator networks.

We thank Martin Folkesson for making the initial simulations in this study.


format: text

text:
# Central Proper-Motion Kinematics of NGC 6752 (catalog )1footnote 11footnote 1 Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.

## 1 Introduction

With HST now having been in orbit for over 10 years, we can now combine high resolution images over a sufficiently long baseline to measure proper motions in the nearer globular clusters. Proper motions give us two further components to each star’s 6-dimensional phase space coordinate. When combined with radial velocities and positions, five-sixths of the phase space can be mapped out, leading to greater understanding of the structure and evolution of globular clusters. Radial velocities have been measured for many globular clusters, although the number where there are sufficient measurements for good statistics is still small. This is, in part, a consequence of the amount of telescope time required to acquire spectra of many stars. Fiber-fed spectrographs such as HYDRA or 2DF go some way to alleviating this problem, but fiber crowding restricts their use in the centers of globular clusters. Fabry-Perot techniques use direct imaging to measure radial velocities in more crowded regions, but their spectral resolution and range is limited, and calibration can be difficult. Proper motions, on the other hand, rely only on direct imaging, and fairly short exposures suffice. The time we are required to wait between exposures to measure a given velocity increases with the distance of the object, but decreases with the telescope’s spatial resolution.

At about 4 kpc, NGC 6752 is a good candidate for a proper-motion study. It is a prominent southern cluster, which has revealed some unusual features in recent years. It is not clear whether NGC 6752 is truly anomalous in some way, or whether the wealth of observational detail available has uncovered features that are present in many other less well-studied clusters as well.

Exposures of NGC 6752 taken with the HST pointed at the position of the cluster center according to Shawl & White (1986) clearly miss the highest density region by about 10″. Examination of Figure 1 of Rubenstein & Bailyn (1997) clearly demonstrates this. The method Shawl & White used is easily misled by bright foreground stars and this appears to be the case here. A bright star some 101″ north-west of the cluster may have biased the determination of the center of the cluster (R. White, 2002, private communication). Instead, we use the center determined by Rubenstein & Bailyn (1997), which is given only in their pixel coordinates. We have used the World Coordinate System information in the HST frame header to get J2000 equatorial coordinates of $`\alpha =19^h10^m52\stackrel{\mathrm{s}}{\mathrm{.}}16,\delta =59\mathrm{°}59\mathrm{}4\stackrel{}{\mathrm{.}}0`$ This is probably good to better than $`0\stackrel{}{\mathrm{.}}2`$.

Djorgovski & King (1986) identify NGC 6752 as a “post-core-collapse” globular cluster on the basis of its density profile, which does not fit a King model. However the profile is also poorly fit by a power law. Furthermore, Rubenstein & Bailyn (1997) have identified a high binary fraction in the core, and dynamical models suggest that a true core-collapse phase cannot be achieved until most of the binaries have been “burned” (Gao et al., 1991). Interestingly, the binary fraction appears to drop sharply beyond 11″ from the cluster center (Rubenstein & Bailyn, 1997). The drop-off is much sharper than would be expected from stars of twice the turnoff mass.

NGC 6752 is among the clusters showing a bimodality in CN abundances on the giant branch (Norris et al., 1981). It was recently found that in 47 Tuc, this bimodality exists on the main sequence itself (Cannon et al., 2002; Harbeck, Smith, & Grebel, 2003). There is some evidence that this is also the case in NGC 6752 (Gratton et al., 2001). The horizontal branch of NGC 6752 is blue as befits its metallicity (\[Fe/H\]$`=1.64`$), but its extent is unusual. This is one of the clusters where the horizontal branch features an extreme blue tail reaching stars with $`T_{eff}32,000`$ K (Momany et al., 2002). The kinematics of the cluster itself are also unusual. Dinescu et al. (1997) find it to be the oldest, most metal-poor cluster with a measured space velocity to have disk-like kinematics. In its kinematics, it is quite similar to 47 Tuc despite being much more metal poor. That an old, metal-poor cluster can have the kinematics of the thick disk has interesting implications for the history of the Galaxy.

A variety of interesting compact objects has been reported in NGC 6752 . A recent Chandra X-ray study (Pooley et al., 2002) has identified 19 sources within the half mass radius of the cluster, of which over half are apparently CVs belonging to the cluster. This contrasts with the results from 47 Tuc (Grindlay et al., 2001) in which most of the faint X-ray sources appear to be millisecond pulsars. However, radio observations have identified five millisecond pulsars in NGC 6752 (D’Amico et al., 2002). Three of these are close to the center of the cluster, and two of these three have negative period derivatives. These set lower limits on the gravitational acceleration of the pulsars by the cluster potential. D’Amico et al. set a lower limit of 10 on the mass-to-light ratio in solar units inside 7″. This, however, was based on the Shawl & White (1986) center. More recently, they have also recognized the deviation of the true center from that of Shawl & White (1986). Using a center very close to ours, the pulsars lie closer to the center and the new lower-limit on the mass-to-light ratio is 6 in solar units (A. Possenti, 2002, private communication). This is still greater than what would be suggested by the only previous velocity study of NGC 6752 in which integrated-light spectra by Dubath, Meylan, & Mayor (1997) yielded a central velocity dispersion of 4.5 km s-1. As we will see below, we obtain a significantly higher velocity dispersion, so this may not be a problem.

The outermost of the millisecond pulsars is in a $`P_{\mathrm{orb}}=0.86`$ day circular binary at least 3.3 half-mass radii (8 pc) away from the cluster center (D’Amico et al., 2002). This is an unusual place to find a binary millisecond pulsar, and the circular orbit sets firm limits on how it came to be there. Colpi, Possenti, & Gualandris (2002) argue that the most likely scenario leading to this system’s present location is one where the pulsar system is scattered by a binary black hole with total mass of 3 to 100 $`M_{\mathrm{}}`$. Such a binary black hole is probably distinct from any black hole which may form part or all of the dark component implied by the high central mass-to-light ratio.

As noted above, it is not clear whether these various anomalies mean that NGC 6752 is unusual, or simply well-studied. But either way, a more detailed understanding of its dynamics would appear to be important. The remainder of this paper is in four parts. In the next section we discuss in detail our observations and reduction procedure. Following that we present our velocity dispersion profiles. The methods for inferring such kinematic conclusions from the data are, in our opinion, unsatisfactory. We discuss this in §4. We sum up our results in the concluding section.

## 2 Observations and Reduction

Our proper motions are based on two epochs of WFPC2 data from HST. The earlier epoch consists of 116 PC-only frames taken as part of program GO 5318 on 1994 August 18-19. These were part of a search for variable stars in the center of the cluster. In this study we have used the 26.0 s F555W frames. These are arranged in a $`3\times 3`$ dither pattern centered on the Shawl & White (1986) center. We shall refer to this as the “E1994” data.

The second epoch consists of 12 WFPC2 frames taken as part of program GO 7469 on 1999 September 12. There are 4 frames at each of 3 positions arranged in an ‘L’-shaped pattern with offsets between them the size of a PC field. The original intention had been to reproduce the position of frames taken with the original PC in 1991 August as part of program GTO 2943. In these data the center of the cluster appears principally on the WF chips. This is due to the error in the assumed center, as well as a measure of confusion in defining the offsets required to put the WFPC2 PC in the same places as in the 1991 data. We shall refer to this epoch as “E1999”.

The result of the combination of these two data sets is that while the area is limited to that of one PC frame by the E1994 data, our proper motions will be limited by the measurement precision of the E1999 positions.

All the frames were reduced by the methods of Anderson & King (2000). This uses the images of stars observed in multiple frames to measure the effective point-spread-function in the frame. This ePSF takes into account the effects of the pixellation on the undersampled WFPC2 images and allows us to measure the position of point sources to a high precision. Our reductions leave us with 116 raw star lists for the 1 field covered by the E1994 data and 48 lists for the over-lapping fields of the E1999 data. What is now required is to correct the positions in these lists for two things: 34-row error and distortion in the optics. The 34-row error (Anderson & King, 1999) was corrected using the model described in the Appendix. Distortion in each WFPC2 chip was corrected using the results of Anderson & King (2003).

For the E1994 data, the corrected positions were transferred into a reference frame by polynomial transformations including terms up to second order in the positions. The quadratic terms were required to allow for residual variation in the distortion terms. The reference frame consisted of the weighted means of the translated positions. These were stable after two iterations. An iterative $`3\sigma `$ clipping algorithm was employed in producing the mean positions. Any observation that deviated by more than $`3\sigma `$ in either dimension was rejected in producing the mean position.

The combination of the E1999 data was complicated by the multiplicity of individual pointings of the various frames. What we ultimately did was to first rigidly transform the corrected positions for each frame to a common one. The data from each CCD was transformed independently. The mean positions in this common frame formed the initial approximation to a reference frame. The rigidly transformed lists were then iteratively transformed onto the reference frame and then combined as discussed in the previous paragraph to form a new reference frame. Several iterations were required before the difference in positions between two successive approximations were reduced to insignificance. Both quadratic and cubic polynomials were used for the final transformations, but the final positions agreed within the estimated errors. We also tried local transformations using the 100 nearest neighbors to each star, but, again, there was no improvement in the quality of the final positions.

The final positions in each epoch were calculated as weighted means of the translated positions in each frame. The uncertainty in the positions were estimated from the standard deviation of the individual positions after the rejection of $`3\sigma `$ outliers.

To measure dispersions, it is important to have confidence in our estimate of the errors in the individual positions. In order to validate our error estimate, we divided the E1999 data into two halves and calculated the positions and errors for each half separately. (Our error budget for the proper motions is dominated by the uncertainty in the E1999 positions. Since the E1994 data consists of 116 PC frames, the uncertainties there are an order-of-magnitude smaller.) For each star we calculated the position differences in each dimension divided by the standard deviation in the difference. If the errors are correct and normal, these should be distributed as a Gaussian with zero mean and unit standard deviation. What we found was that while the mean was consistent with zero, the standard deviation was 0.90, indicating the errors to be underestimated. All errors in the final E1999 data were then divided by 0.90 to correct for this. As will be seen, the measured velocities are large enough that this 10% effect in the error estimation will have no appreciable effect on the final velocity dispersions.

The final E1994 and E1999 star lists were then matched and a global quadratic transformation was calculated between them using the E1994 positions as the reference. Cubic and local transformations give the same results with larger errors for the local transformations. The differences between the transformed positions were then taken as the proper motions in the E1994 coordinate system. A plate scale of 0.04558″ pixel-1 appropriate to the PC1 chip (Biretta et al., 2001) and the epoch interval of 5.07 years were used to convert the proper motion into units of mas yr-1. The proper motions in each dimension were examined for systematic effects with respect to the stellar coordinates in both epochs. None were seen. We also calculated two proper motions for each star by dividing the frames at each epoch into two samples. The proper motions in each half-sample were consistent within the estimated errors. Thus, we are confident that these have been estimated correctly.

In total there are 19888 stars with four or more observations in the E1999 data. We require this number of observations in order to be able to make reliable error estimates. Our experience when dividing the sample into halves suggests that for these data four observations should be sufficient to estimate the errors adequately. 1281 of these stars match stars included in the E1994 sample. For these stars the median combined error of the two proper motion components is 0.31 mas yr-1, with a mode of 0.17 mas yr-1. The distribution of uncertainties drops off rapidly, so we have established a cut-off at 1 mas yr-1 and retained all 1121 stars with errors less than this. It should be noted that the stars rejected are not exclusively those with large proper motions. The proper motions of the retained stars are shown in the upper-right of Figure 1. Note that no field population is apparent in this diagram. This is not unexpected given the relatively high galactic latitude of NGC 6752. There is a similar lack of non-cluster stars in the color-magnitude diagram of Rubenstein & Bailyn (1997).

Also in Figure 1 we show histograms of the $`x`$ and $`y`$ components of the proper motions in E1994 coordinates. Analysis of the distributions show that they are not described well by single Gaussians. In each component there is a central core with the addition of a second, more broadly distributed component. We interpret the latter as likely due to blended stellar images in the E1999 WF data being matched with single stars in the E1994 data and giving spuriously high proper motions. We fit a two-Gaussian mixture model (Nemec & Nemec, 1993) to the proper motions in each component. Then, for each star, we calculated the probability that the star was in each component and a joint probability that the star was in the narrow component. The mixture models had 12 to 18% of the stars in the broad component and, in the joint distribution, 15% had probabilities of being in the narrow component of less than 50%.

As support for this interpretation we show in Figure 2 the differences between the measured $`V`$ magnitudes in the E1999 and E1994 data sets. For the latter we have used the photometry of Rubenstein & Bailyn (1997). The stars with joint probabilities of less than 50% are shown as solid symbols. These are found preferentially on the bright side of the distribution. This is as we would expect if the images in the E1999 data were the combination of two stellar images. The low-probability stars are concentrated to the center as we would expect. A KS test gives a probability of zero that the two groups come from the same radial distribution. Even in the most crowded regions the majority of the stars are well measured, however. Only 19% of the stars within 11″of the center have low probabilities. This rises to 26% within 5″and 37% for the 35 stars within 3″of the center.

Rejecting those stars with joint probabilities less than 50%, we have 950 stars left. Our final step was to resolve the proper motions into radial and tangential components with respect to the radius vector between each star and the Rubenstein & Bailyn (1997) center. The proper motions as a function of radius and azimuth are presented in Figures 3 and 4. The azimuthal sample we use below is a 431 star subsample of the full star list. They have been selected to cover a uniform radial region so as not to introduce any biases arising from any radial variations in the velocity dispersion. The stars selected cover $`160\mathrm{°}`$ about the cluster center to $`16\stackrel{}{\mathrm{.}}7`$. The E1994 frames lie within 4°of North along the $`x`$ and East along the $`y`$ axes, so the azimuths in Figure 4 run roughly from North through East.

## 3 Velocity Dispersion Profile

Extracting kinematic information from a set of velocities is as traditional as it is fraught with difficulty. It is customary to interpret the observed velocities as coming from a Gaussian distribution and use them to measure the dispersion in the distribution, usually as a function of position. There are certainly reasons to believe that this is not actually the case, but, in the absence of sufficient data to attempt to recover the distribution function non-parametrically (as in, for example Merritt, Meylan, & Mayor 1997), or the use of models transferred into the observational domain, it suffices as a first cut at interpreting the velocity data. It is in this spirit we calculate the velocity dispersions for NGC 6752 from our proper motions.

In order to investigate the radial and azimuthal dependence of the velocity dispersion we have binned the stars into equal number bins and used a maximum-likelihood technique to estimate the dispersion in each bin. For a given star with observed velocity $`v_i`$ and uncertainty $`ϵ_i`$, if the systemic velocity is $`\overline{v}`$ and the velocity dispersion in its bin is $`\sigma _r`$, for bins $`r=1,\mathrm{},N_b`$, then the likelihood of measuring $`v_i`$ is given by

$$L(v_i|\overline{v},\sigma _r,ϵ_i)=\frac{1}{\sqrt{2\pi \left(ϵ_i^2+\sigma _r^2\right)}}\mathrm{exp}\left[\frac{\left(v_i\overline{v}\right)^2}{2\left(ϵ_i^2+\sigma _r^2\right)}\right].$$
(1)
For each bin, the total likelihood is the product of eq. 1 for all stars in the bin. From the data $`\overline{v}`$ is consistent with zero and we have assumed this value in what follows. For the radial sample we have used 5 bins of 190 stars each. For the somewhat smaller azimuthal sample we have used 3 bins of 108 stars and one of 107. The $`1\sigma `$ uncertainty in the dispersion estimate has been taken as the half-width of the symmetric region about the maximum likelihood containing 68% of the probability. For the binning used, this uncertainty in the measured dispersion is approximately 4% for the radial bins and 6% for the azimuthal ones.

Figure 5 shows the velocity dispersion profile in the radial direction for the two components. The proper motions have been converted into velocities assuming a distance of 4.17 kpc. This corresponds to a distance modulus $`(mM)_0=13.1\pm 0.1`$, and represents a compromise between Renzini et al. (1996) who found a white dwarf distance modulus of $`13.05\pm 0.1`$ and Carreta et al. (2000) who find $`(mM)_0=13.21\pm 0.04`$ using Hipparcos subdwarf distances. At this distance 1 mas/yr = 19.77 km s-1.

Inside 11″, i.e. the central, binary-dominated core, the dispersions in the two components are consistent and flat with a dispersion of about 12.4 $`\pm 0.5`$ km s-1. This is over twice the dispersion in the radial velocity as measured by Dubath, Meylan, & Mayor (1997) and is inconsistent with their value. If we adopt the surface brightness profile parameters of Trager, King, & Djorgovski (1995), then our velocity dispersion suggests a central $`V`$-band mass-to-light ratio of 9 in solar units. This is consistent with the lower limit from the pulsar timings, but much larger than that typical for galactic globular clusters (McLaughlin, 2000).

Outside the core, the dispersion drops with radius and it appears that the dispersion in the tangential component drops more quickly than does the dispersion in the radial component. The presence of such radial anisotropy so close to the core would be quite unexpected. In the sophisticated N-body models of Baumgardt & Makino (2002), for example, little anisotropy at all is seen in the inner half of the models, and what there is is tangentially anisotropic, not radially as is the case here. The differences between the radial and tangential dispersions are small, at most $`2\sigma `$, but there is a systematic trend that is highly suggestive.

Figure 6 shows the velocity dispersion as a function of azimuth about the center of the cluster. As discussed above, this does not cover the entire $`360\mathrm{°}`$ due to the off-center pointing of HST. Any differences between the radial and tangential velocity dispersions do not appear significant nor systematic.

## 4 Discussion

Our main results are that the central velocity dispersion in NGC 6752 is larger than that measured previously on the basis of radial velocities, and that there appears to be substantial radial anisotropy just outside the core of the cluster. The significance of these results depend on how much credence we give to the velocity dispersion estimates.

Our analysis of the velocities has been fairly conventional, in that we have binned the data to derive a velocity dispersion at a number of distances from the cluster center. The only unusual feature is the use of variable uncertainties, but, in the limit where the errors on the velocities are equal, our method reduces to the quadrature removal of that error from the standard deviation of the measurements. Our advantage over Pryor & Meylan (1993) is that our repeated measurements provide a good estimate of uncertainties on a star-by-star basis. The common use, however, of this technique to explore globular cluster dynamics has a number of serious conceptual problems. Here we list and discuss the difficulties in the standard approach to exploring radial trends in velocity dispersion.

1) Binning is not, in general, an appropriate way of characterizing radial change. There are well-known biases involved in binning, including unequal distribution of data over the bin, changes in expected values of a quantity over the bin, the intrinsic trade-off between resolution (the number of bins) and precision (the number of data points within a given bin). Bayesian techniques (e.g. Drukier et al. 1998) can be used to get around this problem, but care must be taken since one is interested in the radial change of a distribution, rather than of a quantity.

2) In a similar vein, the mass dependence of the velocity distribution must be taken into account. Mass segregation is a clear feature of globular cluster dynamics, and therefore the velocity dispersion depends on the mass of the objects being measured. This mass varies from object to object, and the distribution of masses varies with radius. Even objects with very similar observational characteristics may have quite different masses. Consider, for example, the case of high-mass-ratio, double-main-sequence, binary stars, which may lie almost on top of the isochrone for single stars. Since the fraction of binaries doubtless varies strongly with radius, as in the present case, one may not be measuring the same mass distribution at different radii even when one can control for the luminosity distribution.

3) The measurement of the velocity dispersion in a globular cluster is hampered by the limited number of stars available. This is especially true in the center of a cluster, where the measurement of the velocity distribution may be forever limited by the finite number of stars present. Examination of any of the HST images of the central regions of globular clusters will quickly illustrate this point. Our proper motions already cover nearly every star within 11$`\mathrm{}`$ of the center down to well below the main sequence turn off. Since the mass function in the core of this cluster (like most dense clusters) is flat, or perhaps even inverted (Rubenstein & Bailyn, 1999), we cannot gain more than a factor of two by going fainter on the main sequence. There is a small vignetted section in the E1994 data, some small regions were insufficiently observed in the E1999 data, and some stars were lost due to blending, but the potential gain in total number of stars from those lost in this way is quite small. Unless we can detect large numbers of white dwarfs and neutron stars, the total number of stars in the core of any given cluster is limited.

4) The shape of the velocity distribution is not in general understood. One of the most important assumptions in both binned and Bayesian methods of exploring the velocity distribution in clusters is that this distribution can be characterized by a single number (the dispersion) generally taken to be the standard deviation of the distribution. Unless the distribution is Gaussian, a dispersion measured this way is not the same as the second moment of the true velocity distribution, which could bias further analysis. But this assumption of normally distributed velocities must be wrong at some level. Even in single-mass isotropic King models the velocity distribution is not a Maxwell-Boltzmann distribution, but one with a cutoff. In a realistic treatment, it is far from clear what a measurement of a “dispersion” really means. It is particularly problematic in small samples (in one bin, for example), where a single outlier can greatly change the measured value of any moment of a distribution.

5) Finally, the very concept of a velocity distribution is poorly defined in a globular cluster. Defining a “distribution” imposes an assumption of continuity on the population. In a globular cluster, and especially in the center, the dynamics may be severely influenced by non-fluid effects. There are only a finite number of stars, and the stochastic properties of their particular location in phase space may be of considerable importance. This is particularly true given that the gravitational potential, which gives rise to the velocity distribution, is created by those very same stars. The situation is quite different in the case of the center of a galaxy. In that case, the black hole is so much more massive than the stars that the potential is smooth and well-defined, and there are generally enough stars or gas particles that better observations lead to more independent measurements in an obvious way. Even in the case of a galaxy, however, it is not clear that radial velocity measurements alone are sufficient to independently constrain both the mass-to-light ratio and the mass of any central black hole (Valluri, Merritt, & Emsellem, 2002).

There has been in the past some level of recognition of some of these problems. None of the solutions to date, however, have been entirely satisfactory. Alternatives to standard binning have been advocated by Gebhardt et al. (1994) and Gebhardt et al. (2000) in their LOWESS and ROTATE2 algorithms. These algorithms attempt to provide a non-parametric characterization of the velocity-dispersion profile. While the first explicitly makes a smoothed, non-parametric estimate of the velocity-dispersion profile, the latter amounts to nothing more than a maximum-likelihood method applied to a set of overlapping bins of undisclosed and variable width and population. The derived profiles are frequently characterized by large, high-frequency oscillations. The character and location of these oscillations are easy to replicate using a more standard binning procedure. As the uncertainties are taken from bootstrapped samples based on the inferred profile, the error bands parallel the estimated profile. It is then hard to decide how much credence to give the apparently significant fluctuations in the dispersion.

Drukier et al. (1998) attempted to use a Bayesian approach to measuring the dispersion profile both in bins and in various parametric forms. The underlying method is similar to the maximum-likelihood method we use here, but their method of hypothesis comparison is incorrect. What that attempt at hypothesis comparison does highlight is the wide range of possible parametrizations and the necessity to objectively choose the best amongst them. To take some concrete examples, we might find of use answers to the following sorts of questions: Are our velocities here better represented by a velocity-dispersion profile consisting of a single power-law or one with a core? At what level do the radial and tangential components have profiles that disagree? Ultimately, however, we would wish to get beyond parametrization to allow the data to “speak for themselves”.

Merritt (1993, 1996) and Merritt & Tremblay (1994) have advocated a series of non-parametric inversion techniques to move from the observed surface-density profile and projected velocity-dispersion profile to the more intrinsic density profile, velocity-dispersion profile, and ultimately the potential and distribution function. These methods have much to recommend them and a practical application has been demonstrated in Merritt, Meylan, & Mayor (1997) in the case of $`\omega `$ Cen. Merritt’s techniques find the minimum of a penalized log likelihood of the form

$$\mathrm{log}_p=\underset{\mathrm{data}}{}\left[\frac{g_p𝐀g}{ϵ}\right]^2+\alpha P\left(g\right),$$
(2)
where $`𝐀`$ is a projection operator that brings the desired function, $`g`$, into the observable space where we observe $`g_p`$. The quantity $`ϵ`$ is the measurement error and $`P`$ is a functional that assigns large penalties to noisy solutions for $`g`$. The degree of smoothing depends on the value of $`\alpha `$ and is necessary since the inversion problem is under-constrained. The function $`g`$ can be, for example, the mass-density profile, a velocity or velocity dispersion, or the distribution function. Until one gets to the inverted distribution function one is still relying on the assumption that the velocity distribution is indeed Gaussian in estimating its dispersion and from this the potential. These techniques give results that are generally independent of any functional form imposed by parametrizing assumptions. The resulting functions are those that best match the data, given the degree of smoothing. They do require large numbers of velocities and cannot be applied if a more limited quantity is all that is available. What is lacking is a measure of how convincing a representation that is.

Gebhardt & Fischer (1995) use a philosophically similar approach, although quite differently implemented, to estimate the mass-density profile from the observed surface-density profile and velocity-dispersion profile estimated by the LOWESS algorithm (Gebhardt et al., 1994). They use spline-smoothed estimates of the profiles which are then numerically differentiated and integrated to obtain, via Jeans equation, the mass-density profile. Their method, like Merritt’s, recovers the original model for simulated data. The assumption, though, is that the velocity ellipsoid is isotropic.

A somewhat more direct approach to measuring the potential of a globular cluster is seen in Gerssen et al. (2002, 2003) for M15. They invert the surface-density profile under various assumptions for the mass-to-light ratio to derive, via Jeans’s equation, velocity-dispersion profiles. They then compare these to the observed velocities using a maximum-likelihood technique. The large discrepancies between these profiles and that measured using the Gebhardt et al. (2000) method, despite the supposedly consistently high likelihood, cannot help but leave us feeling that one or both of the dispersion estimator and the likelihood argument are incorrect.

Indeed, this work highlights the risk of over-interpreting velocity data. Gerssen et al. claim to have found evidence for a $`3,000M_{\mathrm{}}`$ black hole in the center of M15. This claim is largely independent of the strength of their claimed fit to the M15 velocity-dispersion profile discussed above. Instead, the evidence in support of a black hole rests on the velocities of about 12 stars within 1$`\mathrm{}`$ of the cluster center, 8 of these being new HST/STIS measurements with median error of 5.7 km s-1. Even in their original presentation (Gerssen et al., 2002), the evidence for a black hole is, at best, a $`1\sigma `$ detection. In the final result Gerssen et al., 2003 correcting for a scale error in Dull et al. (1997) (Dull et al., 2003), even the claim for a black hole becomes ambivalent and dependent on the apparent existence of a 20,000$`M_{\mathrm{}}`$ black hole in the M31 cluster G1(Gebhardt, Rich, & Ho, 2002). But the claim for G1 is even weaker. The method used for estimating the mass of the putative black hole there is known to give false minima in $`\chi ^2`$ in the face of insufficient data even for galaxies (Valluri, Merritt, & Emsellem, 2002). A globular cluster, with its problems of granularity, stellar sparseness, and small radius of influence for a black hole, is probably the worst place to apply these techniques.

We do not suggest that we know what is the correct way to proceed here. What we would like to do is raise the level of awareness of the issues involved. What is needed is to identify techniques which work optimally for a particular quantity and type of data in the face of the limits discussed above. Application of these methods must be tempered by an appreciation of their biases and limitations, and with proper attention to whether the conclusions drawn are actually supported by the data or are just artefacts of the method used.

## 5 Conclusions

Our HST proper motions suggest that the velocity dispersion in the center of NGC 6752 is surprisingly large. At 12.5 km s-1 it is much larger than the measured dispersion along the line of sight. While there is some uncertainty in the distance to NGC 6752 it is certainly known to better than the factor of roughly two which would be required to bring the two measurements into agreement. Our dispersion indicates a V-band mass-to-light ratio of 9. Radial anisotropy appears outside of about 1.5 core radii. Further observations are required to confirm the measurements presented here and to extend them to a larger area. More importantly, the theoretical basis underlying the kinematic interpretation of stellar velocities in globular cluster needs to be reconsidered.

###### Acknowledgements.

We would like to thank J. Anderson for the use of his astrometry code and him, I.R. King, and S. Shaklan for useful discussions. The comments of the anonymous referee and Anderson on the original version of the appendix prompted us to re-examine our original 34-row correction. The result is the much improved version presented below. Drukier and Bailyn were supported by a NASA LTSA grant NAG5-6404. van Altena and Girard acknowledge the long-term assistance of their colleagues on the HST Astrometry Science Team (W. H. Jefferys, P.I. and G. F. Benedict, Co-P.I.) who collaborated in the GTO observations and made it possible for us to obtain the GO observations that provided the second-epoch observations used in this investigation. They also acknowledge financial support from NASA for the GTO observations and of the STScI for the GO observations.

## Appendix A 34-row correction

The so called 34-row error (Anderson & King, 1999) originated in the construction of the WFPC2 CCDs. It manifests itself as a periodic effect where every 34th-row of the CCD is a few percent narrower than the others. In this case we lose uniformity in the $`y`$-coordinate. Every 34th-row occupies less physical space, but is still assigned the same unit distance in coordinate space. The net effect is that the observed $`y`$ positions of stars deviate systematically from their true positions in a series of step functions.

Anderson & King (1999) correct for this using a sawtooth function with period 34.1333 pixels based on measured residuals to a mean position. This combines a step function with a scale stretch. Anderson & King argue that this is sensible since, if the 34th row is narrower than the average row, then the other 33 rows are wider. But this is to conflate the plate scale with the physical source of the defect. And, while the correction suggested by Anderson & King (1999) is correct on average, it strikes us as less than satisfactory. The inter-flaw spacing, we feel, ought to be an integral number of pixels, presumably 34. Further, the flawed rows should be identifiable, allowing for a correction which would be correct in detail, rather than just on average. To that end, we have reconsidered the nature of the 34-row flaw and have reexamined the WFPC2 flat fields to test our reconstruction.

Shaklan et al. (1995) describe the construction of the reticles for the WFPC2 CCDs in some detail. A 0.5 µm by 0.5 µm e-beam was scanned in the $`y`$-direction for 512 µm at a time laying down the desired pattern. The beam was shifted by 0.5 µm in the $`x`$-direction and the scan repeated until the pattern for the full 800-pixel width of the CCD had been laid down on the reticle. The e-beam was then shifted by 512 µm in $`y`$ and the next 512 µm region was scanned. This was repeated until the full 800-pixel height of the CCD had been laid down, completing the pattern. What appears to have happened is that the 512 µm $`y`$-scans actually occupied 512.5 µmon the reticle and, consequently, the CCDs. The shift between vertical scans remained 512 µm so at the edge of each vertical scan a 0.5 µm region was scanned twice by the e-beam with two consecutive regions of the CCD pattern. We postulate that the e-beam was magnified by 0.1% in the $`y`$-direction resulting in slightly rectangular pixels, 15 µm wide as planned, but approximately 0.1% taller. A slight inclination of the reticle plane with respect to the e-beam could account for this.

The result is that 512 µm of CCD pattern were replicated onto 512.5 µm of the reticle. At the edge of each scan we find 1 µm of pattern projects onto 0.5 µm of the physical CCD and mapping from pattern space to physical position on the CCD is double valued. In that 0.5 µm region the reverse mapping goes to two parts of the pattern spaced by $`0.5/\delta `$ µm where $`\delta 1.001`$ is the magnification factor. This puts a 512 µm periodicity into the underlying pattern.

In addition, there is a second periodicity, one of 15 µm, corresponding to the pixel size. This maps into $`15\delta `$ µm on the CCD for the normal pixels. Every approximately 34 pixels you get the pixel that contains the doubled-up flaw and, consequently, this row is narrower by about 0.5 µm.

The details of which rows are damaged depends on the details of the relative phasings of the 15 and 512 µm periods and how these map onto the 800 pixels that define the $`y`$-dimension of the CCD. We can certainly see that the 15 µm cycle has phase zero at the edges of the CCD. It appears, however, that the 512 µm cycle does not start at the same place. If you make some assumptions about these phases, you can easily show that the typical inter-flaw spacing is 34 rows. There is a 15-flaw pattern where, in the 7th and 15th flawed pixels, the spacing is 35 pixels. There are two cases, where the 34th pixel is either normal or where both the 34th and 35th pixels are damaged by a doubled-up region that contains a pixel boundary. Thus, on average, as noted by Anderson & King (1999), the periodicity is $`512/15=34.1\dot{3}`$ pixels, but this comes from a mode of 34 pixels in 13/15 cases and 35 in the other 2/15. The correction should take into account the detailed position of the flawed rows rather than only dealing with the problem on average.

In order to test our reconstruction, we examined the WFPC2 flat fields. To do this, we took each flat and divided it into a $`1\times 25`$ pixel median filter of itself. We averaged together each row and all four WFPC2 CCDs since the flaw pattern is identical on each. Figure 7 shows the F343N, F502N, and F953N flats after they have been processed in this manner and the vignetted region trimmed off. The flawed rows are readily identifiable by their low intensity. In this figure, the vertical axis is the relative response function of the rows due to the 34-row anomaly. The flawed rows are listed in Table 1. The flaw number is given by $`(15r/512)`$, where $`r`$ is the flawed row. As expected, the modal inter-flaw spacing is 34 rows. As Anderson & King (1999) reported in proof from a private communication from R. Gilliland, there are a number of flaws which occupy two adjacent rows. Our model predicts that this would occur when the 1µm doubling region contains a pixel boundary, but this should only happen every 15th flaw, i.e. twice across the CCD, not the six times seen. We note that the flaw pairs do obey our expected 15 flaw periodicity as does the relative response at all wavelengths.

We identify flaws 6 and 21 (rows 235-6 and 747-8) with the expected situations where the flaw region contains a pixel boundary. The two rows in each flaw appear to have the same size, roughly half the usual size. Flaws 13 and 14, half-way between as expected, also have a doubled structure, but, in these cases one row is clearly more flawed than the other. We refer to the row with the deeper signal as primary and the other as secondary. The primary flaws are in rows 474 and 509, with the expected 35 row spacing. Other flaws in the vicinity of these two types of double flaws also show the primary/secondary structure, particularly flaws 7 and 22. Without understanding the details of the flawed regions on the level of the CCD electronics, it is impossible to fully explain this pattern, but it appears to us that the flaw impacts a region around it to a range of several µm. If a pixel boundary falls within the affected region, then the neighboring row will also appear slightly flawed in the flat fields as seen in Figure 7. The variation in the relative response with wavelength we also attribute to the electronic structure of the CCD pixels. Longer-wavelength light penetrates deeper into the CCD, so the variations with filter are probably telling us something about the physical properties of the flawed silicon. Note, however, that the pipeline flat-field frames used here are combined from flat fields in several filters (Koekemoer, Biretta, & Mack, 2002). The weakness, or indeed absence, of flaws 2, 10, and 17 at the bluer wavelengths where they are strongest in the infrared, is particularly interesting. Our model predicts that these flaws (note the 15-flaw period again) are the most likely to have the doubled region lying in the center of the pixel. All this confirms the report by Gilliland in the Anderson & King (1999) note. As noted there, these variations will affect the photometric 34-row correction, but should not apply to the astrometric correction.

As Shaklan et al. (1995) note, the CCDs, to a high precision, are square. This is the result of the excess height of the normal pixels ($`15\delta \mathrm{\mu m}`$) being compensated for by the narrower, flawed, rows. The 800 pixels of the CCD are defined by just over 23 repetitions of the 512 µm period, resulting in 23 flawed rows. The 777 normal rows each exceed their expected size by $`15(\delta 1)`$, so each flawed row is $`1777(\delta 1)/(23\delta )`$ tall in units of the $`y`$-dimension of the pixels. This is 0.966 pixels for $`\delta =1.001`$ or 0.967 pixels for $`\delta =512.5/512`$. We adopt 0.967 for our correction. The difference in $`x`$ and $`y`$ scales will be taken care of in the plate solution.

Our correction is as follows. For each star with centroid not lying in a flawed row as listed in Table 1 we reduce the $`y`$-coordinate by 0.033 pixels per flaw with lower $`y`$ value. For stars with centroids lying in a flawed row, the $`y`$-coordinate is first rescaled in that pixel to occupy a pixel of size 0.967 rather than 1.000 before the offset is applied for the previous flawed rows. For the double-flaws 6 and 21 each row is treated has having half the flaw size. For the remaining flaws with primary/secondary structure, only the primary row is considered as being short. A comparison with the correction of Anderson & King (1999) shows that, apart from the scale difference, the maximum difference in the corrected positions is 0.015 pixels in the flawed rows. In keeping with our philosophy, the 0.1% scale difference between $`x`$ and $`y`$ is taken care of in the plate solution.

The data in this paper were corrected with an earlier version of this correction where a strict 34-pixel spacing was assumed and the differences in detected positions of stars in offset frames were used to infer the size of the narrow rows and the phase of the 34-pixel period with respect to the CCD. The result was a slightly larger narrow row and slightly different flawed rows at the top and bottom of the CCD. We have compared the final proper motions with the new correction and our results are unchanged.


text:
# On hydrodynamic shear turbulence in Keplerian disks: via transient growth to bypass transition

## 1 Introduction

The emergence of X-ray astronomy exposed the importance of accretion phenomena in many astrophysical systems, such as binary stars, quasars and active galactic nuclei. The possible mechanisms governing accretion have been schematically understood during the very first years of investigations (Shakura & Sunyaev 1973; Pringle 1981): the inward transport of matter and the outward transport of angular momentum in accretion disks was ascribed to a turbulent (anomalous) viscosity. The consequent development in understanding the physics of this phenomenon and the causes of that turbulence has been irregular and has taken considerable time. But substantial progress has been achieved in the nineties with the discovery of a linear instability in magnetized disks (Balbus & Hawley 1991, 1992, 1998; Hawley & Balbus 1991, 1992; Hawley et al. 1995; Stone et al. 1996). In contrast, the solution of the accretion problem in the non-magnetized case has not yet reached sufficient maturity. Moreover, the very occurrence of turbulence in non-magnetized disks has been questioned by Balbus et al. (1996) and Balbus & Hawley (1998). The reason for this situation is that cylindrical flows with Keplerian profile belong to the class of smooth shear flows, i.e. which present no inflection point; it is well known that these flows are spectrally stable, although they may become turbulent in the laboratory. The explanation of this behavior remains one of the fundamental problems of fluid mechanics. But the situation has changed in recent years, with the recognition that the so-called nonmodal approach provides an adequate, and even optimal mathematical formalism for describing the dynamics of perturbations in plane-parallel shear flows. This led to the emergence of the bypass concept for the onset of turbulence in spectrally stable flows. This concept has triggered much interest among fluid dynamicists, but it had little impact so far on the astrophysical community. That is why feel it necessary to present the salient results of that novel approach in a journal which is read by the astrophysicists, to demonstrate its relevance for the dynamics of Keplerian disks, and to join a bibliography of the subject.

Our approach parallels that of Longaretti (2002), who analyzed the existing experimental and numerical results, and who also insisted, as we shall do here, on the similarity of subcritical shear turbulence in plane and rotating flows.

Our paper is organized as follows. In Sec. 2 we outline the fundamental problems of fluid mechanics bearing on the description of plane-parallel smooth shear flows. We present the recent developments in shear flow analysis made by the hydrodynamic community, and the novel bypass concept for the onset of turbulence in shear flows. In Sec. 3 we analyze the quadratic forms of the dynamical equations and discuss the essential role of pressure perturbations, which counteract the effect of the Coriolis force in 2D flows. In Sec. 4 we describe the transient growth of two dimensional perturbations. This will allow us to stress the similarity between the behavior of Keplerian and Cartesian flows, and the negligible role of the Coriolis force in the 2D perturbation kinematics and energetics. Finally we conclude in Sec. 5 that the bypass concept provides a plausible scenario leading to turbulence in astrophysical disks. In the Appendix we present a simple sketch of the bypass scenario.

## 2 From spectral decay to transient growth and bypass transition

### 2.1 On the onset of turbulence in spectrally stable shear flows

From the theoretical viewpoint there are flows that are spectrally stable at all Reynolds numbers (e.g. plane Couette or pipe Hagen-Poiseuille flows), while some others become spectrally unstable at high enough Reynolds numbers (e.g. plane Poiseuille or Blasius flows). In the latter case the flow is characterized by the critical Reynolds number that is the marginal value of this parameter above which the spectral instability occurs. In this sense the critical Reynolds number of the spectrally stable flow is infinity.

Here we shall consider smooth shear flows (without inflection point in the velocity profile), like in Keplerian disks. These are linearly stable according to classical fluid mechanics. As shown by Rayleigh (1880), the existence of an inflection point (more precisely of a vorticity extremum, as pointed out by Fjørtøft 1950) in the equilibrium velocity profile is a necessary condition for the occurrence of a linear (spectral) instability in hydrodynamic flows. Thus smooth shear flows (flows without a vorticity extremum) are “relaxed” in this context: they are spectrally stable, meaning that exponentially growing solutions are absent. However, it is well known from laboratory experiments and from numerical simulations that finite amplitude perturbations may cause a transition from laminar to turbulent state at moderate, less then critical Reynolds number.

This has led to the development of the concept of nonlinear instability in hydrodynamics (cf. Bayly 1986; Bayly et al. 1988; Herbert 1988; Orszag & Kells 1980; Orszag & Patera 1980, 1983). Until about ten years ago, the predominant view of this laminar-turbulent transition was centered around the slow linear amplification of exponentially growing perturbations (the familiar T-S waves), which modify the flow profile and thereby allow a secondary instability, further nonlinearity and finally a breakdown to turbulent flow. According to this concept of nonlinear instability, the perturbations (and the turbulent state itself) are energetically sustained by nonlinear processes. This concept of nonlinear instability has been borrowed by astrophysicists, who still use it to explain turbulent processes in smooth astrophysical flows, where no spectrally unstable solution is known, and in particular in Keplerian disks flows (see Balbus & Hawley 1998). However, there are subcritical transition phenomena that can not be attributed to the nonlinear instability.

During the last decade of the 20th century, another viewpoint emerged in the hydrodynamic community on understanding the onset of turbulence in spectrally stable shear flows, labeled as bypass transition (cf. Boberg & Brosa 1988; Butler & Farrell 1992; Reddy & Henningson 1993; Trefethen et al. 1993; Morkovin 1993; Gebhardt & Grossmann 1994; Henningson & Reddy 1994; Baggett et al. 1995; Waleffe 1997; Grossmann 2000; Reshotko 2001; Chagelishvili et al. 2002; Chapman 2002; Rempfer 2003). Although the bypass transition scenario involves nonlinear interactions – which intervene once the perturbations have reached finite amplitude – the dominant mechanism leading to these large amplitudes appears to be linear. This concept of the onset of turbulence is based on the linear transient growth of vortex mode (aperiodic) perturbations. The potential for the transient growth has been recognized for more then a century (see Kelvin 1887, Orr 1907a,b). However only recently has the importance of the phenomenon been better understood (cf. Moffatt 1967; Marcus & Press 1977; Gustavsson & Hultgren 1980; Craik & Criminale 1986; Farrell & Ioannou 1993; Reddy & Henningson 1993; Chagelishvili et al. 1997).

The bypass concept implies that the perturbation energy extracted from the basic flow by linear transient mechanisms causes the increase of the total perturbation energy during the transition process. The nonlinear terms are conservative and only redistribute the energy produced by the linear mechanisms. (A simple sketch of the bypass scenario in a Keplerian flow is given in the Appendix.)

Thus, according to this concept, the transient growth of perturbations (i.e. the linear process) is the key element in the transition to turbulence in spectrally stable flows. The importance of this linear process has been stressed in the titles of several seminal papers: Henningson & Reddy (1994) (“On the role of linear mechanisms in transition to turbulence”); Baggett et al. (1995) (“A mostly linear model of transition to turbulence”) and Reshotko (2001) (“Transient growth: A factor in bypass transition”).

### 2.2 On the mechanism of perturbation amplification in spectrally stable shear flows

The scheme of investigation of shear flow dynamics implies the following steps: introduction of perturbations into a mean flow, linearization of the governing equations, and description of the dynamics of perturbations and flow using the solutions of the initial value problem (i.e. following temporal balances in the flow). In principle this can be done, but in practice it is a formidable task. Therefore the mathematical approach was changed. This was done by assuming that the solution is separable in eigenmodes, and then establishing the existence of at least one unstable eigensolution. This approach became canonical in time, and resulted in a shift of ground: attention was directed to the asymptotic stability of the flow and no attention was paid to any particular initial value or to the finite time period of the dynamics. Indeed, this phase of the evolution was not thought to have any significance – it was left to speculation. But recently the early transient period for the perturbations has been shown to reveal “rich” and complicated behavior leading to different consequences.

It was found in 1990s that smooth shear flows are crowded by intense processes of mean flow energy extraction by perturbations, energy exchange between perturbations, etc. even in the linear approximation, while following the classical theory they are spectrally stable and consequently “relaxed.” Especially it has been shown that a superposition of decaying normal modes may grow initially, but will eventually decay as time goes on – a new linear transient channel of energy exchange between the mean flow and perturbations has appeared. Moreover, it has been shown that transient growth can be significant even for subcritical values of the Reynolds number and that its interplay with nonlinear processes can result in transition to turbulence without any “nonlinear instability” of the flows.

### 2.3 On the nonorthogonality of linear operators and the nonmodal approach

In fact an exact resonance is not necessary to obtain the transient growth of perturbation. This is the consequence of the non-normal character of operators that describe the linear dynamics of perturbations in flows (see, e. g. Reddy et al. 1993). The fact that the eigenfunctions of the linearized Navier-Stokes equations are not orthogonal (i.e. the operator is non-normal) is enough to allow for solutions that exhibit transient growth, depending on the initial conditions, before finally decaying (see Criminale & Drazin 1990; Trefethen et al. 1993). The mechanism of the transient growth is essentially inviscid – the operators are highly non-normal for large Reynolds numbers and the transient growth is asymptotically large in $`Re`$. (This fact is extremely important in the case of Keplerian accretion disks, where the Reynolds number is literally astronomical: $`Re>10^{10}!`$)

These developments provoked a change of paradigm in the study of linear processes in the considered flows. It was the so-called nonmodal approach that became extensively used and even canonized in the 1990s on these grounds. As a result, substantial progress in the understanding of the shear flow phenomena has been achieved. The nonmodal analysis – some modification of the initial value problem – implicates the change of independent variables from a laboratory to a moving frame and the study of temporal evolution of spatial Fourier harmonics (SFH) of perturbations without any spectral expansion in time (see Sec. 4).

## 3 The similarity between Keplerian flow and plane parallel shear flow

Since our aim is to elucidate the basic similarity of the dynamics of plane shear flows and of Keplerian disks, we approximate the Keplerian flow in the two dimensional (2D) limit by its tangent plane parallel flow, while retaining the effect of rotation through the Coriolis force. This model, which ignores the purely geometrical complications, is known as the shearing sheet model; it has been has been used for the study of the linear dynamics of both wave and vortex mode 2D perturbations (cf. Goldreich & Lynden-Bell 1965; Goldreich & Tremaine 1978; Drury 1980; Nakagawa & Sekiya 1992 for wave mode perturbations, and Lominadze et al. 1988; Fridman 1989 and Ioannou & Kakouris 2001 for vortex mode perturbations).

The dynamical equations are written in the local co-moving Cartesian co-ordinate system:

$$xrr_0;yr_0(\varphi \mathrm{\Omega }_0t),$$
(1)
where ($`r,\varphi `$) are standard cylindrical co-ordinates and $`\mathrm{\Omega }_0`$ is the local angular velocity at $`r=r_0`$:

$$\mathrm{\Omega }(r)=\mathrm{\Omega }_0+\frac{\mathrm{\Omega }}{r}(rr_0)\mathrm{\Omega }_0+A\frac{x}{r_0}.$$
(2)
$`A(r_r\mathrm{\Omega })|_{r=r_0}`$ is the shear parameter (in Keplerian disks $`A=3\mathrm{\Omega }_0/2<0`$).

For the present purpose, we shall consider only 2D perturbations, independent on the axial coordinate $`z`$. The resulting dynamical linear equations for the perturbations of the radial velocity ($`u_r`$), the azimuthal velocity ($`u_\varphi `$) and the pressure ($`p`$) take the form:

$$\frac{u_r}{t}+Ax\frac{u_r}{y}2\mathrm{\Omega }_0u_\varphi =\frac{p}{x}+\nu \mathrm{\Delta }u_r,$$
(3)
$$\frac{u_\varphi }{t}+Ax\frac{u_\varphi }{y}+2\mathrm{\Omega }_0u_r+Au_r=\frac{p}{y}+\nu \mathrm{\Delta }u_\varphi ,$$
(4)
$$\frac{u_r}{x}+\frac{u_\varphi }{y}=0.$$
(5)
The incompressible limit is taken here to leave out wave mode perturbations and to keep only the vortex mode perturbations, which are the basic ingredient of the bypass scenario. In this case the energy density of the perturbation depends only on its kinematic characteristics and not on the thermodynamics, such as pressure.

Pressure terms are also absent when the dynamical equations are written in vorticity form (for $`\omega =\text{curl}𝐮`$), taking the curl of eqs. (3, 4). This is probably the reason why the pressure perturbations have been often ignored when discussing the dynamics of shear flows. However they play a very important role, even in the incompressible case, because they are the mediators of momentum exchange between fluid particles, which results in the transient growth of the vortex mode. (The physics of the “mediator activity” of the pressure perturbations is described in detail in Chagelishvili et al. 1993, 1996.) The importance of the pressure terms is clearly seen in the following analysis.

We multiply eqs. (3) and (4) respectively by $`u_r`$ and $`u_\varphi `$, in order to put them in kinetic energy form, and average them over a domain which is symmetrical in $`x`$. This procedure is similar to that performed in papers by Balbus et al. (1996) and Balbus & Hawley (1998).

$$\frac{}{t}\frac{u_r^2}{2}=2\mathrm{\Omega }_0u_ru_\varphi u_r\frac{p}{x}\nu |u_r|^2,$$
(6)
$$\frac{}{t}\frac{u_\varphi ^2}{2}=[2\mathrm{\Omega }_0+A]u_ru_\varphi u_\varphi \frac{p}{y}\nu |u_\varphi |^2.$$
(7)
The source of instability is represented by the shear term $`Au_ru_\varphi `$. But the Coriolis term $`2\mathrm{\Omega }_0u_ru_\varphi `$ is explicitly present too, and in eq. (7) it overbalances the source term, since $`A=3/2\mathrm{\Omega }_0`$ in the Keplerian disk. It thus would seem at first sight that the Coriolis force has a profound influence on the shear flow stability, as it has been argued by Balbus & Hawley (1996). But this conclusion is contradicted by the fact that the Coriolis terms disappear when summing eqs. (6) and (7), meaning that the growth rate of the total kinetic energy is independent of the Coriolis force.

A more detailed analysis will show that the Coriolis terms may be eliminated also on the level of the dynamical equations, by a suitable renormalization of the pressure perturbation. In this two-dimensional and incompressible case, the perturbation velocity field derives from the stream function $`\psi `$:

$$u_r=\frac{\psi }{y};u_\varphi =\frac{\psi }{x}.$$
(8)
Renormalizing the pressure perturbation

$$p^Rp2\mathrm{\Omega }_0\psi ,$$
(9)
eqs. 35, 6 and 7 can be rewritten as follows:

$$\frac{u_r}{t}+Ax\frac{u_r}{y}=\frac{p^R}{x}+\nu \mathrm{\Delta }v_r,$$
(10)
$$\frac{u_\varphi }{t}+Ax\frac{u_\varphi }{y}+Au_r=\frac{p^R}{y}+\nu \mathrm{\Delta }v_\varphi ,$$
(11)
$$\frac{u_r}{x}+\frac{u_\varphi }{y}=0.$$
(12)
$$\frac{}{t}\frac{u_r^2}{2}=u_r\frac{p^R}{x}\nu |u_r|^2,$$
(13)
$$\frac{}{t}\frac{u_\varphi ^2}{2}=Au_ru_\varphi u_\varphi \frac{p^R}{y}\nu |u_\varphi |^2.$$
(14)
In Sec. 2 we outlined the bypass concept in Cartesian/planar shear flows (a simple sketch of the bypass scenario is presented in the Appendix). Since we wish to apply a similar scenario to differentially rotating disks, we shall write the Cartesian counterparts of the averaged equations for comparison. If $`y`$ is the streamwise variable, and $`x`$ is the shearwise variable as in the considered case (i.e., $`x`$ is analog to $`r`$, and $`y`$ to $`\varphi `$), then the equations are:

$$\frac{}{t}\frac{u_x^2}{2}=u_x\frac{p^c}{x}\nu |u_x|^2,$$
(15)
$$\frac{}{t}\frac{u_y^2}{2}=A^cu_xu_yu_y\frac{p^c}{y}\nu |u_y|^2,$$
(16)
where $`A^c`$ is the shear parameter of the Cartesian flow $`U_0^c(0,A^cx)`$; $`u_x`$, $`u_y`$ and $`p^c`$ are the perturbations of velocity and pressure, respectively. If the flow is set up with $`U_0^c`$ decreasing with $`x`$ (i.e. $`A^c<0`$), then eq. (16) shows that steady flow is marked by transport in the direction of increasing $`x`$, as in the Keplerian disk.

These equations are identical to the normalized equations derived above for the Keplerian case, which stresses the similarity of these two flows. But they deal with spatially averaged quantities, and one may wish to establish that the similarity is even more profound, namely that the local kinematics and energy of 2D vortex perturbations of plane and disk flows are the same in time and space, starting from identical initial conditions.

## 4 Transient growth in 2D

Mainly for the purpose of illustrating the bypass mechanism, we shall describe now the evolution in time of an initial perturbation, in two dimensions. We perform a spatial Fourier transform of all relevant variables, as shown here for the pressure fluctuation:

$$p(x,y,t)=\widehat{p}(k_x(t),k_y,t)\mathrm{exp}\left[\mathrm{i}k_x(t)x+\mathrm{i}k_yy\right].$$
(17)
In a shearing flow, the perturbations cannot keep the form of a simple wave, since the wave-number of each spatial Fourier harmonics (SFH) depends on time (see Criminale & Drazin 1990 for a rigorous mathematical interpretation). In our geometry

$$k_x(t)=k_x(0)Ak_yt,$$
(18)
thus the wave-number of each SFH varies in time along the flow shear: it “drifts” in $`𝐤`$-space. Therefore $`p`$, in eq. (17) above, depends on time through $`k_x(t)`$.

Equations (5), (6) and (7) take the following form in Fourier space:

$$\frac{\mathrm{d}}{\mathrm{d}t}\widehat{u}_r(k_x(t),k_y,t)2\mathrm{\Omega }\widehat{u}_\varphi (k_x(t),k_y,t)=$$

$$=\mathrm{i}k_x(t)\widehat{p}(k_x(t),k_y,t),$$
(19)
$$\frac{\mathrm{d}}{\mathrm{d}t}\widehat{u}_\varphi (k_x(t),k_y,t)+(2\mathrm{\Omega }+A)\widehat{u}_r(k_x(t),k_y,t)=$$

$$=\mathrm{i}k_y\widehat{p}(k_x(t),k_y,t),$$
(20)
$$k_x(t)\widehat{u}_r(k_x(t),k_y,t)+k_y\widehat{u}_\varphi (k_x(t),k_y,t)=0.$$
(21)
For simplicity, we ignored here the viscous forces; their action may be easily incorporated in the analysis, as we shall see below.

One readily shows that these equations possess a time invariant

$$=k_y\widehat{u}_r(t)k_x(t)\widehat{u}_\varphi (t),$$
(22)
which expresses the conservation of vorticity in Fourier space. Making use of this invariant, we may write the solution of the system (1921) as follows:

$$\widehat{u}_r(k_x(t),k_y,t)=\frac{k_y}{k_x^2(t)+k_y^2},$$
(23)
$$\widehat{u}_\varphi (k_x(t),k_y,t)=\frac{k_x(t)}{k_x^2(t)+k_y^2},$$
(24)
$$\widehat{p}(k_x(t),k_y,t)=\frac{\mathrm{i}\left[Ak_y^2+2\mathrm{\Omega }_0(k_x^2(t)+k_y^2)\right]}{\left[k_x^2(t)+k_y^2\right]^2}.$$
(25)
We see that the perturbed quantities vary aperiodically in time (as is natural of vortex mode perturbations), and that they undergo transient amplification. Energy is exchanged between the perturbations and the background flow, and that represents the basis of the bypass transition scenario. It is important to note that the kinematic characteristics of the perturbations do not depend on the rotation rate $`\mathrm{\Omega }_0`$ – they are identical to those of Cartesian shear flows. Only the pressure perturbation depends on $`\mathrm{\Omega }_0`$. Consequently, the energetics of the SFH is identical to that of the Cartesian case:

$$\widehat{u}^2(k_x(t),k_y,t)\widehat{u}_{r}^{}{}_{}{}^{2}(k_x(t),k_y,t)+\widehat{u}_{\varphi }^{}{}_{}{}^{2}(k_x(t),k_y,t)=$$

$$=\frac{1}{k_x^2(t)+k_y^2}^2=\widehat{u}^2(k_x(0),k_y,0)\frac{k_x^2(0)+k_y^2}{k_x^2(t)+k_y^2}.$$
(26)
We see that the kinetic energy reaches maximum amplitude for $`k_x(t)=0`$. Similarly, the renormalized pressure perturbations do not involve $`\mathrm{\Omega }_0`$:

$$\widehat{p}^R(k_x(t),k_y,t)=\widehat{p}(k_x(t),k_y,t)2\mathrm{\Omega }_0\widehat{\psi }(k_x(t),k_y,t)=$$

$$=\frac{\mathrm{i}Ak_y^2}{[k_x^2(t)+k_y^2]^2}.$$
(27)
This expression is identical to that describing the evolution of the pressure perturbation in a Cartesian flow, which proves the similarity of transient growth in both types of flows, in two dimensions.

The above analysis has been done in the inviscid case. Accounting for the viscous forces is a straightforward procedure (see Fridman 1989): it may be done by multiplying the obtained equations by the factor $`\mathrm{exp}\left[\nu _0^tdt^{}(k_x^2(t^{})+k_y^2)\right]`$. For instance, eq. (19) and (26) will read as follows:

$$\widehat{u}_r(k_x(t),k_y,t)=\frac{k_y}{k_x^2(t)+k_y^2}e^{\nu _0^tdt^{}[k_x^2(t^{})+k_y^2]},$$
(28)
$$\widehat{u}^2(k_x(t),k_y,t)=\frac{^2}{k_x^2(t)+k_y^2}e^{2\nu _0^tdt^{}[k_x^2(t^{})+k_y^2]}.$$
(29)
To estimate the maximum amplification which is achieved in this transient growth at $`k_y^2k_x^2(0)`$, we assume that the largest perturbation wave-number corresponds to the Kolmogorov dissipation scale: $`|k_x(0)|r_0Re^{3/4}`$, where the Reynolds number is defined as $`Re=\mathrm{\Omega }_0r_0^2/\nu `$. Then

$$\frac{\widehat{u}_{\mathrm{max}}^2}{\widehat{u}^2(k_x(0),k_y,0)}\frac{k_x^2(0)}{k_y^2}\left(\frac{Re}{(k_yr_0)^2}\right)^{3/2},$$
(30)
which can reach a huge value in astrophysical disks.

## 5 Discussion and conclusions

We have seen that the concept of nonlinear instability in smooth spectrally stable shear flows has undergone substantial revision in the hydrodynamic community. Nowadays the concept of the bypass transition to turbulence has become favorite and is under intensive development.

Let us summarize the main features of this concept:
– the onset of turbulence and the turbulent state itself in smooth spectrally stable shear flows is supported energetically by the linear transient growth of vortex mode perturbations – the key ingredient of the turbulence is the vortex mode (eddy) perturbation, and the key phenomenon is the linear transient growth of the perturbation;
– nonlinear processes do not contribute to any energy growth, but regenerate vortex mode perturbations that are able to extract shear flow energy. Doing so, nonlinear processes only indirectly favor the energy extraction by the vortex mode perturbations;
– the non-orthogonal nature of the linearized Navier-Stokes equations is the formal basis of the transient growth;
– the non-orthogonal nature increases with increasing Reynolds number; thus the operators are highly non-normal for the huge Reynolds numbers of Keplerian disks $`(Re>10^{10})`$ and the transient growth is asymptotically large in $`Re`$.

In this paper we wanted to contribute to the revival of hydrodynamic shear turbulence as a possible explanation for the “anomalous” viscosity in non-magnetized Keplerian accretion disks. Specifically, we wished to draw the attention of astrophysicists to the bypass concept.

We showed that in disk flows there exist vortex mode perturbations that are similar to those which are held responsible for the onset of turbulence in the Cartesian shear flow. The key point of our analysis is the interpretation of the important role of the pressure perturbations in the dynamical processes.

In fact, the kinematics and energetics of the vortex mode perturbations are identical in the rotating disk and the Cartesian shear flows in the 2D case. The Coriolis force only causes deviation of the pressure perturbations from the Cartesian case. By focusing on the epicyclic motions, and underestimating the action of the pressure terms, one is led to the false conclusion that the Coriolis force suppresses hydrodynamic turbulence in Keplerian flows (Balbus et al. 1996; Balbus & Hawley 1998).

This property of the pressure perturbations has been established here in the 2D case. In the more general case of 3D perturbations, the dynamics of vortex mode perturbations is somewhat more complicated, as it will be shown in a forthcoming paper (Tevzadze et al.). But the role of the pressure perturbations is still to counteract the Coriolis force. In the 3D disk case, the wave-number domain where the perturbation undergoes transient growth is smaller in comparison to 3D Cartesian flows, but this is compensated by the very large Reynolds number characterizing astrophysical disks.

Other scenarios which may lead to hydrodynamic turbulence in astrophysical disks have been presented recently: they invoke a linear (spectral) instability arising from the stratification perpendicular to the disk (Dubrulle et al. 2002) or due to deviations from cylindrical rotation (Urpin 2002; Klahr & Bodenheimer 2002). Such instabilities may well compete with the bypass mechanism presented here, and it is not possible to conclude presently which is the best candidate for rendering astrophysical disks turbulent.

In laboratory experiments the field narrows, because there is no stratification in the fluid (other than imposed on purpose), and there the bypass mechanism provides an attractive explanation for the turbulence detected in flows which are linearly stable (angular momentum increasing outwards). We refer to Couette-Taylor experiments performed by Wendt (1933), Taylor (1936), Coles (1950) and Van Atta (1966). Very recently such turbulence has been observed also in rotation profiles which share the properties of Keplerian disks, namely with their angular velocity decreasing outwards (Richard 2001).

Decisive conclusions about the self-sustenance of the turbulent state need to be supported by numerical simulations. To our knowledge, the simulations reported so far have failed to detect hydrodynamic shear turbulence in rotating flows with angular momentum increasing outwards. In a recent paper, Longaretti (2002) has discussed the possible explanations of this negative result. The main reason is probably the lack of spatial resolution, which still prevents from reaching even the Reynolds numbers at which turbulence is detected in the laboratory. At this stage, we hope that the analysis presented here will help to build a fruitful background for the nonlinear numerical stability analysis of rotating astrophysical shear flows. In any way, the definite answer will result from the convergence of theoretical, numerical and experimental work.

###### Acknowledgements.

This work is supported by the International Science and Technology Center grant G-553. G.D.C. would like to acknowledge the hospitality of Observatoire de Paris (LUTH). The authors would like to thank the referee for his helpful comments on the early version of the paper.

## Appendix A Bypass scenario of the transition to turbulence

We present here a simple sketch of the bypass scenario applied to Keplerian disk flow, in the wave-number plane $`(k_x,k_y)`$ (see Fig. 1). We shall define as “active domain” the region where viscous dissipation may be neglected, i.e. where $`k_x^2+k_y^2<k_\nu ^2`$, with $`k_\nu Re^{3/4}/r_0`$.

The linear dynamics of the perturbation may be described by following each of its spatial Fourier harmonics (SFH) in this wave-number plane. We single out one that is located initially at some point 1 in the amplification part of the “active domain”, which meets the condition $`k_x/k_y<0`$. According to eq. (18), as $`k_x(t)`$ varies in time, the SFH drifts in the direction marked by the arrows, since $`A<0`$. (We present the drift of the SFH only in the upper half-plane; since the perturbation is real, there is a counterpart in the lower half-plane.) Initially, as $`|k_x(t)|`$ decreases, the energy of the SFH grows. This growth lasts until the wave vector crosses the line $`k_x=0`$ (point 2). Then, while undergoing attenuation, the SFH continues its drift until it reaches point 3, where it is dissipated through viscous friction. The same will occur with all other Fourier harmonics. Consequently, if the nonlinear interaction between different Fourier harmonics is inefficient, the perturbation disappears eventually. Permanent extraction of shear energy by the perturbations is necessary for their maintenance, which is possible when quadrants II and IV where $`k_x/k_y<0`$ are being repopulated through nonlinear interaction between Fourier harmonics of quadrants I and III which have reached sufficient amplitude. This may be achieved through three-wave processes $`𝐤_\mathrm{𝟏}+𝐤_\mathrm{𝟐}𝐤`$; four wave processes $`𝐤_\mathrm{𝟏}+𝐤_\mathrm{𝟐}+𝐤_\mathrm{𝟑}𝐤`$; five wave processes, etc. In the figure we present an example of three wave process $`𝐤_\mathrm{𝟏}+𝐤_\mathrm{𝟐}𝐤`$ that contributes to the regeneration of an SFH in the amplification area, transferring perturbation energy to it from the attenuation areas. The bypass scenario implies the dominance of this regeneration tendency of nonlinear processes, i.e. predominant transfer of perturbation energy from quadrants I and III to quadrants II and IV, in other words nonlinear positive feedback. To what extent the reproduction of the Fourier harmonics in quadrants II and IV is sustained, even in the case of positive feedback, depends both on the amplitude and on the spectrum of the initial perturbation. Nonlinear decay processes are weak at low amplitudes and are not able to compensate the linear drift of SFH in k-plane. As a result, weak perturbations are damped without any trace, and without inducing transition to turbulence. The higher the amplitude of initial perturbation, and the stronger are the nonlinear effects. At a certain amplitude (which, of course, depends on the initial vortex perturbation spectrum in k-plane and on the Reynolds number), the nonlinear processes are able to compensate the action of the linear drift, thus ensuring permanent return of SFH to the amplification areas (quadrants II and IV). This eventually ensures a permanent extraction of energy from the background flow and maintenance of the perturbations, and thus of turbulence.

Therefore the bypass scenario can be realized only in the case of finite amplitude perturbations and in each case it has a threshold that depends on the perturbation spectrum and the Reynolds number.

Let us summarize. According to the bypass concept, vortex mode (eddy) perturbations are the basic ingredient of hydrodynamic shear turbulence. The bypass scenario involves the interplay of four basic phenomena:
– the linear “drift” of SFH in the k-plane;
– the transient growth of SFH;
– the usual viscous dissipation;
– nonlinear processes that close the feedback loop of the transition by mixing – by the angular redistribution of SFH in k-plane.