text:
# Structural Parameters and Dynamical Masses for Globular Clusters in M33 1footnote 11footnote 1Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration.

## 1 Introduction

In the Local Group, by far the best studied globular cluster (GC) systems are those of the Milky Way and M31. Significant amounts of data have been gathered for globular clusters in these galaxies, with extensive compilations in Harris (1996) (Milky Way, 147 clusters) and Barmby et al. (2000) (M31, 435 clusters). Both of these large spirals are known to contain two distinct GC subpopulations: a metal-poor halo population and a more metal-rich population, associated with the bulge or thick disk (Kinman, 1959; Zinn, 1985; Minniti, 1995; Jablonka et al., 1998; Perrett et al., 2002). The presence of two GC subpopulations presumably indicates different formation processes, and may hold clues to the formation and early evolution of the host galaxies. A particularly controversial issue is whether or not age differences exist between globular clusters of different metallicities. In M31, Barmby & Huchra (2000) found evidence for age differences of 4–8 Gyrs between the metal-rich and metal-poor clusters and Sarajedini, Chaboyer, & Demarque (1997) concluded that a similar age spread exists among GCs in the Milky Way. Other recent studies have proven less conclusive, but age differences of a few Gyrs cannot be ruled out (Stetson, Vandenberg, & Bolte, 1996; Buonanno et al., 1998; Rosenberg et al., 1999).

The GC system of the third spiral galaxy in the Local Group, M33, is more sparse and consequently more difficult to identify and study. Early work on the star clusters in this galaxy was done by Hiltner (1960) and Kron & Mayall (1960), who noted that some of them had unusually blue colors. Photometry for about 130 cluster candidates was given by Christian & Schommer (1982, 1988), who tentatively identified 27 old GC candidates in M33. They also noted that M33 appears to contain a number of “blue populous” clusters similar to those in the Large Magellanic Cloud (Hodge, 1961), but without the pronounced age gap that exists between $`3`$ Gyrs and 12–15 Gyrs in the LMC (Girardi et al., 1995). The old GCs exhibit halo-like kinematics, while the younger clusters follow the rotation of the disk (Schommer et al., 1991). Color-magnitude diagrams for 10 halo GCs in M33, obtained with WFPC2 on board HST, were presented by Sarajedini et al. (1998), who noted that most of the clusters showed no blue horizontal branches in spite of having low metallicities. They suggested that this might be due to a second-parameter effect, indicating that many of the “old” M33 clusters may be several Gyrs younger than their Galactic counterparts. Additional clusters have been identified on WFPC2 images by Chandar, Bianchi, & Ford (2001), and again Chandar et al. (2002) suggested a larger age spread among the halo GCs compared to the Milky Way.

It is important to establish how much the old globular clusters in different galaxies have in common. A larger age spread among halo GCs in M33, for example, might indicate that this galaxy assembled over a longer period of time, while the structure of individual clusters may provide information about the gas clouds out of which they formed (Murray & Lin, 1992; Elmegreen & Efremov, 1997). In the Milky Way, a number of fundamental correlations are known to exist between parameters such as surface brightness, core radius and velocity dispersion of globular clusters (Kormendy, 1985; Djorgovski, 1995). McLaughlin (2000) found that all Milky Way GCs are consistent with a constant $`V`$–band core mass-to-light ratio (within the error margins) of 1.45 $`\text{M}\text{}\text{L}\text{}^1`$. He also found a tight relation between the binding energy of individual clusters and the total luminosity, with a weak decrease in binding energy as a function of Galactocentric distance. The dependence on Galactocentric distance may hold information about pressure gradients in the early Galactic halo.

Detailed information about the structure of extragalactic GCs, even within the Local Group, is more limited for obvious reasons. At the distance of M31 and M33, a typical ground-based resolution of $`1\mathrm{}`$ corresponds to a linear scale of about 4 pc, roughly similar to the typical half-light radii of globular clusters. Thus, information about the detailed structure of individual clusters is not easily obtained. However, with the HST the situation is greatly improved, with one pixel on the Planetary Camera (PC) corresponding to about 0.2 pc. This is of great importance not only for studies aiming at a morphological characterization of the clusters, but also for dynamical studies in which velocity dispersions must be tied to knowledge about the cluster structure in order to provide accurate mass-to-light ratios. So far, the mass-to-light ratios and other properties of globular clusters in M31 appear to be identical to those observed in Milky Way GCs (Djorgovski et al., 1997; Dubath & Grillmair, 1997; Barmby, Holland, & Huchra, 2002). Similarity in the M/L ratios (modulo metallicity effects) is, of course, to be expected if the clusters have similar ages, unless there are variations in the stellar mass function. No information about dynamical masses has yet been published for M33 GCs, but could potentially be used as a tool to independently check whether age differences exist.

In this paper we study 4 globular clusters in M33, selected from the sample of Sarajedini et al. (1998). The clusters were originally selected based on their halo-like kinematics and red colors ($`\mathrm{B V}>0.6`$), and should be as close an analogy to the halo clusters in the Milky Way as possible. We use new high-dispersion spectroscopy from the HIRES spectrograph on the Keck I telescope to measure velocity dispersions for the individual clusters, and combine these with structural information from the same HST images used by Sarajedini et al. (1998) to derive mass-to-light ratios and compare with data for Milky Way globular clusters. Throughout the paper we will adopt a distance modulus of 24.84 and the reddenings for each cluster determined by Sarajedini et al. (2000).

## 2 Data

Observations of the four clusters (R12, H38, M9 and U49) were carried out on Oct 24 and 25, 1998 with the HIRES spectrograph (Vogt et al., 1994) on the Keck I telescope. A slit width of 0$`\stackrel{}{\mathrm{.}}`$725 was used, providing a spectral resolution of $`\lambda /\mathrm{\Delta }\lambda 54000`$. The wavelength range was 3730–6170 Å, distributed over 38 echelle orders, although a useful S/N was achieved only for $`\lambda 4000`$ Å. For each cluster, 7–8 individual exposures of 1800 s each were obtained. In addition, the star HD 1918 (G9 III) was observed as a template star for the velocity dispersion measurements.

Extraction of the spectra from the CCD images was done with the highly automated MAKEE package, written by T. Barlow and specifically tailored for reduction of HIRES data. MAKEE automatically performs bias and flatfield corrections, cleans the images of cosmic-ray (CR) hits and then locates and extracts each echelle order from the CCD images. Wavelength calibration was done using spectra of ThAr calibration lamps mounted within HIRES, with zero-point corrections based on sky lines. The individual 1d spectra of each cluster were co-added using a sigma-clipping algorithm to eliminate any residual CR events, producing a S/N for the co-added spectra of 15–25 per resolution element.

### 2.1 Velocity dispersions

A variety of techniques have been developed to measure velocity dispersions from integrated spectra in cases like ours where the expected dispersions (5–10 km/s) are comparable to the instrumental resolution. The most direct method is to simply convolve the template star spectrum with Gaussian profiles corresponding to different velocity dispersions, and then find the velocity dispersion that minimizes the residuals when the (scaled and smoothed) template star spectrum is subtracted from the cluster spectrum. One appealing aspect of this method is that the quality of the fit and the match between template and object spectra can be readily inspected, but the velocity dispersions are quite sensitive to an accurate determination of the velocity difference between the two objects as well as a good match between the cluster and template spectra. Alternatively, the fitting can be done in the Fourier domain (Illingworth, 1976), utilizing the fact that the slope of the power spectrum is strongly sensitive to the spectral resolution. A third method is to cross-correlate the cluster spectrum with the template spectrum and use the width of the cross-correlation peak as an indicator of the velocity dispersion (Tonry & Davis, 1979; Ho & Filippenko, 1996). Yet another technique was used by Dubath & Grillmair (1997), who convolved the cluster spectrum with a suitably designed mask to obtain a “mean” spectral line whose width is related to the velocity dispersion.

The method of directly fitting a smoothed template spectrum to the cluster spectrum is illustrated in Figure Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration., where we compare two echelle orders of the spectrum of cluster H38 with the template star. Each panel shows the cluster spectrum, the template star spectrum, the best-fitting smoothed template star spectrum and the residuals. Even though the smoothing results in a much better fit than the raw template spectrum, indicating that the broadening of lines in the cluster spectrum due to velocity dispersion is well resolved, significant residuals are still visible for many of the stronger lines. There are no systematic trends in the residuals (some are positive, others negative), but the template star is evidently not an ideal match to the cluster spectrum.

We also measured velocity dispersions with the cross-correlation method, using the FXCOR task in the RV package in IRAF222 IRAF is distributed by the National Optical Astronomical Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc. under contract with the National Science Foundation . This method is less dependent on an exact match between the cluster and template spectra, because a poorer match will tend to change the amplitude rather than the width of the cross-correlation peak. In Fig. Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. we show the cross-correlation functions (CCFs) for cluster H38 (order 28) vs. the template star (solid line), as well as the CCFs for the template star spectrum broadened with three different Gaussians (0, 6 km/s, 12 km/s) vs. the template star spectrum itself (dashed lines). The best fit is obtained for a velocity dispersion around 6 km/s. As a useful “byproduct”, the cross-correlation technique automatically provides radial velocities for the clusters, which are needed as input for the other methods.

Line-of-sight velocity dispersions ($`\sigma _x`$) are listed in Table 1 for both the direct fitting and cross-correlation techniques. The $`\sigma _x`$ values are averaged over all fitted echelle orders, with the number of echelle orders used for the fits listed in column (2). We give both the rms scatter of the $`\sigma _x`$ values measured on the individual echelle orders, as well as the formal standard errors on the mean. However, many of the uncertainties, related to choice of template star, fitting technique etc. are systematic rather than random, and simply estimating the uncertainties from the scatter between the various echelle orders may not give a realistic estimate of the true uncertainties. Table 1 also lists the radial velocities of the clusters and measurements from Schommer et al. (1991) and Chandar et al. (2002) for comparison. The dominating uncertainty on our radial velocity measurements is actually the radial velocity of the template star itself, given as $`+35.7\pm 2`$ km/s in SIMBAD. For all clusters except H38, our radial velocities agree with those of Chandar et al. (2002) within the errors. The mean difference between our measurements and theirs is $`19\pm 13`$ km/s.

As seen in Table 1, the cross-correlation technique gives systematically lower velocity dispersions than the direct fitting, with differences of 0.4 – 0.9 km/s. We tested the direct fitting method by artificially broadening a template star spectrum, adding noise and applying the fitting algorithm to the resulting spectrum. The dispersions of the Gaussians used for the smoothing were reproduced with an accuracy of better than 0.1 km/s for a S/N similar to that of the cluster spectra. For even lower S/N the accuracy decreased, but no systematic trends were evident. However, we note that one pixel corresponds to a velocity difference of about 2.0 km/s, so the difference between the two methods may well be due to centering and/or binning effects. Direct fitting is more sensitive to such effects than the cross-correlation technique and we may assume that the latter is more accurate. In the following discussion we therefore adopt the velocity dispersions obtained by the cross-correlation method.

### 2.2 Structural parameters

To obtain structural parameters for the clusters we used images from the Wide Field Planetary Camera 2 on board the Hubble Space Telescope, described in Sarajedini et al. (2000). The clusters are in all cases roughly centered on the Planetary Camera (PC) chip. Each dataset contained exposures in the F555W ($`V`$) and F814W ($`I`$) bands, allowing for two independent sets of measurements for each cluster.

From visual inspection of the WFPC2 images (Figs 1, 2, 4 and 6 in Sarajedini et al. 2000) it is clear that none of the clusters are strongly elongated. We used the ELLIPSE task in the STSDAS package in IRAF to fit elliptical isophotes to the cluster images, after filtering the images with an $`11\times 11`$ pixels median filter to create the smooth profiles required by ELLIPSE. Fig. Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. shows the position angle (N through E) and ellipticity as a function of semi-major axis for each cluster and confirms that all clusters have small ellipticities. For radii less than about $`0\stackrel{}{\mathrm{.}}5`$, the smoothing and random fluctuations due to individual stars make the fits meaningless. For H38 and M9, strong variations in the position angle are seen also at larger radii, most likely because the ellipticities are so small that the PA is very poorly constrained. U49 is the only one of the four clusters where visual inspection of the images hints at some elongation, and the ELLIPSE fits yield an ellipticity of $`ϵ`$ 0.10–0.15 for this cluster, with a relatively stable PA around 30 degrees. For the remaining clusters the ellipticities are less than $`0.1`$, consistent with the mean value of $`0.07\pm 0.01`$ for Milky Way globulars (White & Shawl, 1987) and also in good agreement with the mean $`ϵ=0.11\pm 0.01`$ for M31 GCs (Barmby, Holland, & Huchra, 2002).

Additional structural parameters for the clusters were derived by fitting single-mass King–Michie models (King, 1966, hereafter simply “King” models) to the WFPC2 images. Because the ellipticities are generally small and the position angles ill-determined and possibly varying with radius, we simply assumed circularly symmetric profiles for the model fits. We then carried out a least-squares fit directly to the images, solving for the core radius $`r_c`$, concentration parameter $`W`$ and central surface brightness $`\mu _0`$. The minimization was done with the downhill simplex (“amoeba”) algorithm described in *Numerical Recipes* (Press et al., 1992). Although the clusters are quite well resolved, crowding near the center and strongly variable completeness functions as a function of radial distance made direct starcounts unfeasible. For the same reasons, attempts to subtract individual stars from the images might also lead to systematic errors in the derived profiles. Assuming that field stars are uniformly distributed throughout the frames, we decided that the approach that would be least likely to suffer from systematic effects due to partial resolution of the clusters was to simply fit the King models to the images directly. The background level was determined as the average of all pixel values in an annulus centered on the cluster, with an inner radius of 200 pixels (9 arcsec) and 100 pixels wide. This is about the largest possible background annulus that could be consistently used, considering variations in the exact position of the clusters on the PC chip and the limited field size. The fit itself was carried out for radii less than 200 pixels. We also attempted to solve for the centroid of the profile simultaneously, but these fits turned out to be unstable and we found that a more reliable approach was to determine the center of the profiles by carrying out the fits for various central positions and inspecting the residuals. Even for fairly small centering errors ($`1`$ pixel), asymmetries in the residuals were clearly visible.

As noted by other authors (Chandar, Bianchi, & Ford, 2001; Barmby, Holland, & Huchra, 2002), the effect of the HST PSF is not entirely negligible when measuring structural parameters for star clusters in M31 or M33. Our modeling included a convolution of the 2-d King profiles with the PC PSF, generated by the TINYTIM software333 Available at the URL http://www.stsci.edu/software/tinytim/tinytim.html  and, like previous studies, we found small but measurable differences in the fits when the PSF convolution was included. Core and effective radii were typically $`0.1`$ pc ($`0\stackrel{}{\mathrm{.}}02`$) smaller, while the central surface brightnesses $`\mu _0`$ increased by 0.1–0.2 mag with PSF convolution.

Fits to the F555W images are shown in Fig. Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration.. The error bars indicate the actual cluster profiles while the solid lines are the King model fits. The light profiles were measured in concentric apertures with the PHOT task within DAOPHOT in IRAF, measuring the sky background in the same annuli that were used for the fits. Note that, because of the 2-d fitting procedure, the solid lines are *not* direct fits to the actual datapoints in Fig. Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration., but are measured on the best-fitting convolution product of a 2-d King profile and the PC PSF, using the same procedure as for the cluster profiles. The King models generally provide excellent fits to the globular cluster profiles. In particular, we see no evidence for central cusps that might indicate core-collapsed clusters. The fitted parameters are listed in Table 2, with separate fits for the F555W and F814W images. Uncertainties were estimated by repeating the fitting procedure 18 times using different initial values for $`r_c`$, $`W`$ and $`\mu _0`$. We have adopted the standard deviation on each fitted parameter, rather than the standard error on the mean, as our uncertainty estimates. Alternatively, the difference between the F555W and F814W fits may be used as an estimate of the uncertainties, and provides similar uncertainty estimates to those in Table 2.

Cluster R12 was also measured by Chandar, Bianchi, & Ford (1999). Although they fitted the empirical version (King, 1962) of the King models rather than the theoretical King–Michie models used here, the two types of King models are structurally similar enough that a comparison of the derived core radii is illustrative. Fitting a King (1962) profile to the F555W image, Chandar et al. found a core radius of $`0\stackrel{}{\mathrm{.}}234`$ for R12 which is only slightly larger than our value of $`0\stackrel{}{\mathrm{.}}225`$ (corrected for the PSF). Without correction for the PSF we get $`r_c=0\stackrel{}{\mathrm{.}}253`$, now slightly larger than the Chandar et al. value. The small differences between our measurements and those of Chandar et al. can easily be ascribed to the different choices of King profiles and fitting procedure, small centering errors etc.

In addition to the three basic parameters $`r_c`$, $`W`$ and $`\mu _0`$, Table 2 also includes the parameter $`C`$, which is the ratio of the tidal radius $`r_t`$ to the core radius $`r_c`$, as well as the effective radius $`r_\mathrm{e}`$, defined as the radius containing half the cluster light in projection. These two additional parameters both follow directly from the King model $`r_c`$ and $`W`$ parameters, but are listed for convenience.

One potential concern is that the tidal radius for one of the clusters, U49, is larger than the inner radius of the sky annulus, leading to possible systematic errors in the sky background determination. For the best-fitting model, a sky background level of about 0.2 ADU is expected at $`r=200`$ pixels, decreasing rapidly outwards. We tried varying the background level by $`\pm 0.5`$ ADU and redid the fits for U49. The core and effective radii were affected by less than $`0\stackrel{}{\mathrm{.}}05`$ (0.2 pc) while the dimensionless $`C`$ parameter changed by about 1, i.e. the changes were in both cases comparable to the uncertainties in Table 2.

The core and effective radii for the M33 clusters are well within the range spanned by Milky Way globular clusters. The $`r_c`$ distribution for Milky Way GCs peaks at about 1 pc, but has a tail extending up to above 10 pc. Excluding core-collapsed clusters, the median value is 1.3 pc. The median half-light radius $`r_\mathrm{e}`$ for Milky Way GCs is 3.3 pc (using values from Harris, 1996), but the distribution is again strongly asymmetric. It may be worth noting that the reddest M33 cluster (R12, see sect 2.3) also has the smallest effective radius, an effect that has been observed in many other galaxies including the Milky Way and M31 (e.g. Barmby, Holland, & Huchra, 2002; Larsen et al., 2001). Another point of interest is that previous studies of Milky Way and M31 globular clusters have found that low-concentration clusters tend to be more elliptical (Barmby, Holland, & Huchra, 2002), which is consistent with U49 being the most elongated as well as the largest of the clusters in our sample. However, the 4 clusters studied here clearly constitute a too small sample to draw conclusions as to whether or not the same effects are generally present in M33.

### 2.3 Integrated photometry

Fig. Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. shows the integrated $`V`$ magnitudes of the clusters as a function of radius. The curves-of-growth are drawn with dashed lines beyond the tidal radii obtained from the King profile fits. The background level was determined in the same way as for the King profile fits and the number counts in F555W and F814W in each aperture were converted to $`V`$ magnitudes and $`VI`$ colors using the transformations in Holtzman et al. (1995). Because the Holtzman et al. zero-points include an implicit $`0.1`$ mag aperture correction from their $`r=0\stackrel{}{\mathrm{.}}5`$ reference aperture to infinity, which does not apply in our case, we have added 0.1 mag to the magnitudes. Apart from small color terms, this gives exactly the same results as using the photometric zero-points in the image headers directly.

The curve-of-growth for M9 has a “bump” at $`r=200`$ pixels, due to a bright star that is located just within the outer aperture radius, but the curves are otherwise smooth and monotonically increasing over most of the radial range. For a couple of clusters (R12 and H38), the curves decrease slightly at large radii, reflecting the inherent uncertainties in the background determination. A small decrease of 0.5–1.0 ADU in the background levels would make the curves-of-growth for R12 and H38 increase monotonically out to their tidal radii.

Both Christian & Schommer (1988, hereafter CS88) and Chandar, Bianchi, & Ford (2001, CBF01) have published integrated photometry for the 4 clusters, using ground-based CCD imaging and WFPC2 data, respectively. In Table 3 we compare our photometry with CS88 and CBF01. CS88 used an aperture radius of $`3\stackrel{}{\mathrm{.}}7`$, corresponding to 82 pixels on the PC camera, while CBF01 used an aperture of $`2\stackrel{}{\mathrm{.}}2`$ (48 pixels) for H38, M9 and R12 and $`2\stackrel{}{\mathrm{.}}7`$ (59 pixels) for U49. We redid the photometry in the same apertures used by CBF01 and CS88 as well as in a larger aperture. For R12 and H38 this larger aperture coincided with the tidal radii of the clusters, for M9 we used an aperture of $`r=190`$ pixels (to avoid the bright star) and for U49 we used $`r=200`$ pixels. For the comparison with CBF01 we used the same background annulus as they did ($`3\stackrel{}{\mathrm{.}}5<r<5\mathrm{}`$), but for the other measurements the same annulus was used as for the King model fits and the curves-of-growth.

While our $`VI`$ color measurements agree fairly well with those of CBF01, our $`V`$ magnitudes are systematically fainter than theirs by about 0.18 mag. Part of this discrepancy is due to the fact that no correction was applied to the Holtzman et al. zeropoints in CBF01 (Chandar, priv. comm.). By experimenting with the PHOT options, we found that the remaining 0.08 mag can be accounted for by differences in the background level estimates. Because the resolution into individual stars makes the sky histogram highly asymmetric, the background measurements in these fields are very sensitive to the choice of algorithm. Specifically, if we use the MODE option (as did CBF01) instead of simply calculating the mean, the discrepancy is reduced to a few times 0.01 mag. However, because the contribution from individual, resolved stars is part of the “sky” background within the photometric apertures and therefore needs to be subtracted, a simple mean appears to be the most appropriate background estimator for our purpose.

Our photometry in the $`3\stackrel{}{\mathrm{.}}7`$ aperture agrees closely with that of CS88, being only 0.02 mag fainter on the average. For some clusters the curves-of-growth continue to rise well beyond $`3\stackrel{}{\mathrm{.}}7`$, and the measurements in Table 3 suggest that this aperture may underestimate the luminosity of the most extended cluster, U49, by $`0.2`$ mag. However, measurements in the largest apertures are quite sensitive to uncertainties in the background level, and for R12 and H38 the integrated magnitudes in the largest aperture are actually slightly fainter than in the $`3\stackrel{}{\mathrm{.}}7`$ aperture, in accordance with the behavior of the curves-of-growth pointed out above. The small adjustments in the background level that would make the curves-of-growth monotonic would affect the photometry at $`r=50`$ pixels ($`2\stackrel{}{\mathrm{.}}3`$) by only $`0.02`$ mag, but at the tidal radii the changes amount to 0.1–0.2 mag. This may provide a more realistic estimate of the true uncertainty on the integrated photometry than pure photon statistics, which give formal errors of less than 0.002 mag. In practice, the uncertainty on the background level is dominated by stochastic variations in the numbers of individual, bright stars within the aperture annulus, as well as true variations due to dust obscuration etc. A larger field of view would help reduce the stochastic variations, but would increase the sensitivity to large-scale variations in the background. For the following discussion we will use our photometry in columns (10) and (11) of Table 3 as the best approximation to the total integrated magnitudes of the clusters.

## 3 Masses and mass-to-light ratios

Once the velocity dispersion and structural parameters of a cluster are known, the total mass can be estimated. The simplest way is to use the virial theorem which implies

$$M_{\mathrm{vir}}=a\frac{\sigma _{3D}^2r_h}{G}$$
(1)
where the constant $`a`$ has a value of 2.5 (Spitzer, 1987, p. 11). In this expression, $`\sigma _{3D}`$ is the three-dimensional velocity dispersion, $`\sigma _{3D}^2=3\sigma _x^2`$, and $`r_h`$ is the 3-dimensional half-*mass* radius which is generally larger than the projected half-*light* radius $`r_\mathrm{e}`$ by a factor of 4/3. In the present case, however, the King model fits allowed $`r_h`$ to be directly evaluated for each individual cluster, without using $`r_\mathrm{e}`$ as an intermediate step. An alternative method is to use “Kings formula” (Richstone & Tremaine, 1986; Queloz, Dubath, & Pasquini, 1995),

$$M_{\mathrm{king}}=\frac{9}{2\pi G}\frac{\mu r_c\sigma _0^2}{\alpha p}$$
(2)
Here $`\sigma _0`$ is the *central* (projected) velocity dispersion and $`r_c`$ is the core radius of the cluster. The remaining parameters, $`\mu `$ (a dimensionless mass, not to be confused with the central surface brightness $`\mu _0`$), $`p`$ (ratio of central surface density to central 3-d density) and $`\alpha `$ (ratio of half-width at half maximum to core radius, usually close to 1) are all functions of the concentration parameter and were again evaluated using the King model fits.

The velocity dispersion $`\sigma _{3D}`$ used in (1) is the mean value averaged over the entire cluster, while the $`\sigma _0`$ used in (2) is the central value. Neither of these is the actual observed value $`\sigma _x`$, since a slit of finite size was used for the observations. Using the King model fits, which provide the velocity dispersion in the cluster as a function of radius, we estimated the aperture corrections from our slit to the central and global values, assuming a slit width of 0.72″ (3.2 pc at the distance of M33) and infinite length. The aperture corrections are relatively modest, generally amounting to less than 10% in either direction.

It should be noted that both methods implicitly assume a constant M/L ratio as a function of radius and isotropic velocity distributions. In practice, some mass segregation is expected because of equipartition of energy, and has indeed been detected observationally in some cases (e.g. King, Sosin, & Cool, 1995; Ferraro et al., 1997; Sosin, 1997). The King models used here are single-mass models and certainly an oversimplification of the true dynamical situation in the clusters, but in practice the difference between dynamical masses obtained by fitting single-mass King models and more sophisticated methods tends to be relatively small (e.g. Djorgovski, 1995; Meylan et al., 2001). Anisotropy is expected to occur mostly in the outer parts of the clusters (Gunn & Griffin, 1979; Lupton, Gunn, & Griffin, 1987; Takahashi & Lee, 2000) and should not have a strong impact on M/L ratios derived from the integrated light.

The computed central and global ($`\sigma _{\mathrm{}}`$) velocity dispersions are listed in Table 4, col. (2) and (3). Columns (4) and (5) give the corresponding King and virial masses, using the structural parameters measured on the F555W frames. Considering that neither the rms nor the standard errors on the mean are likely to provide realistic estimates of the true errors on the velocity dispersions, we have adopted 0.5 km/s for the errors on $`\sigma _x`$, corresponding roughly to the typical difference between the measurements based on direct fit and cross-correlation. We have omitted errors on $`\sigma _{\mathrm{}}`$ and $`\sigma _0`$ in Table 4 since these scale, at least to first order, with the errors on $`\sigma _x`$. Column (6) lists the central density in $`\text{M}\text{}\mathrm{pc}^3`$ which, per definition, is $`\rho _0=M_{\mathrm{King}}/(\mu r_c^3)`$. The $`V`$-band (King) mass-to-light ratios, $`M/L_V`$, are given in col. (7) and column (8) lists the binding energy of each cluster, calculated as in McLaughlin (2000). For the mass-to-light ratios we have used the $`V`$ magnitudes in Table 3 and the distance modulus of 24.84 and the reddenings of individual clusters given in Sarajedini et al. (2000). Other estimates of the distance modulus of M33 are $`24.52\pm 0.14`$(random)$`\pm 0.13`$(systematic) (Lee et al., 2002, based on Cepheids), $`24.81\pm 0.04`$(random)$`\pm 0.05`$(systematic) (Kim et al., 2002, red giants), and $`24.64\pm 0.09`$ (Freedman, Wilson, & Madore, 1991, Cepheids). The derived M/L ratios are inversely proportional to the assumed distance, and could be underestimated by up to 14% if the smallest of the above distances is correct.

As seen from Table 4, the virial and King masses are generally very similar. This agreement is not entirely fortuitous, since the two estimates are not independent. In fact, the ratio $`M_{\mathrm{King}}/M_{\mathrm{vir}}`$ is

$$\frac{M_{\mathrm{King}}}{M_{\mathrm{vir}}}\frac{\mu }{\alpha p}\frac{r_c}{r_h}\frac{\sigma _0^2}{\sigma _{3D}^2}$$
(3)
where the constants have been excluded. Since all the factors in Eq. 3 are functions of the concentration parameter only, this is also true for $`M_{\mathrm{King}}/M_{\mathrm{vir}}`$. This ratio is plotted in Fig. Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. and is evidently very close to 1 for the relevant range of concentration parameters.

The average of the $`M/L_V`$ ratios in Table 4, $`M/L_V=1.53\pm 0.18`$, is in excellent agreement with McLaughlin’s $`1.45\pm 0.1`$ value for the core $`M/L_V`$ of Milky Way GCs. Cluster R12 formally has a somewhat lower $`M/L_V`$ than the other three clusters, but we note that this is also the most compact of the clusters and it is possible that the effects of mass segregation and departures from the single-mass King model approximation are more pronounced here. For 6 clusters in M31 with HST imaging, Dubath & Grillmair (1997) obtain a mean $`M/L_V`$ ratio of $`M/L_V=1.2\pm 0.21`$ or $`M/L_V=2.0\pm 0.44`$, based on King model fits and use of the virial theorem, respectively. These estimates nicely bracket our mean value for the M33 clusters. If we had used the velocity dispersions based on the direct fitting instead of those obtained by cross-correlation, the masses (and resulting $`M/L_V`$ ratios) would have increased by $`20`$%, still in good agreement with the Milky Way value.

## 4 Fundamental plane relations

We finally explore how the M33 clusters fit into the “fundamental plane” relations obeyed by Milky Way globular clusters. Figure Structural Parameters and Dynamical Masses for Globular Clusters in M33 111Based on data obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. shows our M33 data compared with data for Milky Way GCs, taken from Pryor & Meylan (1993) and Harris (1996). Core-collapsed Milky Way clusters have been omitted from the plots. Panels (a) and (b) show the bivariate core parameter relations (Djorgovski, 1995), with Milky Way GCs plotted as $`+`$ symbols and M33 clusters as filled circles. As shown by Djorgovski, both of these relations are consistent with GCs being virialized systems with a constant $`M/L`$ ratio. The M33 clusters fit remarkably well onto the the fundamental plane relations obeyed by the Milky Way clusters, consistent with their similar mass-to-light ratios. We note that Barmby, Holland, & Huchra (2002) found slight deviations between the fundamental plane relations for their sample of M31 globular clusters and the Milky Way data, but it is unclear to what extent those differences are due to systematic measurement errors or real physical differences between the GCs in the two galaxies.

Panel (c) shows the cluster binding energies as a function of absolute $`V`$ magnitude, $`M_V`$. McLaughlin (2000) found a tight relation between these two quantities for Milky Way GCs, with an additional weaker dependence on Galactocentric distance $`r_{\mathrm{gc}}`$. He gave the relation as

$$\mathrm{log}(E_b/\mathrm{ergs})=[(39.86\pm 0.40)0.4\mathrm{log}(r_{\mathrm{gc}}/8\mathrm{k}\mathrm{p}\mathrm{c})]+(2.05\pm 0.08)\mathrm{log}(L/L_{}).$$

For the four M33 globular clusters we find a mean value of $`\mathrm{log}(E_bL^{2.05})=\mathrm{\hspace{0.17em}39.81}`$, with a scatter of only 0.13 dex. Ignoring the term in $`r_{\mathrm{gc}}`$, we thus find once again that the M33 globular clusters fall nearly exactly on the same relation as the Milky Way globular clusters.

The similarity to Milky Way GCs in general, and in the mass-to-light ratios in particular for M33 clusters is interesting in the context of the results by Sarajedini et al. (1998) and Chandar et al. (2002), since a younger age should generally imply lower M/L ratios unless there is an excess of low-mass stars. For a Miller-Scalo IMF, for example, population synthesis models by Bruzual & Charlot (priv. comm.) predict an increase in $`M/L_V`$ of about a factor of 2 from 5 Gyrs to 12 Gyrs, and a factor $`5`$ from 2 Gyrs to 12 Gyrs. Indeed, Fischer, Welch, & Mateo (1992) find a much lower M/L ratio ($`M/L_V0.2`$) for the $`2`$ Gyr old cluster NGC 1978 in the LMC than for any of the clusters studied here. However, a younger age for the M33 clusters is still possible if a steeper IMF slope conspires with a young age to produce the same M/L. Of our four GCs, Chandar et al. (2002) list spectroscopic age estimates for U49 (12–17 Gyrs) and H38 (5 Gyrs), but we note that absolute ages based on integrated spectra are still uncertain (e.g. Schiavon et al., 2002). On the other hand, if H38 and U49 are both old clusters ($`10`$ Gyr), then alternative explanations will be required for their excessively red horizontal branch morphologies, the presence of stars above their RGB tips, and their anomalously bright red clump magnitudes (Sarajedini et al., 1998, 2000). Clearly, more work needs to be done before the question of a real, large age spread among M33 halo GCs can be answered with confidence.

## 5 Summary

We have measured velocity dispersions for 4 old “halo” globular clusters in M33 and combined these with structural parameters from WFPC2 images to obtain mass-to-light ratios for the clusters. Our analysis shows that the M33 GCs have $`M/L_V`$ ratios essentially identical to those of GCs in the Milky Way and M31, as well as similar structural parameters (core- and effective radii). The M33 clusters also fall on the same “fundamental plane” and binding energy vs. luminosity relations as Milky Way and M31 clusters (Djorgovski, 1995; Dubath & Grillmair, 1997; McLaughlin, 2000). Unless the stellar IMFs in these clusters have an excess of low-mass stars relative to GCs in the Milky Way, the similarity in M/L ratios and the fundamental plane relations suggest similar high ages. It would be very interesting to measure velocity dispersions for M33 clusters spanning a wide range of ages and look for the expected variations in M/L ratio.

This work was supported by National Science Foundation grant number AST9900732 and by NSF Career Grant No. AST0094048. JPH was supported by the Smithsonian Institution. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. We thank Rupali Chandar for helping to sort out the differences between our photometry and hers, and the Keck staff for their assistance with the observations.


format: text

text:
# On angular momentum transport in convection-dominated accretion flows

## 1 Introduction

A phenomenon of underluminous accreting solar mass black holes in X-ray binaries and supermassive black holes in galactic nuclei has stimulated recent investigations of radiatively inefficient accretion flows (see Narayan, Mahadevan & Quataert 1998 and Narayan 2002 for reviews). Contrary to a standard Shakura-Sunyaev accretion disk model, which successfully explains relatively soft and luminous X-ray sources, models of radiatively inefficient flows are used to explain the significant deficit of radiation observed in some hard X-rays sources. A particular example of such underluminous sources is the Galactic center, Sagittarius $`A^{}`$, which hosts a $`2\times 10^6M_{}`$ black hole. The Galactic center has a luminosity that is well below estimates based on the models of Shakura-Sunyaev accretion disk or spherical Bondi accretion flow (Melia & Falcke 2001). Other examples of underluminous supermassive black holes have been provided by recent Chandra observations of ellipticals NGC 1399, NGC 4472, NGC 4636 and NGC 6166 (Loewenstein et al. 2001; di Matteo et al. 2001).

Numerical simulations of radiatively inefficient accretion flows, made in the framework of viscous hydrodynamical (HD) approach, revealed that low viscosity flows are convectively unstable and convection strongly influences structure of such flows (Igumenshchev, Chen, & Abramowicz 1996; Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle, & Begelman 1999; Igumenshchev, Abramowicz, & Narayan 2000; McKinney & Gammie 2002). An analytical model, which reproduces basic features of the HD simulations, was called convection-dominated accretion flow (CDAF, see Narayan, Igumenshchev, & Abramowicz 2000; Quataert & Gruzinov 2000). CDAFs consist of a hot plasma at about virial temperatures and have a flatten time-averaged radial density profile, $`\rho R^{1/2}`$, where $`\rho `$ is the density and $`R`$ is the radius. In CDAFs, a significant fraction of binding energy, most of which is released in the innermost region of accretion flows, is transported outward by convection motions. HD CDAFs rotate with about Keplerian angular velocity (somewhat below the Keplerian value), and angular momentum is simultaneously transported inward by convection and outward by anomalous viscosity, so the net angular momentum flux is nearly zero (Narayan et al. 2000).

HD models of rotating accretion flows use ad hoc assumption on the nature of angular momentum transport. Specifically, one must assume a non-zero viscosity $`\nu `$ to provide outward angular momentum transport in the flows; without such a transport inward spiraling of matter into the black hole is not possible. An important disadvantage of the HD approach is that functional dependences of $`\nu `$ on other parameters of the problem can not be self-consistently found. Moreover, it is not even clear that a HD viscous diffusion model of angular momentum transport could correctly describe real flows. It is believed that the problem can be self-consistently solved in the frame-work of MHD approach. A small-scale MHD turbulence (with characteristic length-scale of magnetic field $`\mathrm{}<R`$) initiated by the magneto-rotational instability (Balbus & Hawley 1991) and/or large-scale magnetic loops ($`\mathrm{}R`$) could be responsible for angular momentum transport in rotating flows.

The problem of MHD accretion involves non-linear physical effects, requiring numerical simulations. Global MHD numerical simulations of radiatively inefficient accretion flows have been started recently by several research groups, thanks to development of high performance computers and advances in numerical technics. Two groups have used similar designs of numerical experiments, in which calculations have been started with an equilibrium torus orbiting a black hole. Despite of this similarity, the designs are importantly different in initial configuration of magnetic field. One group ( Machida, Matsumoto, & Mineshige 2001) assumed that the torus has only toroidal component of magnetic field at the beginning of simulations. Evolution of the torus after several orbital periods was resulted in development of small-scale MHD turbulence and convection motions; inner parts of the torus formed a turbulent accretion flow. Machida et al. (2001) noted that the accretion flow has radial structure similar to CDAFs found in the low viscosity HD simulations. Another group (Stone & Pringle 2001; Hawley 2000; Hawley, Balbus, & Stone 2001; Hawley & Balbus 2002) assumed poloidal magnetic field inside the tori as an initial condition of their numerical models. Resulting accretion flows were significantly different from the results of Machida et al. (2001). These flows formed a powerful bipolar magnetized outflows, which sandwich a relatively narrow region of equatorial inflow. The bipolar outflows efficiently remove angular momentum from the inflowing material, providing accretion. Convection was not observed in these simulations. The numerical results of these two groups have demonstrated that topology of magnetic field in rotating accretion flows is an important factor determined the flow structure, and mechanisms and efficiency of angular momentum transport.

Study of spherical (non-rotating) MHD accretion flows provided further insight on the role of convection in radiatively inefficient accretion flows. It was shown with help of 3D MHD numerical simulations (Igumenshchev & Narayan 2002) that magnetized Bondi flows experience a transition to convection-dominated regime of accretion. This regime (convection-dominated Bondi flow or CDBF) has no HD counterparts among non-rotating accretion flows, but it is similar in many aspects to the rotating HD CDAFs. Unless different physical mechanisms are responsible for development of convection in HD CDAF and CDBF, simulations showed that convection itself governs the flow structure independently on details of those mechanisms. This explains a formal similarity of analytic self-similar solutions described HD CDAFs and CDBFs (see Narayan et al. 2001; Igumenshchev & Narayan 2002).

In this Letter I will discuss properties of CDAF with magnetic field, and will concentrate on the problem of angular momentum transport. This study is motivated by recent works by Balbus and collaborators (Hawley et al. 2001; Balbus & Hawley 2002), who argued that there is a fundamental difference in dynamics of MHD and viscous HD flows. Because of this difference and consequent difference in mechanisms of angular momentum transport, they concluded that CDAFs can not be formed with magnetized plasma (for discussions, see Narayan et al. 2002). I will show, based on analytic arguments and results of previous numerical simulations, that CDAFs can exist independently of the changes of details of angular momentum transport, provided that the changes are due to magnetic field.

## 2 Self-similar solutions

Analytical theory of CDAF is based on a self-similar solution of a simplified set of equations described radiatively inefficient accretion flows. I briefly reproduce here main results of the theory following to Narayan et al. (2000) and Igumenshchev & Narayan (2001). Consider stationary height-integrated equations described an accretion flow onto the black hole of mass $`M`$. In absence of mass outflows, the continuity equation reads

$$\dot{M}=4\pi RH\rho v_R,$$
$`(2.1)`$
where $`\dot{M}`$ is the accretion rate, $`H=(c_s/v_K)R`$ is the scale height, $`c_s=\sqrt{P/\rho }`$ is the sound speed, $`P`$ is the pressure, $`v_R`$ is the accretion velocity, $`v_R<0`$, $`v_K=\sqrt{GM/R}`$ is the Keplerian velocity and $`G`$ is the gravitational constant.

In the radial momentum equation, I ignore the convection term $`v_Rdv_R/dR`$, assuming that $`v_R^2v_K^2`$. In addition, I ignore the turbulent pressure and consider the total pressure $`P`$ consisted of the gas pressure $`P_g=R_g\rho T`$ and magnetic pressure $`P_m=P_g/\beta `$, $`P=P_g+P_m=(\beta +1)P_g/\beta `$, where $`T`$ is the temperature, $`R_g`$ is the gas constant, and $`\beta `$ is a parameter. Then the equation takes the form

$$\mathrm{\Omega }_K^2R\mathrm{\Omega }^2R=\frac{1}{\rho }\frac{dP}{dR},$$
$`(2.2)`$
where $`\mathrm{\Omega }`$ is the angular velocity and $`\mathrm{\Omega }_K=v_K/R`$ is the Keplerian angular velocity. Equation (2.2) describes a balance between gravitational, centrifugal and pressure gradient forces.

The angular momentum equation can be written in the form of balance of advection and diffusion transport terms,

$$\frac{1}{R}\frac{d}{dR}\left(R\mathrm{\Sigma }v_R\right)=\frac{1}{R}\frac{d}{dR}\left(\mathrm{\Sigma }\nu _{out}R^3\frac{d\mathrm{\Omega }}{dR}\right)+\frac{1}{R}\frac{d}{dR}\left(\mathrm{\Sigma }\nu _{in}R\frac{d}{dR}\right),$$
$`(2.3)`$
where $`\mathrm{\Sigma }=2\rho H`$ is the surface density and $`=\mathrm{\Omega }R^2`$ is the specific angular momentum. The term on the left hand side of (2.3) describes inward advection of angular momentum, whereas the two terms on the right hand side correspond to two different mechanisms of diffusion angular momentum transport. The first term on the right hand side of (2.3), involved the viscosity coefficient $`\nu _{out}`$, describes the standard viscous transport, when the corresponded angular momentum flux is proportional to the gradient of angular velocity. In stationary accretion flows, $`\mathrm{\Omega }`$ is always decreasing function of radius, so the term describes outward transport of angular momentum. In the second term on the right hand side of (2.3), the angular momentum flux is proportional to the gradient of angular momentum with the viscosity coefficient $`\nu _{in}`$. This term describes a tendency of axisymmetric turbulence to uniform distribution of angular momentum in accretion flows (e.g. Ryu & Goodman 1992). In flows with nearly Keplerian rotation, $`\mathrm{\Omega }R^{3/2}`$, the latter term represents inward transport of angular momentum.

The energy equation is

$$\rho v_RT\frac{ds}{dR}\rho v_R\left(\frac{1}{\gamma 1}\frac{dc_s^2}{dR}\frac{c_s^2}{\rho }\frac{d\rho }{dR}\right)=\frac{1}{R^2}\frac{d}{dR}\left(R^2F\right)+Q,$$
$`(2.4)`$
where $`s`$ is the specific entropy, $`\gamma `$ is the adiabatic index, $`F`$ is the outward energy flux due to convection, and $`Q`$ is the energy dissipation rate per unit volume. Equation (2.4) states that the locally released internal energy (the term $`Q`$) can be advected inward (the term involved $`ds/dR`$) and/or transported outward by convection (the term involved $`F`$).

In the considered theory, simple parametric expressions are used to represent the viscosity coefficients $`\nu _{out}`$ and $`\nu _{in}`$,

$$\nu _{out}=\alpha _{out}c_sH,$$
$`(2.5)`$
$$\nu _{in}=\alpha _{in}c_sH,$$
$`(2.6)`$
the convective flux,

$$F=\alpha _cc_sR\rho T\frac{ds}{dR},$$
$`(2.7)`$
and the dissipation rate per unit volume,

$$Q=\alpha _d\frac{v_R}{R}\rho \frac{GM}{R},$$
$`(2.8)`$
where $`\alpha _{out}`$, $`\alpha _{in}`$, $`\alpha _c`$ and $`\alpha _d`$ are positive dimensionless parameters.

Equations (2.1)-(2.4) with (2.5)-(2.8) have two types of self-similar solutions, which correspond to two values of the self-similar index $`a=3/2`$ and $`1/2`$. The solutions can be written as follows,

$$\rho (R)=\rho _0R^a,$$
$`(2.9)`$
$$v_R(R)=v_0v_K\left(\frac{R_G}{R}\right)^{3/2a}R^{a2},$$
$`(2.10)`$
$$c_s(R)=c_0v_KR^{1/2},$$
$`(2.11)`$
$$\mathrm{\Omega }=\mathrm{\Omega }_0\mathrm{\Omega }_KR^{3/2},$$
$`(2.12)`$
where $`R_G=2GM/c^2`$ is the gravitational radius and $`c`$ in the speed of light. The dimensionless coefficients $`v_0`$, $`c_0`$, and $`\mathrm{\Omega }_0`$ in (2.10)-(2.12) can be expressed algebraically through the parameters $`\alpha _{out}`$, $`\alpha _{in}`$, $`\alpha _c`$, $`\alpha _d`$, $`\beta `$ and $`\gamma `$, substituting (2.9)-(2.12) into (2.2)-(2.4). Note, that not all values of the parameters lead to a consistent solution. The coefficient $`\rho _0`$ in (2.9) is proportional to the accretion rate $`\dot{M}`$ and scales out of the problem. The solution with $`a=3/2`$ corresponds to advection-dominated accretion flow (e.g. Narayan & Yi 1994; Abramowicz et al. 1995), whereas the solution with $`a=1/2`$ represents CDAF (Narayan et al. 2000; Quataert & Gruzinov 2000). Convection-dominated solution with $`a=1/2`$ and $`\mathrm{\Omega }_0=0`$ describes CDBF (Igumenshchev & Narayan 2002).

## 3 Angular momentum transport

Consider the CDAF self-similar solution. In this case one can neglect the advection terms (the terms involved $`v_R`$) in (2.3)-(2.4). Taking into account (2.5) and (2.6), equation (2.3) can be re-written in the form,

$$\frac{d\mathrm{ln}\mathrm{\Omega }}{d\mathrm{ln}R}=\frac{2\alpha _{in}}{\alpha _{out}+\alpha _{in}}.$$
$`(3.1)`$
HD CDAFs only exists in the form of flows with nearly Keplerian rotation, $`\mathrm{\Omega }\mathrm{\Omega }_K`$. From this fact and from equation (3.1), one obtains the following relation between parameters $`\alpha _{in}`$ and $`\alpha _{out}`$,

$$\alpha _{out}=\frac{1}{3}\alpha _{in}.$$
$`(3.2)`$
The relation (3.2) significantly limits properties of convection turbulence in rotating CDAFs: an intensity of axisymmetric convection motions, which transport angular momentum inward, must be supported at a some specific level, given by (3.2), to compensate viscous outward spread of angular momentum. Could this be truth for real flows? HD simulations of radiatively inefficient accretion flows revealed that only CDAF-type solutions exist in the low viscosity limit. These solutions have nearly Keplerian rotation, $`\mathrm{\Omega }\mathrm{\Omega }_K`$, and convection in these flows is almost axisymmetric due to ability of the Keplerian shear motion to uniform convective perturbations in $`\varphi `$-direction (azimuthal direction). It was found with help of 2D and 3D HD simulations (Igumenshchev & Abramowicz 2000; Igumenshchev et al. 2000) that the time-averaged $`R\varphi `$-component of Reynolds stress tensor has negative sign in CDAFs, which confirms that axisymmetric convection transports angular momentum inward. Moreover, strengths of inward and outward angular momentum fluxes were estimated to be close to each other in the numerical models. Thus, results of numerical simulations of HD CDAFs are in a good qualitative and quantitative agreement with results of the self-similar theory (for more detailed comparison of the theory and simulations, see Igumenshchev & Abramowicz 2000).

Magnetic fields can affect properties of CDAFs found in HD simulations. As it was discussed in introduction, there is an important difference between HD and MHD models of radiatively inefficient accretion flow that MHD models have a non-rotating counterpart of CDAFs, whereas HD models do not have it. So, in MHD flows there is no longer the requirement of nearly Keplerian rotation, and in principle, the accretion flow can consist of rotating and non-rotating parts with transition regions between them. Solutions described such transition regions can not be self-similar and they are out of scope of this Letter. However, one can better understand properties of MHD CDAFs considering flows of extreme rotation rates.

Consider slowly rotating ($`\mathrm{\Omega }\mathrm{\Omega }_K`$) MHD CDAFs, in which rotation has negligible dynamical effect. Radial structure of such flows can be described by the CDBF self-similar solution. The law of rotation in the flows can be found with help of (3.1), if one knows the parameters $`\alpha _{out}`$ and $`\alpha _{in}`$. Relative importance of the parameters can be estimated as follows. In slowly rotating flows, a weak rotation shear is not able to uniform non-axisymmetric perturbations in $`\varphi `$-direction, as it is in the case of fast rotating HD CDAFs, and convection motions are unlikely to be axisymmetric. Thus, both non-axisymmetric convection and small scale magnetic turbulence transport angular momentum outward much more efficiently, than inward, and one gets $`\alpha _{out}\alpha _{in}`$. In this limit, using (3.1), one finds that the flow must rotate as a rigid-body, $`\mathrm{\Omega }(R)=\mathrm{\Omega }_B`$, where $`\mathrm{\Omega }_B`$ is the angular velocity at the outer boundary (according to our assumption, $`\mathrm{\Omega }_B\mathrm{\Omega }_K`$). Relative dynamical importance of rotation in such slowly rotating CDAFs is decreased with decreasing of the radius.

Next, consider MHD CDAFs which rotates with about Keplerian angular velocity at the outer boundary, $`\mathrm{\Omega }_B\mathrm{\Omega }_K(R_B)`$. Two possible structures of the flow can be proposed in this case, depending on the mechanism of angular momentum transport. The first structure is identical to the structure of HD CDAF, in which the flow rotates with nearly Keplerian angular velocity everywhere down to the inner boundary, $`\mathrm{\Omega }\mathrm{\Omega }_K`$. This requires that convection turbulence in MHD CDAFs has the same properties as in the case of HD CDAFs: axisymmetric convection transports angular momentum inward and almost compensates outward angular momentum flux due to small scale MHD turbulence. In such rotating MHD CDAFs, $`\alpha _{in}`$ and $`\alpha _{out}`$ satisfy to (3.2). The second structure is formed if the outward angular momentum flux exceeds the inward one. In this case one has (see eq. \[3.1\])

$$\alpha _{out}>\frac{1}{3}\alpha _{in}\mathrm{and}\frac{d\mathrm{ln}\mathrm{\Omega }}{d\mathrm{ln}R}>\frac{d\mathrm{ln}\mathrm{\Omega }_K}{d\mathrm{ln}R}.$$
$`(3.3)`$
The excessive flux results in formation of a region inside $`R_B`$, in which the flow slows down rotation with respect to Keplerian one, as it follows from (3.3). In the limit $`RR_B`$, when the rotation is reduced significantly, one obtains the slowly rotating MHD CDAF considered above. Thus, the second possible structure of MHD CDAFs can consist of two parts: the outer transition part, where the flow rotation is slowed down from the boundary level at $`R_B`$ till significantly sub-Keplerian level at $`RR_B`$, and the inner part, where the flow rotates slowly, almost as a rigid body.

## 4 Conclusion

The convection-dominated accretion flow (CDAF) model consistently represents radiatively inefficient accretion flows into black holes in the frame-work of viscous hydrodynamical (HD) approach. In HD CDAFs, rotation is nearly Keplerian, supported by simultaneous and balanced transport of angular momentum inward by convection and outward by viscosity. It is not yet clear how the MHD approach could change the CDAF model. Results of direct numerical simulations of magnetized radiatively inefficient accretion flows indicate that solutions strongly depend on assumed topology of magnetic field in the flow (Machida et al. 2001; Hawley & Balbus 2002). At present, however, absence of detailed MHD simulations makes it difficult to describe actual structure of magnetized radiatively inefficient accretion flows.

In this Letter I have shown that the existence and properties of CDAF do not crucially depend on direction of convection angular momentum transport in magnetized plasma, contrary to what was thought previously based on results of HD simulations. MHD CDAF can exist even if convection presumably transports angular momentum outward, rather than inward. In this case of dominant outward transport the flow at small radii becomes similar to non-rotating convection-dominated Bondi flow (CDBF). I have considered two possible structures of MHD CDAF, which are determined by the relative strengths of inward and outward angular momentum transport. One possibility is that the radial structures and rotation profiles of MHD and HD CDAFs are similar, supposing that convection in magnetized accretion flows transports angular momentum inward with the same efficiency as it does in viscous HD flows. If, however, the inward transport is less efficient in MHD flows, than the outward one, then the radial structure of CDAFs can be different. In the latter case the flow consists of two parts. In the outer part, the angular velocity is reduced from the outer boundary value (supposed to be close to Keplerian one) till significantly sub-Keplerian value. In the inner part, the flow rotates very slowly, so its radial structure is represented by the self-similar CDBF solution. A choice between these two possible structures of MHD CDAFs can be made on the base of high-resolution 3D numerical simulations.

###### Acknowledgements.

Author thanks Ramesh Narayan and Marek Abramowicz for helpful discussions and support. The study was supported by the U.S. Department of Energy (DOE) Office of Inertial Confinement Fusion under Cooperative Agreement No. DE-FC03-92SF19460, the University of Rochester, the New York State Energy Research and Development Authority, and RFBR grant 00-02-16135. The support of the DOE does not constitute an endorsement by the DOE of the views expressed in this article.


text:
# 1 Introduction

## 1 Introduction

The long term monitoring of microquasar GRS 1915+105 has shown broad correlation between the X-ray and non-thermal radio emission . The radio emission in GRS 1915+105 can be classified into three classes; (i) the relativistic superluminal radio jets of flux density $``$ 1 Jy with decay time-scales of several days (, ), (ii) the baby jets of 20 $``$ 40 min durations with flux density of 20 $``$ 200 mJy both in infrared (IR) and radio (, ), and (iii) the plateau state with persistent radio emission of 20 $``$ 100 mJy for extended durations . In the case of superluminal jets, the radio emission has steep spectra and are observed at large distances (400 $``$ 5000 AU) from the accretion disk , . The radio emission during other two classes has flat spectra and they occur close to the accretion disk (within a few tens of AU). The multi-wavelength studies have indicated strong disk-jet connection for class 2(, , ). In this paper, we present our radio observations at 1.28 and 15 GHz during June 18 to July 3, 2001 (just prior to when source went to the plateau state on July 3). These observations include faint radio emission, isolated baby radio flares riding above smooth as well as slowly rising/decaying radio emission (broad flares), and large radio flare ($``$ 200 mJy on June 20).

## 2 Observations and discussion

The 1.28 GHz radio observations were carried out with a bandwidth of 16 MHz using the GMRT at Pune, India  on 2001 June 18, 22, 23, 27, 28, 29, 30 and July 1. The flux density scale is set by observing primary calibrator 3C286 or 3C48. The integration time of 32 s was chosen. The data recorded from GMRT has been converted to FITS and was analysed using Astronomical Image Processing System (AIPS). The details of observations/analysis is given elsewhere . The radio observations at 15 GHz are from the Ryle telescope at Cambridge. Details of the instrument and analysis can be found elsewhere .

The RXTE/ASM flux and X-ray hardness ratio (5-12 keV/1.5-5 keV) are shown in the top panel of Figure 1. The radio lightcurves at 1.28 and 15 GHz produced at 5 min interval are shown in the bottom panel of Figure 1. GRS 1915+105 exhibited significant radio emission on all days, except on June 27, when the source showed much weaker radio emission ($``$ 5 mJy). Figure 1 brings out following important points:

1. The radio fluxes at 1.28 and 15 GHz are consistent and suggest a flat radio spectrum during these observations except during the large radio flare on June 20-24. The flux and decay time of this flare suggest that it belongs to the class of relativistic jets. The higher flux at 1.28 GHz is consistent with steep spectrum during the decay of this flare.
2. The X-ray hardness ratio changed from $``$ 0.7 to $``$ 0.85 around June 30 prior to when source went to the plateau state on July 3. The RXTE/ASM flux becomes stable at higher hardness ratio. The X-ray hardness ratio suggests that the source remains in the low-hard state during this period.

The radio light curve observed on June 28 is shown in Figure 2. The flux density was $``$ 5 mJy at about UTC 17 hour and reached gradually to 50 mJy at UTC 24 hour. When the observations resumed on 29 (not shown here), the source was ”caught” at 70 mJy at UTC 16 hour, and the flux started decaying slowly to the value of 10 mJy at UTC 24 hour. This is probably the first detail observation of the precursor flares to the plateau state with rise and decay times of over six hours. It is also important to note that these broad flares and the change in the hardness ratio around June 30 set the stage for the plateau state.

These broad flares consist of mini (baby) radio flares riding above rising/decaying radio emission. In contrast, the radio lightcurve observed on June 30 shows mini (baby) flares riding above almost smooth radio emission (Figure 3).

These isolated baby flares are modeled as adiabatically expanding synchrotron clouds ejected from the accretion disk . One such flare along with model fit data is shown in Figure 4. This model provides estimate of the spectrum index from single frequency observations and successfully explains the observed delay times between different frequencies. The radio emission during these observations is consistent with flat radio spectrum which is in agreement with observed flux at 15 GHz.

## Acknowledgments

We thank the staff of the GMRT that made these observations possible. GMRT is run by the National Center of Astrophysics of the Tata Institute of Fundamental Research. We also thank RXTE/ASM team for making their data public.


text:
# UBVRI photometric sequences for symbiotic stars

## 1. Introduction

Availability of suitable UBVRI photometric sequences around symbiotic stars is an essential step in promoting a large-scale documentation effort of present time multi-band lightcurves as well as reconstruction of the historical behavior.

To serve for measuring archival plates, the sequence stars have to be well isolated from neighbors to avoid blending on short focal length patrol plates, have to be grouped close to the symbiotic variable so as to fall within the same eyepiece field of view, and must extend over a wide range in brightness to easily cover both faint states in quiescence (like eclipses) as well as bright outbursts.

To assist present time CCD photometry (carried out mostly by amateurs with telescopes of 1-2 m focal lengths), the sequences must also cover a wide range in colors (to allow simultaneous derivation of transformation coefficients), be compact enough to easily fit onto a small CCD ($`<`$5 arcmin in diameter) and be accurately placed on the UBVRI system.

The UBVRI sequences on the Johnson-Cousins system that we have calibrated around more than 80 symbiotic stars follow the above prescriptions and are discussed in this paper. Their use on Asiago archive plates is described by Jurdana-Sepic and Munari (2002, and references therein).

## 2. The sequences

Figure 1 and Table 1 show a typical UBVRI photometric sequence among the 82 we have calibrated. The increasing errors at the brightest magnitudes are due to saturation effects, and those at the faintest magnitudes are dominated by the photon statistics.

The symbiotic stars for which we have already obtained the multi-epoch UBVRI observations suitable for the calibration of the photometric sequences are listed in Table 2. Of the 82 listed objects, 40 have their sequences published by Henden and Munari (2000, 2001); the remaining 42 will follow shortly. The program stars have been mainly selected from the catalogues of symbiotic stars of Allen (1984) and Belczyński et al. (2000)

All observations were made with the 1.0 m Ritchey-Chrétien telescope of the U.S. Naval Observatory, Flagstaff Station. A Tektronix/SITe 1024$`\times `$1024 thinned, backside-illuminated CCD was used, along with proper Johnson UBV and Kron-Cousins RI filters. The telescope scale is 0.67 arcsec/pix and typical seeing was $``$2 arcsec. The color transformation equations for this instrumental combination are particularly good, with slopes deviating from unity by only a few hundredths of a magnitude (see Eq.(1)-(5) in Henden and Munari 2000). Only data collected under photometric conditions (transformation errors under 0.02 mag) have been used, and the calibration has been made against the Landolt equatorial fields (Landolt 1983, 1992).

Each photometric sequence is composed of 10 to 15 stars, which have been selected and ordered on the basis of their $`B`$ band magnitude. The latter is reproducible by most filter-equipped CCD cameras, it is the closest one to the $`m_{pg}`$ band of historical observations, and the $`B`$ band is particularly well suited to investigate the variability of symbiotic stars. Stars of varying color are included in the sequence to allow simultaneous determination of transformation coefficients. The comparison sequences are tighter around the usual brightness of each symbiotic star and become looser away from it.

The stars included in the comparison sequences have been checked on at least three nights for variability. We cannot obviously rule out beyond doubt that some of them are indeed variable, but the agreement at a few milli-mag level of their magnitude measured on different nights over an interval of months gives some confidence in their use.

## References

Allen, D.A. 1984, Proc. A.S.A., 5, 369

Belczyński K., Mikolajewska J., Munari U., Ivison R.J. & Friedjung M. 2000, A&AS, 146, 407

Jurdana-Sepic, R., Munari, U., 2002, this volume

Henden A., Munari U. 2000, A&AS143, 343

Henden A., Munari U. 2000, A&A372, 145

Landolt, A. U., 1983, AJ88,439.

Landolt, A. U., 1992, AJ104,340.


text:
# Low-frequency characterization of quantum tunneling in flux qubits

## I Introduction

Josephson-junction flux qubits are known to be candidates for solid-state quantum computing circuits.Makhlin  This qubit variety has good tolerance to external noise, especially to dangerous background-charge fluctuations.Blatter  A flux qubit is a superconducting loop, the two lowest-energy states of which differ in the direction of circulating persistent current. For many flux qubits, these two states become degenerate when the external flux $`\mathrm{\Phi }_\mathrm{x}`$ threading the loop equals $`\mathrm{\Phi }_0/2`$ ($`\mathrm{\Phi }_0=h/2e`$ is the flux quantum), and quantum tunneling between them becomes possible. Moving $`\mathrm{\Phi }_\mathrm{x}`$ away from $`\mathrm{\Phi }_0/2`$ lifts the degeneracy and applies a bias between the two states. When the biasing energy exceeds the tunneling amplitude $`\mathrm{\Delta }`$ the tunneling stops, but the relative phase between the two states will still evolve in time. This, together with coherent tunneling, provides single-bit quantum gate operations. To have a universal set of gates, necessary for quantum computing, one needs to be able to couple two qubits. The methods of coupling two flux qubits and performing gate operations are beyond the present scope. Instead, we propose a method to characterize the quantum behavior of a flux qubit by coupling it to a tank circuit. The discussion will be quite general and can be applied to different types of flux qubit such as rf-SQUID,rfsquid  three-Josephson-junction (3JJ),Mooij ; Orlando  multi-terminal,amin02  etc. We will use the example of the 3JJ qubit, where quantum superposition of the macroscopic current states has been observed.Mooij

Due to the loop self-inductance, the total qubit flux $`\mathrm{\Phi }`$ may differ from $`\mathrm{\Phi }_\mathrm{x}`$, depending on the direction of the persistent current. Figure 1 shows the $`\mathrm{\Phi }`$$`\mathrm{\Phi }_\mathrm{x}`$ curve for a typical flux qubit. The solid lines correspond to classical behavior. Near the degeneracy point, the diagram is hysteretic, a signature of the qubit’s bistability. This has been observed for the 3JJ in Refs. Tanaka, ; APL7, . In the quantum regime, tunneling between the states at degeneracy may eliminate the hysteresis (dashed line in Fig. 1). This phenomenon will be discussed in detail below.

In general, one can plot the classical (local) minimum energies of a flux qubit as in Fig. 2a. The left (right) branch then corresponds to (counter-)clockwise flow of the spontaneous current. The hysteresis is also evident from this diagram. In the quantum regime, there will be discrete local states in each of the qubit’s bistable potential wells. From now on we denote the lowest-lying such states as $`\mathrm{\Psi }^\mathrm{l}`$ and $`\mathrm{\Psi }^\mathrm{r}`$, corresponding to “left” and “right” directions of the persistent current respectively. At $`\mathrm{\Phi }_\mathrm{x}=\mathrm{\Phi }_0/2`$, resonant tunneling will render the lowest eigenstates of the full Hamiltonian as superpositions $`(\mathrm{\Psi }^\mathrm{l}\pm \mathrm{\Psi }^\mathrm{r})/\sqrt{2}`$. A small splitting equal to $`2\mathrm{\Delta }`$ will appear between their energies (Fig. 2b). Starting with the qubit in its ground state (lower band in Fig. 2b), adiabatically changing $`\mathrm{\Phi }_\mathrm{x}`$ will keep it in the ground state. This means that by passing through the degeneracy point, the qubit will continuously transform from $`\mathrm{\Psi }^\mathrm{l}`$ to $`\mathrm{\Psi }^\mathrm{r}`$. This pure quantum behavior is shown by the dashed line in Fig. 1. On the other hand, if $`\mathrm{\Phi }_\mathrm{x}`$ changes rapidly, there is a considerable probability to excite the qubit and therefore continue on the same classical branch (left or right). This so-called Landau–Zener effect can be used to distinguish the classical from the quantum energy curves.

The curvature of the energy profile is related to the qubit’s effective inductance and is therefore important for measurement. Figure 3 displays the second derivative of the curves in Fig. 2. In the classical regime (Fig. 3a), the hysteretic behavior is the same as for the energy. On the other hand, in the quantum regime the hysteresis is replaced by a sharp spike, due to the level anticrossing. The appearance of this spike can be ascribed to enhanced susceptibility of the system due to tunneling. Its size and width can provide information about $`\mathrm{\Delta }`$.

A simple experimental implementation is to inductively couple the qubit to an $`LC`$ tank circuit with known inductance $`L_\mathrm{T}`$, capacitance $`C_\mathrm{T}`$, and quality $`Q`$ through a mutual inductance $`M`$ (Fig. 4). The resonant characteristics of the tank circuit (frequency, phase shift, etc.) will then be sensitive to the qubit inductance and therefore to its energy curvature. In particular, the spike in Fig. 3b appears as sharp dips in both phase shift and tank voltage as a function of $`\mathrm{\Phi }_\mathrm{x}`$ (see Section III).

This method, known as impedance measurement technique (IMT), has been used for current–phase measurements of Josephson junctions. It originates from the pioneering work of Rifkin and Deaver,Rifkin-Deaver  and is analyzed in detail in Ref. IMT-2, . IMT has also successfully been applied to a 3JJ qubit in the classical regime,APL7  and the hysteretic dependence of the ground-state energy on $`\mathrm{\Phi }_\mathrm{x}`$ (cf. Fig. 2a) was observed as predicted in Ref. Orlando, . The method has also been used for the investigation of quantum transitions in an rf-SQUID (Ref. Clark,  and references therein).

First of all, in Section II we calculate the two qubit energies in more detail than in Ref. Orlando, . In Section III, we study the qubit’s interaction with a high-quality resonant tank, showing that low-frequency IMT yields useful information about the qubit’s quantum behavior. Finally, in Section IV, the effect of noise is considered.

## II Quantum dynamics of the 3JJ qubit

The 3JJ qubitOrlando  consists of three Josephson junctions in a loop with very small inductance $`L`$, typically in the pH range. This insures effective decoupling from the environment. Two junctions have equal critical current $`I_\mathrm{c}`$ and (effective) capacitance $`C`$, while those of the third junction are slightly smaller: $`\alpha I_\mathrm{c}`$ and $`\alpha C`$, with $`0.5<\alpha <1`$. If the Josephson energy $`E_\mathrm{J}=I_\mathrm{c}\mathrm{\Phi }_0/2\pi `$ is much larger than the Coulomb energy $`E_C=e^2/2C`$, the Josephson phase is well defined. Near $`\mathrm{\Phi }_\mathrm{x}=\mathrm{\Phi }_0/2`$, this system has two low-lying quantum states.Orlando ; Mooij1  The energy splitting between them in the presence of a small flux bias has been given in Ref. Mooij, , but only for a particular choice of, e.g., $`\alpha `$ and $`gE_\mathrm{J}/E_C`$. In this section we derive the splitting with its explicit dependence on the qubit parameters. The energy levels are derived from the Hamiltonian \[see Eq. (12) in Ref. Orlando, \]

$$H_0=\frac{P_\phi ^2}{2M_\phi }+\frac{P_\theta ^2}{2M_\theta }+U(f_\mathrm{x},\phi ,\theta ),$$
(1)
where $`\phi =(\phi _1+\phi _2)/2`$, $`\theta =(\phi _1\phi _2)/2`$ with $`\phi _{1,2}`$ the phase differences across the two identical junctions, $`P_\phi =i\mathrm{}_\phi `$, $`P_\theta =i\mathrm{}_\theta `$, $`M_\phi =(\mathrm{\Phi }_0/2\pi )^22C`$, $`M_\theta =(1+2\alpha )M_\phi `$, and

$$U(f_\mathrm{x},\phi ,\theta )=E_\mathrm{J}[\alpha 2\mathrm{cos}\phi \mathrm{cos}\theta +\alpha \mathrm{cos}(2\pi f_\mathrm{x}+2\theta )].$$
(2)
In contrast to Ref. Orlando, , we define the flux bias $`f_\mathrm{x}=\mathrm{\Phi }_\mathrm{x}/\mathrm{\Phi }_0\frac{1}{2}`$ as a small deviation from degeneracy.

Since the qubit is assumed to have small $`L`$ and $`I_\mathrm{c}`$ (typically $`L10`$ pH, $`I_\mathrm{c}100`$ nA), the shielding factor $`LI_\mathrm{c}/\mathrm{\Phi }_00.001`$. Hence, in (1) we have neglected the shielding current, considering $`\mathrm{\Phi }`$ as an external flux.

At $`f_\mathrm{x}=0`$, the potential (2) has two minima at $`\phi =0`$, $`\theta =\pm \theta _{}`$, with $`\mathrm{cos}\theta _{}=1/2\alpha `$ ($`\theta _{}>0`$). Tunneling lifts their degeneracy, leading to energy levels $`E_\pm =\epsilon _0\pm \mathrm{\Delta }`$. To find the levels for $`|f_\mathrm{x}|1`$ we expand Eq. (2) near its minima, retaining linear terms in $`f_\mathrm{x}`$ and quadratic terms in $`\varphi ,\theta `$. Define $`\theta _{}^{\mathrm{r}/\mathrm{l}}`$ as the minima, shifted due to $`f_\mathrm{x}`$:

$$\theta _{}^{\mathrm{r}/\mathrm{l}}=\pm \theta _{}+2\pi f_\mathrm{x}\frac{12\alpha ^2}{4\alpha ^21};$$
(3)
that is, the upper (lower) sign refers to the right (left) well. The potential energy then reads

$`{\displaystyle \frac{U}{E_\mathrm{J}}}`$ $`={\displaystyle \frac{1}{2\alpha }}f_\mathrm{x}{\displaystyle \frac{\pi }{\alpha }}\sqrt{4\alpha ^21}+{\displaystyle \frac{\phi ^2}{2\alpha }}\left(1\pm 2\pi f_\mathrm{x}{\displaystyle \frac{2\alpha ^21}{\sqrt{4\alpha ^21}}}\right)`$
$`+(\theta \theta _{}^{\mathrm{r}/\mathrm{l}})^2\left(2\alpha {\displaystyle \frac{1}{2\alpha }}\pm f_\mathrm{x}{\displaystyle \frac{\pi }{\alpha }}{\displaystyle \frac{2\alpha ^2+1}{\sqrt{4\alpha ^21}}}\right).`$ (4)
Near degeneracy, the eigensolutions of $`H_0\mathrm{\Psi }_\pm =E_\pm \mathrm{\Psi }_\pm `$ can be written as superpositions $`\mathrm{\Psi }_\pm =a_\pm \mathrm{\Psi }^\mathrm{l}+b_\pm \mathrm{\Psi }^\mathrm{r}`$, yielding the well-known eigenenergies $`E_\pm =(\epsilon ^\mathrm{l}+\epsilon ^\mathrm{r})/2\pm \sqrt{(\epsilon ^\mathrm{l}\epsilon ^\mathrm{r})^2/4+\mathrm{\Delta }^2}`$, with $`\epsilon ^{\mathrm{r}/\mathrm{l}}=\mathrm{\Psi }^{\mathrm{r}/\mathrm{l}}|H_0|\mathrm{\Psi }^{\mathrm{r}/\mathrm{l}}`$. The matrix element $`\mathrm{\Delta }`$ cannot accurately be found in terms of $`\mathrm{\Psi }^{\mathrm{r}/\mathrm{l}}`$. In what follows it is assumed constant,

$$\begin{array}{c}\mathrm{\Delta }=2E_\mathrm{J}\sqrt{\frac{2\alpha 1}{\alpha g}}\hfill \\ \hfill \times \mathrm{exp}\left[\sqrt{\frac{g(2\alpha +1)}{\alpha }}(\mathrm{arccos}\frac{1}{2\alpha }\sqrt{4\alpha ^21})\right],\end{array}$$
(5)
neglecting its dependence on $`f_\mathrm{x}`$.

To find the dependence of $`E_\pm `$ on $`f_\mathrm{x}`$, we take $`\mathrm{\Psi }^{\mathrm{r}/\mathrm{l}}`$ to be oscillator ground states in their respective wells:

$$\begin{array}{cc}\hfill \mathrm{\Psi }^{\mathrm{r}/\mathrm{l}}& =\frac{1}{\sqrt{\mathrm{}\pi }}\left(M_\phi \omega _\phi ^{\mathrm{r}/\mathrm{l}}M_\theta \omega _\theta ^{\mathrm{r}/\mathrm{l}}\right)^{1/4}\hfill \\ & \times \mathrm{exp}\left(\frac{M_\phi \omega _\phi ^{\mathrm{r}/\mathrm{l}}}{2\mathrm{}}\phi ^2\frac{M_\theta \omega _\theta ^{\mathrm{r}/\mathrm{l}}}{2\mathrm{}}(\theta \theta _{}^{\mathrm{r}/\mathrm{l}})^2\right),\hfill \end{array}$$
(6)
corresponding to

$$\epsilon ^{\mathrm{r}/\mathrm{l}}=E_\mathrm{J}\left(\frac{1}{2\alpha }f_\mathrm{x}\frac{\pi }{\alpha }\sqrt{4\alpha ^21}\right)+\frac{\mathrm{}\omega _\phi ^{\mathrm{r}/\mathrm{l}}}{2}+\frac{\mathrm{}\omega _\theta ^{\mathrm{r}/\mathrm{l}}}{2},$$
(7)
where

$`\mathrm{}\omega _\phi ^{\mathrm{r}/\mathrm{l}}`$ $`=E_\mathrm{J}\sqrt{{\displaystyle \frac{4}{\alpha g}}}\left(1\pm \pi f_\mathrm{x}{\displaystyle \frac{2\alpha ^21}{\sqrt{4\alpha ^21}}}\right),`$ (8)
$`\mathrm{}\omega _\theta ^{\mathrm{r}/\mathrm{l}}`$ $`=E_\mathrm{J}\sqrt{{\displaystyle \frac{4(2\alpha 1)}{\alpha g}}}\left(1\pm \pi f_\mathrm{x}{\displaystyle \frac{2\alpha ^2+1}{\left(4\alpha ^21\right)^{3/2}}}\right).`$ (9)
Combining the above, one finds the eigenenergies

$$E_\pm =\epsilon _0\pm \sqrt{E_\mathrm{J}^2f_\mathrm{x}^2\lambda ^2(\alpha )+\mathrm{\Delta }^2},$$
(10)
where

$$\epsilon _0=E_\mathrm{J}\left(\frac{1}{2\alpha }+\frac{1+\sqrt{2\alpha 1}}{\sqrt{\alpha g}}\right),$$
(11)
$$\begin{array}{cc}\hfill \frac{\alpha }{\pi }\lambda (\alpha )& =\sqrt{\frac{\alpha }{g}}\left(\frac{2\alpha ^21}{\sqrt{4\alpha ^21}}+\frac{2\alpha ^2+1}{\sqrt{2\alpha +1}(4\alpha ^21)}\right)\hfill \\ & \sqrt{4\alpha ^21}.\hfill \end{array}$$
(12)
The splitting given by Eq. (10) differs from that of Eq. (1) in Ref. Mooij,  by a factor $`\lambda (\alpha )`$ which explicitly accounts for the dependence of $`E_\pm `$ on $`\alpha `$ and $`g`$.

For stationary states, the current in the qubit loop can be calculated either as the average of the current operator $`\widehat{I}_\mathrm{q}=I_\mathrm{c}\mathrm{sin}(\phi +\theta )`$ over the eigenfunctions, or as the derivative of the energy over the external flux:

$$I_\mathrm{q}=\mathrm{\Psi }_\pm |\widehat{I}_\mathrm{q}|\mathrm{\Psi }_\pm =\frac{E_\pm }{\mathrm{\Phi }}=\pm I_\mathrm{c}f_\mathrm{x}\frac{\lambda ^2(\alpha )}{\pi }\frac{E_\mathrm{J}}{\mathrm{}\omega _0},$$
(13)
where $`\mathrm{}\omega _0=E_+E_{}`$. In equilibrium at finite temperature $`T`$, Eq. (13) readily generalizes to

$$\begin{array}{cc}\hfill I_\mathrm{q}& =\mathrm{\Psi }_+|\widehat{I}_\mathrm{q}|\mathrm{\Psi }_+\rho _{++}^{\mathrm{eq}}+\mathrm{\Psi }_{}|\widehat{I}_\mathrm{q}|\mathrm{\Psi }_{}\rho _{}^{\mathrm{eq}}\hfill \\ & =I_\mathrm{c}\frac{E_\mathrm{J}f_\mathrm{x}\lambda ^2(\alpha )}{\pi \mathrm{}\omega _0}\mathrm{tanh}\left(\frac{\mathrm{}\omega _0}{2k_\mathrm{B}T}\right)\hfill \end{array}$$
(14)
with the density matrix elements $`\rho _{++}^{\mathrm{eq}}=e^{E_+/k_\mathrm{B}T}/Z`$ and $`\rho _{}^{\mathrm{eq}}=e^{E_{}/k_\mathrm{B}T}/Z`$, where $`Z=e^{E_+/k_\mathrm{B}T}+e^{E_{}/k_\mathrm{B}T}`$.

## III Qubit–tank interaction

We propose here to extract information about the quantum dynamics of a flux qubit with the aid of a classical linear high-quality tank circuit, coupled to the qubit via a mutual inductance $`M`$. The tank consists of a capacitor $`C_\mathrm{T}`$, inductor $`L_\mathrm{T}`$, and a resistor $`R_\mathrm{T}`$ which are connected in parallel and driven by a current source $`I_\mathrm{b}(t)`$ (Fig. 4). The problem of coupling a quantum object to a dissipative classical one has no unique theoretical solution. However, if we assume that the classical object is much slower than the quantum one we may solve for the latter’s motion, accounting for the coupling coordinates of the former as mere external parameters.Clark  Here, the characteristic frequency $`\mathrm{\Delta }/h`$ of the qubit is in the GHz range, while the resonances $`\omega _\mathrm{T}`$ of our tank circuit lie below 100 MHz. There exist two different schemes of coupling a tank circuit to the qubit. First we consider direct biasing, where a current $`I_\mathrm{b}(t)=I_0\mathrm{cos}\omega t`$ is fed directly into $`L_\mathrm{T}`$ (Fig. 4).

### III.1 Direct biasing scheme

The voltage across the tank circuit evolves as

$$\ddot{V}+\frac{\omega _\mathrm{T}}{Q}\dot{V}+\omega _\mathrm{T}^2V=M\omega _\mathrm{T}^2\dot{I}_\mathrm{q}+\frac{1}{C_\mathrm{T}}\dot{I}_b(t).$$
(15)
Here, $`Q=\omega _\mathrm{T}R_\mathrm{T}C_\mathrm{T}1`$ and $`\omega _\mathrm{T}=1/\sqrt{L_\mathrm{T}C_\mathrm{T}}`$; $`I_\mathrm{q}`$ is given by (13) or (14), and depends on the qubit flux $`\mathrm{\Phi }=\mathrm{\Phi }_\mathrm{x}+MI_L`$, where $`I_L=L_\mathrm{T}^1V𝑑t`$ is the current in $`L_\mathrm{T}`$ and $`\mathrm{\Phi }_\mathrm{x}`$ is time independent. Below we study the simplest case $`k_\mathrm{B}T\mathrm{\Delta }`$, so that the qubit is definitely in its ground state $`E_{}`$. Then, Eq. (15) takes the form

$$\ddot{V}+\frac{\omega _\mathrm{T}}{Q}\dot{V}+\omega _\mathrm{T}^2V=k^2L\omega _\mathrm{T}^2\frac{d^2E_{}}{d\mathrm{\Phi }^2}V+\frac{1}{C_\mathrm{T}}\dot{I}_b(t)$$
(16)
where $`k^2M^2/LL_\mathrm{T}`$,

$$\frac{d^2E_{}}{d\mathrm{\Phi }^2}=\frac{I_\mathrm{c}^2\mathrm{\Delta }^2\lambda ^2(\alpha )}{(2\pi )^2\left(E_\mathrm{J}^2f^2\lambda ^2(\alpha )+\mathrm{\Delta }^2\right)^{3/2}},$$
(17)
and $`f=[\mathrm{\Phi }_\mathrm{x}+MI_L(t)]/\mathrm{\Phi }_0\frac{1}{2}`$ ($`|f|1`$). Thus, (16) is nonlinear in $`V`$. Since the coupling to the qubit is small one may apply the method of harmonic balance, which is well known in rf-SQUID theory.Likharev  Accordingly, if $`\omega \omega _\mathrm{T}`$, then $`V`$ oscillates with frequency $`\omega `$, while its amplitude $`v`$ and phase $`\chi `$ are slow functions of time: $`V(t)=v(t)\mathrm{cos}[\omega t+\chi (t)]`$. From (16) we obtain

$$\dot{v}=\frac{\omega _\mathrm{T}v}{2Q}+\frac{I_0\mathrm{cos}\chi }{2C_\mathrm{T}}$$
(18)
$$\dot{\chi }=\omega _\mathrm{T}\xi _0\frac{I_0\mathrm{sin}\chi }{2vC_\mathrm{T}}\frac{k^2\omega _\mathrm{T}LI_\mathrm{c}^2}{2\mathrm{\Delta }}\left(\frac{\lambda (\alpha )}{2\pi }\right)^2F(v,f_\mathrm{x})$$
(19)
with the detuning $`\xi _0=(\omega _\mathrm{T}\omega )/\omega _\mathrm{T}`$, and where

$$F(v,f_\mathrm{x})=\frac{1}{\pi }_0^{2\pi }𝑑\varphi \frac{\mathrm{cos}^2\varphi }{\left[1+\eta ^2\left(f_\mathrm{x}+\gamma \mathrm{sin}\varphi \right)^2\right]^{3/2}},$$
(20)
with $`\eta =E_\mathrm{J}\lambda (\alpha )/\mathrm{\Delta }`$ and $`\gamma =Mv/\omega _\mathrm{T}L_\mathrm{T}\mathrm{\Phi }_0`$.

Setting $`\dot{v}=\dot{\chi }=0`$ in (18) and (19) one obtains the stationary tank voltage and phase,

$$v^2\left(1+4Q^2\xi ^2(v,f_\mathrm{x})\right)=I_0^2\omega _\mathrm{T}^2L_\mathrm{T}^2Q^2$$
(21)
$$\mathrm{tan}\chi =2Q\xi (v,f_\mathrm{x}),$$
(22)
where we introduced a flux-dependent detuning

$$\xi (v,f_\mathrm{x})=\xi _0k^2\frac{LI_\mathrm{c}^2}{2\mathrm{\Delta }}\left(\frac{\lambda (\alpha )}{2\pi }\right)^2F(v,f_\mathrm{x}).$$
(23)
We have used Eqs. (21)–(23) to find voltage–flux $`v(f_x)`$, phase–current $`\chi (I_0)`$, and phase–flux $`\chi (f_\mathrm{x})`$ characteristics at resonance $`\omega =\omega _\mathrm{T}`$. We take the qubit parameters $`I_\mathrm{c}=400`$ nA, $`\alpha =0.8`$, $`L=15`$ pH, $`g=100`$, a tank with $`L_\mathrm{T}=50`$ nH, $`Q=1000`$, $`\omega _\mathrm{T}/2\pi =30`$ MHz, and $`k=10^2`$. The $`\chi (I_0)`$ curves for several $`f_\mathrm{x}`$ are shown in Fig. 5. The $`\chi (f_\mathrm{x})`$ and $`v(f_x)`$ curves are shown in Figs. 6 and 7 for various $`\mathrm{\Delta }`$. The sharp dips in Figs. 6 and 7 correspond to the spike in the second derivative of the energy profile in Fig. 3b. Clearly, the width of the dips is correlated with $`\mathrm{\Delta }`$: with the increase of $`\mathrm{\Delta }`$ the width of the dips also increases. The $`\chi (f_\mathrm{x})`$ curves for different $`I_0`$ are shown in Fig. 8. The shape and the value of $`\chi `$ are seen to be very sensitive to $`I_0`$. The dependence of the voltage modulation $`\delta vv(f_\mathrm{x}=0)v(f_\mathrm{x}=10^3)`$ on $`I_0`$ is shown in Fig. 9.

### III.2 Scheme with separate driving coil

In this scheme, a bias flux $`\mathrm{\Phi }_\mathrm{b}(t)=\mathrm{\Phi }_{\mathrm{ac}}\mathrm{sin}\omega t`$ is applied to the qubit loop from a separate coil (Fig. 10). The tank response is similar to (16):

$$\ddot{V}+\frac{\omega _\mathrm{T}}{Q}\dot{V}+\omega _\mathrm{T}^2V=\left[M\omega _\mathrm{T}^2\frac{d^2E_{}}{d\mathrm{\Phi }^2}\mathrm{\Phi }_{\mathrm{ac}}+\stackrel{~}{\mathrm{\Phi }}_{\mathrm{ac}}\right]\omega \mathrm{cos}\omega t,$$
(24)
where $`\stackrel{~}{\mathrm{\Phi }}_{\mathrm{ac}}`$ is the flux which the external coil couples directly into the tank and $`d_\mathrm{\Phi }^2E_{}`$ is given by (17) with $`f=(\mathrm{\Phi }_\mathrm{x}+\mathrm{\Phi }_{\mathrm{ac}}\mathrm{sin}\omega t)/\mathrm{\Phi }_0\frac{1}{2}f_\mathrm{x}+f_{\mathrm{ac}}\mathrm{sin}\omega t`$. Rewriting the first term on the rhs of (24) as

$$k\frac{I_\mathrm{c}^2}{\mathrm{\Delta }}\sqrt{L_\mathrm{T}L}\left(\frac{\lambda }{2\pi }\right)^2\omega _\mathrm{T}^2\omega \mathrm{\Phi }_0f_{\mathrm{ac}}G(t)$$
(25)
makes its time dependence manifest:

$$G(t)=\frac{\mathrm{cos}\omega t}{\left[1+\eta ^2\left(f_\mathrm{x}+f_{\mathrm{ac}}\mathrm{sin}\omega t\right)^2\right]^{3/2}}.$$
(26)
The advantage of a separate driving coil is that one can effectively decouple the tank from the fundamental harmonic of the bias, since the qubit signal $`G(t)`$ contains not only $`\omega `$ but also $`2\omega `$, $`3\omega `$, etc. At $`f_\mathrm{x}=0`$, $`G(t)`$ contains only odd harmonics. This can be used to find the degeneracy point in practice, e.g., by tuning $`f_\mathrm{x}`$ so that the tank response vanishes (reaches its maximum) at frequency $`2\omega `$ ($`3\omega `$). We have studied the higher harmonics by solving Eq. (24) numerically with $`\stackrel{~}{\mathrm{\Phi }}_{\mathrm{ac}}=0`$ for $`\omega =\omega _\mathrm{T}/2`$ and $`\omega =\omega _\mathrm{T}/3`$, see Fig. 11. Since the full amplitudes contain contributions from all harmonics, at $`\omega =\omega _\mathrm{T}/2`$, $`f_\mathrm{x}=0`$ one observes a finite dip rather than a zero.

## IV Requirements on noise sources

Figures 6 and 7 clearly reveal the quantum nature of the flux qubit within a range $`|\delta f_\mathrm{x}|510^4`$ from the degeneracy point $`f_\mathrm{x}=0`$. Therefore, the unavoidable external flux noise coupled to the qubit must be much smaller than this value. The most important sources are the Nyquist noise $`I_\mathrm{n}=\sqrt{4k_\mathrm{B}T/R_\mathrm{T}}`$ and the current noise $`I_\mathrm{a}`$ of the preamplifier. The former generates the qubit-flux noise $`\mathrm{\Phi }_\mathrm{n}=MI_\mathrm{n}Q\sqrt{B}`$, where $`B=\omega _\mathrm{T}/2\pi Q`$ is the tank bandwidth. With $`T=20`$ mK and the tank parameters of Section III, one gets $`\mathrm{\Phi }_\mathrm{n}810^6\mathrm{\Phi }_0`$. For $`I_\mathrm{a}=10^{14}`$A$`/\sqrt{\mathrm{Hz}}`$, we estimate the corresponding flux noise as $`\mathrm{\Phi }_\mathrm{a}=MI_\mathrm{a}Q\sqrt{B}710^6\mathrm{\Phi }_0`$. Thus, the noise these sources couple to the qubit is at least two orders smaller than the peak widths in Figs. 6 and 7. On the other hand, these sources give rise to directly detected voltage noise across the tank circuit. The thermal tank noise is $`V_\mathrm{n}=I_\mathrm{n}\omega _\mathrm{T}L_\mathrm{T}Q\sqrt{B}17.6`$ nV. The noise due to $`I_\mathrm{a}`$ is $`V_{\mathrm{a1}}=I_\mathrm{a}\omega _\mathrm{T}L_TQ\sqrt{B}16`$ nV. And finally, if we take $`V_{\mathrm{a2}}=40`$ pV$`/\sqrt{\mathrm{Hz}}`$ for the preamplifier’s own voltage noise, we get $`V_{\mathrm{a2}}\sqrt{B}7`$ nV for its contribution in the tank bandwidth. Comparing these values with the voltage modulation in Figs. 7, 9, and 11, we see they are at least several times smaller than the qubit signal.

## V Conclusion

We have shown that IMT can be used for low-frequency characterization of the ground (in general: equilibrium) state of a flux qubit. The method allows determining the tunnel splitting between qubit states for a broad class of devices; with the term “flux (as opposed to phase) qubit” we stress that the two states must differ not only in Josephson phase, but in actual magnetic flux visible to the outside.Makhlin  The design exclusively employs present-day technology, and the expected noise levels have been shown not to disrupt the measurement. On the qubit time scale, the method is a quasi-equilibrium one; as such, it can determine $`\mathrm{\Delta }`$ but not, e.g., the “$`T_2`$” dephasing time.amb  Efforts to adapt IMT to this and related problems such as qubit readout and control are underway.

###### Acknowledgements.

We thank A.Yu. Smirnov and A.M. Zagoskin for detailed discussions.