text:
# Searching for Variability in the Globular Cluster Messier 4Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with proposal GO-8679.

## 1 Introduction

It is well known that binary stars are important for the determination of stellar masses. In addition, knowledge of the population of these objects in globular clusters (GCs) is fundamental to understanding the evolution and dynamical history of GCs. Even a small fraction of binary stars in a cluster can contribute greatly in preventing or delaying core collapse of the cluster. A full treatment of this topic can be found in Hut et al. (1992).

This paper presents the results of a search for variability in faint objects in a time-series dataset taken with the Hubble Space Telescope (HST) in a field of the globular cluster M4 (NGC 6121; $`\alpha =16^h23^m54\stackrel{\mathrm{s}}{\mathrm{.}}61`$, $`\delta =26\mathrm{°}32\mathrm{}23\stackrel{}{\mathrm{.}}93`$ (J2000)). Even though the the data were not primarily collected to search for variable stars, the potential for the discovery of such objects did exist, since it included 246 separate observations of the same field over about two and a half months. These observations represent the deepest ever imaging of a globular cluster, reaching the faintest part of the white dwarf sequence ever observed. In the context of this project, it allowed for the study of variability within magnitude ranges that have seldom been explored. Another recent study was performed with HST by Albrow et al. (2001) in which nearly continuous observations for 8.3 days were taken in the core of the globular cluster 47 Tucanae for the primary purpose of searching for planetary transits. Although no such transits were found, the study did identify 114 variables, most of which are believed to be binary systems.

In the color-magnitude diagrams of globular clusters, observations of a binary star sequence 0.75 magnitudes brighter than the main sequence are rare. Even in cases where one is found, many stars in the sequence are actually optical doubles, due to the effects of crowding in globular cluster images. However, binary stars have been found in globular clusters (see, e.g., Hodder et al. (1992), Yan and Mateo (1994), Rubenstein and Bailyn (1996), and Kaluzny et al. (1997)). In fact, W UMa type eclipsing binaries are among the most common type of short-period variable stars found in GCs (Rucinski, 2000).

In addition to binary systems, the data were probed for other types of variable objects, such as pulsating stars, planetary transits, as well as nonperiodic objects such as flare stars, and supernovae in background galaxies. The study of these objects is vital to research in diverse areas such as stellar populations, cosmology, and the evolution of stars, clusters, and galaxies.

The following section deals with the observations of M4 and the reduction of the data. Section 3 discusses what one may expect to find in these data based on previous studies of M4, as well as other globular clusters observed with HST. Section 4 describes the process of selecting candidate variable stars, and section 5 discusses how these candidates were searched for periodicity. Section 6 explains how the reliability of results obtained from the period searches were tested. Section 7 describes searches for nonperiodic variable stars, planetary transits, and supernovae. Finally, section 8 presents and discusses results of the analysis.

## 2 Observations and Data Reduction

The data are taken from Hubble Space Telescope (HST) observations of M4 in the program GO 8679 (cycle 9), using the Wide Field Planetary Camera 2 (WFPC2). Observations were taken over 123 orbits, covering a period of 68.2 days: 1 January - 9 April, 2001, in the natural HST photometric filters F606W and F814W (Holtzman et al., 1995), which will be simply referred to here as $`V`$ and $`I`$, respectively. The dataset consists of $`98\times 1300`$ s exposures in the F606W filter ($`V`$ band) and $`148\times 1300`$ s exposures in the F814W filter ($`I`$ band), all taken in the same field of M4, about 6 core radii ($`5\mathrm{}`$) from the cluster center.

### 2.1 Preprocessing and Combined-Image Photometry

The images were preprocessed according to recipes given in Stetson (1987) and Stetson et al. (1998). Crowding was not an issue on these images, due to their distance from the cluster center, as well as the high resolution of the HST images. A master star list, frame-to-frame coordinate transformations, and photometry of the combined image of the field were completed as discussed in Richer et al. (2002), which resulted in color-magnitude diagrams that reach apparent magnitudes of $`V30`$ and $`I28`$.

### 2.2 Single-Image Photometry

To create time-series data for each star found in the M4 field, it was necessary to perform photometry on each individual exposure. To accomplish this, the following steps were taken: (1) the known frame-to-frame transformations were used to create a combined median image of the field for each chip on WFPC2. This was performed using the MONTAGE2 software written by Peter Stetson (1994), which transforms the coordinates of each individual image to match a reference coordinate system, and then produces a combined image. (2) Once the combined median image was created, aperture photometry was performed on all objects using DAOPHOT. Subsequently, point spread function (PSF) photometry was performed using the ALLSTAR program in order to identify stellar centroid positions. This was done by employing PSFs built specifically for each chip and filter on WFPC2 by Stetson (2001). (3) The output photometry and frame-to-frame transformations were then provided to ALLFRAME, again written by Stetson, in order to perform PSF photometry on each individual exposure (Stetson, 1994; Turner, 1997). Once this single-image photometry was completed, the location of any given star and its accompanying photometry was specified for every single exposure taken in the M4 field.

Since HST is above the Earth’s atmosphere, it is not subject to seeing effects experienced by its ground-based counterparts. However, it is susceptible to slight focus changes due to the telescope moving in and out of the Earth’s shadow as it orbits, which causes temperature fluctuations (i.e. breathing). As a result, the full width half-maximum (FWHM) of the stellar PSFs on the frames can vary slightly between individual images. The result is that small offsets in the magnitudes calculated from performing PSF photometry will exist from frame to frame. These must be corrected before performing any analysis on the data. To accomplish this, the magnitudes that are output by ALLFRAME for every image on a given chip were compared to those of a reference image, which was chosen to be the first image taken in the dataset111Use of the combined median image as a reference was also explored, but there was a negligible difference found in the final magnitudes and magnitude errors obtained, compared with those found through use of a single exposure as the reference image.. The magnitude offset between frames was calculated by taking the weighted mean of the magnitude difference for each star common to both frames.

The usual way of rejecting data points consists of fitting the distribution of magnitudes of a given star with a Gaussian profile, and throwing out any data points at a magnitude greater than $`3\sigma `$ or so from the mean of the distribution, and perhaps reiterating this process a number of times. When looking for variability, objects searched may include those which vary greatly in magnitude, such as eclipsing binary stars with deep minima, or supernovae. By performing this so-called “$`k\sigma `$ clipping” type of data point rejection, data points that may be rejected could actually be caused by real effects due to a significantly varying object. In order to avoid the loss of potentially useful data, we based the rejection of a given data point on its magnitude *error* rather than its magnitude. A “bad data point” is considered to be one which has a magnitude error that is significantly larger than the rest of the data points for the corresponding star’s light curve. To illustrate how “significantly large” was defined, Figure 1 shows plots of root-mean-square magnitude (RMS or $`\sigma `$) vs. mean magnitude in $`I`$ and $`V`$ for each chip. It can been seen that the scatter of data points begins to increase significantly at a well defined “elbow” in the distribution, fainter than which the RMS increases rapidly. With a few exceptions, the maximum scatter of the light curve of any star found in the dataset is approximately $`0.1`$ magnitudes. To be conservative, this is the value chosen as the cutoff magnitude error (as returned by DAOPHOT for a single measurement), so that any individual data point with an error above this value was rejected, and any data point with an error below this value was retained.

## 3 Expectations

In a past search for binary stars in M4, seven W UMa binary stars thought to be cluster members were found by Kaluzny et al. (1997) in an $`8\stackrel{}{\mathrm{.}}8\times 8\stackrel{}{\mathrm{.}}8`$ field centered on the cluster with apparent magnitudes reaching as faint as $`V=17.7`$. The field observed in the current dataset is well within this area. However, unlike the above study, the magnitude range that can be probed for variability in the current study ranges from $`20V26`$, so all variables found in that study are saturated in the current dataset. However, there is no reason to expect that the existence of binaries is exclusive to bright magnitudes. In a recent study of the core of the globular cluster 47 Tucanae (47 Tuc), which probed magnitudes as faint as $`V=25`$, 11 detached eclipsing binaries and 15 W UMa systems were discovered in a sample of 46422 main sequence stars, resulting in an observed binary frequency of $`0.056\%`$ (Albrow et al., 2001). In this study of M4, 2102 stars have been analyzed. If no assumptions are made about the dependence of the distribution of binary systems with cluster radius, one would expect that only one or two binary stars should be discovered in the field observed based on the observed frequency in 47 Tuc. This already small number is decreased if it is assumed that the majority of binaries are expected to reside close to the cluster center due to mass segregation, since the M4 field is located at about 6 core radii from the cluster center. The picture becomes even less optimistic considering that only 9 of the 15 W UMa eclipsing binaries found in the 47 Tuc study are below the main sequence turn-off point. However, 71 BY Draconis variables, which are thought to be members of binary systems, were also discovered, making up $`0.15\%`$ of the sample. This converts to a slightly more optimistic expectation of 3 or 4 of these objects being discovered in the M4 field.

As a way of gauging the number of supernova events that are expected to be seen in the current dataset, a comparison is made to a recent study of the Hubble Deep Field (HDF), in which two supernovae were discovered (Gilliland et al., 1999). As with M4, the HDF was observed with the WFPC2 camera on HST, and just as deep, with a limiting magnitude of $`I28`$ (Flynn et al., 1996; Richer et al., 2002). However, unlike the HDF, the M4 images contain zodiacal light, in addition to many stars, due to the cluster’s proximity to the Galactic plane and ecliptic. This results in much higher background levels compared to the HDF, significantly decreasing the number of observable faint galaxies in the field. As will be discussed in section 7.2, 170 galaxies were identified in the M4 images. By contrast, $`2100`$ galaxies were identified in the HDF (Williams et al., 1996). Scaling the discovery of 2 supernovae in the HDF to the current dataset, it is estimated that $`0.08`$ extragalactic SNe would be expected in this set of images.

There also exists the possibility of discovering variable objects which have not yet been discussed, such as cataclysmic variables (CVs), flare stars, transiting planetary-sized objects, or perhaps a new species of variable star. CVs are very unpredictable, making it difficult to quantify how many are expected to be seen in the M4 field. However, in globular clusters, which have a high frequency of stellar interactions, it is expected that an abundant population of CVs would be found (Di Stefano and Rappaport, 1994; Grindlay et al., 2001; Knigge et al., 2002; Townsley and Bildsten, 2002), and indeed, CVs have been found in globular clusters (see, e.g., Shara et el. (1996); Cool et al. (1998); Neill et al. (2002)). A relatively small sample of stars may limit this population in the current observations, but the possible existence of CVs in this dataset is not discounted.

The main science driver for the above-mentioned study of 47 Tuc was to define the frequency of close-in gas giant planets orbiting main sequence stars (Gilliland et al., 2000). This study reported a null result, even though simulations that were carried out clearly showed excellent sensitivity to planetary transits. However, the observations were centered on the core of 47 Tuc, which may be too hostile an environment for planets or protoplanetary disks to survive. The M4 images were taken far from the core of the cluster, where the number density of stars is much smaller, and thus may provide a more suitable environment for planet formation. On the other hand, most extrasolar planets are found to be orbiting stars with relatively high metallicity. M4 is a more metal poor globular cluster than 47 Tuc, perhaps indicating that a smaller fraction of planets will be found in the current dataset. If the null result from 47 Tuc is any indication, it is expected that in the present study of M4, which contains only $`4.5\%`$ of the number of stars in the 47 Tuc dataset, no planetary transits will be found. In addition, planetary transits found by Charbonneau et al. (2000), and very recently by Konacki et al. (2003), show eclipse depths of $`1\%`$ for close-in Jupiter-size planets. For a relatively bright, non-variable star in this dataset, the scatter in its light curve is typically $`0.014`$ magnitudes. This is most likely too large to confirm a $`1\%`$ eclipse depth should one exist in the light curve. On a more positive note, however, M4 is known to contain a millisecond pulsar with a planetary companion (Lyne et al., 1988; Backer, 1993; Rasio, 1994; Sigurdsson, 1995; Thorsett et al., 1999; Sigurdsson et al., 2003). The presence of a planet in this unusual system may be suggesting that planets are rather common in at least this globular cluster.

A further consideration is that the stars for which planetary transits have been observed are stars of solar-like mass and radius. However, the magnitude range of the stars that we have examined for variability is $`19V26`$. This corresponds to a range in stellar masses of $`0.6M_{\mathrm{}}M_{}0.1M_{\mathrm{}}`$. These stars have smaller radii than the Sun, and one would expect that for planetary occulations of such stars, the eclipse depth that would be seen in the star’s light curve would be significantly deeper than for the two transiting planets found to date. In fact, once the stellar mass reaches $`0.1M_{\mathrm{}}`$, close to 100% eclipse depths could be expected for a system than is close to edge-on, assuming such a system contains a planet with Jupiter’s radius.

## 4 Candidate Selection

After reducing the data and rejecting bad data points as described in section 2.2, variability searches were conducted on the light curves of the stars in the dataset. Described in the subsections that follow are various criteria that were employed in order to narrow down the number of candidate variable stars in this dataset. After satisfying the criteria, any star that contained in its light curve less than six data points was immediately rejected. This may seem like a rather generous cutoff value, but when searching for possible supernova events, six points could be enough to determine if it is indeed a viable candidate. Stars were also excluded from the candidate list if:

* The image of the star is saturated on the CCD;
* The star is close enough to the edge of the image so that it can partly appear and disappear over time due to the observational dithering pattern;
* The stellar image is affected by light from a nearby saturated star; or
* The star is partially or exactly aligned with a diffraction spike of a bright or saturated star.

All remaining candidates that were not excluded due to the above criteria are presented and included in tables accompanying each subsection. In total, 2102 stars were analyzed for variability.

### 4.1 Statistical Outliers

A simple identification was first made of outliers in standard deviation (RMS magnitude) as a function of magnitude. The light curves of these stars have a larger scatter than is expected for stars of similar magnitude, and are thus potential variable star candidates. These RMS magnitude outlier stars, not including those rejected due to a lack of data as described in section 4, are shown in Figure 2, and candidates chosen by using this method are listed in Table 1. In this and other tables, coordinates of the objects listed are based on work on M4 by Ibata et al. (1999), and coordinates for stars not included in that study were calculated using the CCMAP and CCTRAN tasks on the IRAF analysis package.

While the above method of choosing variability candidates based simply on the scatter of the light curves of stars is a good first pass in the candidate search, it is not terribly robust against data corrupted by cosmic rays and warm pixels, which are quite common in HST images. A more robust candidate selection criterion that attempts to account for these effects is the variability index statistic $`J`$, a magnitude-independent measure of the confidence that a given star shows real signs of variability (Stetson, 1996). Another attractive feature of this method is that it makes use of data in *all* filters employed in the dataset when determining the variability index.

Values of $`J`$ were computed for every star, and those with values of $`J`$ that were obvious outliers in the variability index distribution as a function of magnitude were considered variable star candidates. These distributions are shown in Figure 3, where stars selected as candidates are plotted as open circles. None of these stars passed the additional criteria described in section 4, and so no variable star candidates arose as a result of using this method.

### 4.2 Main Sequence and White Dwarf Sequence Outliers

As explained in Section 1, not many binary star sequences are found parallel to the main sequence in globular cluster CMDs. Looking at the color-magnitude diagram of M4 (Figure 4, left; Richer et al. (2002)), one can see no obvious second sequence in the magnitude range observed, other than the expected faint white dwarf (WD) sequence. Nevertheless, there are stars that appear to be outliers to the main sequence, lying above it, or between it and the white dwarf sequence. The CMD on the right hand side of Figure 4 is the same as that on the left, but stars which obviously do not lie on the MS or WD sequence are plotted as open circles. From these stars, period searches were performed on those which were not excluded according to the criteria listed in section 4. These remaining candidates are listed in Table 2.

The current dataset also includes white dwarf stars that are coincident with the ZZ Ceti instability strip. Figure 5 shows the M4 CMD, with the ZZ Ceti instability strip bordered by a black box. Unfortunately, the exposure times of the individual images, which are 1300 seconds for each single observation, are on the order of the typical ZZ Ceti period. Any variability in these stars would be “washed out”, since the integration time of each observation is too long to capture the star at various stages in its pulsation cycle. It is useful however to list the stars which are thought to lie within the ZZ Ceti instability strip, as future observations can be made specifically to study these objects. Table 3 lists these ZZ Ceti variable star candidates as well as their coordinates, magnitudes and colors, and Figure 6 presents finder charts for these three candidate variables.

## 5 Periodicity Analysis

Ever since the landmark paper by Stellingwerf (1978) describing the technique, phase dispersion minimization (PDM) has been used extensively and successfully in the discovery of countless numbers of variable stars and binary systems. Unlike Fourier transform analysis, which is particularly well-suited for sinusoidal-type variations, it is not sensitive to the shape of a star’s light curve. In addition, its implementation is very straightforward. However, determination of the statistical significance of results obtained from use of PDM is not so simple, as is the case for Fourier techniques. It has been shown that the PDM probability distribution does not follow an $`F`$ distribution (Schwarzenberg-Czerny, 1997), as originally thought. Methods used here to estimate the reliability of findings made by PDM analysis will be discussed in the next section. PDM analysis was performed on all stars in the dataset, and particular attention was paid to candidate variable stars.

## 6 Reliability Tests

In an effort to understand the reliability of the PDM analysis performed on candidate variable stars, a number of artificial stars were created by sampling different types of variable star light curves at the same observation times as the M4 dataset. These artificial stars were then added to the WF2 images, and PDM analysis was performed on the light curves extracted from the resulting images.

Sinusoidal light curves were simulated with a range of periods, amplitudes, and magnitudes. The range of magnitudes simulated were $`17I28`$ and $`19V28`$. The bright end of these ranges correspond to the saturation limits of the dataset in each bandpass. The faint end corresponds to the faintest magnitude reached in the combined $`I`$ image. Six periods were simulated: 0.085, 0.16, 0.42, 2.1, 5.6, and 15.6 days. The shortest period simulated corresponds to a value just above the Nyquist frequency for this dataset, and the longest period simulated corresponds to approximately the duration of data collection before the occurence of the 44.4-day gap in the observations. Amplitudes of 0.002, 0.02, and 0.2 magnitudes were simulated. In all, 20181 sinusoidal light curves were simulated for each filter.

In addition, two types of eclipsing binary star light curves were simulated by employing the widely-used Wilson-Devinney (W-D) code (Wilson and Devinney, 1971). Using this program, Algol and W UMa-type systems were simulated, each with a 0.42-day period. The simulated Algol system has an orbital inclination ($`i`$) of $`67^{}`$, a relative orbital ellipse semimajor axis length ($`a`$) of $`351.6R_{}`$, and a primary eclipse depth of $`0.3`$ magnitudes. The simulated W UMa system has $`i=47.5^{}`$, $`a=1.88R_{}`$, and a primary eclipse depth of $`0.2`$ magnitudes. A period of 0.42 days was chosen because it exhibits a case of full phase coverage for this dataset. In addition, it is within the range of periods typical of close binary stars. In all, 1922 binary star light curves were simulated and added to the images in each filter. This light curve data was used in conjuction with the ADDSTAR program in DAOPHOT (Stetson, 2001) to add artificial stars to the actual M4 images at different magnitudes. This ensured that the output light curves would have noise levels typical of real stars in the image with similar magnitudes. Photometry was performed on these artificial stars in the same way as on the original images (see section 2.2).

For the artificial sinusoids, a recovered period within 1% of the correct input period was considered to be a successful recovery. Results from the simulation of the artificial sinusoidal light curves are shown in Figure 7. For the curves in both bandpasses with an amplitude of 0.2 magnitudes, PDM analysis recovers the input period almost $`90\%`$ of the time down to a magnitude of 24 in $`I`$ and 26 in $`V`$ (at which point the recovery rate drops to $`64\%`$ at its lowest). There is no dependence of the recovery rate on input period apart from the 0.085 and 15.4-day cases. This gives an idea of the periods to which the PDM analysis is sensitive at this amplitude given the time sampling of this dataset. For the case of light curves with 0.02 and 0.002-magnitude amplitudes, the ability of PDM analysis to recover the correct input period is significantly lessened. For an amplitude of 0.02 magnitudes, the recovery rate in $`I`$-band is very poor regardless of period or magnitude. For $`V`$-band, the period search is also discouraging, although the 0.42 and 2.1-day period recovery rate is slightly higher for brighter stars. In the case of the 0.002-magnitude amplitude curves, the PDM analysis recovery rate is even less successful; the input periods are generally not recovered, regardless of the period or magnitude of the light curve. This is not surprising, since the scatter of the light curves is on the order of this amplitude or greater (see Figure 1), even for bright stars.

The results of the tests for the artificial binary star light curves are found in Table 4. These tables show the number of stars added at each magnitude, the number recovered, and the percentage of those stars successfully recovered in each filter. A recovery was considered to be successful in this case if the period found through PDM analysis is within 1% of either the correct input period (0.42 days), half the input period, or $`\frac{3}{2}`$ the input period. It was found that these are the three major groups within which the periods of the artificial binary stars were recovered. This is because the light curve folded at half or three halves the input period will cause the primary and secondary eclipses to coincide. This caused the PDM analysis to occasionally choose these fractions of the input period. For these artificial binary star reliability tests, the PDM analysis is highly successful in recovering the input periods of the Algol-type light curves down to limiting magnitudes of $`I23`$ and $`V24`$. For light curves fainter than this, less than half are successfully recovered. For the W UMa-type light curves, these limits are $`I24`$ ($`65\%`$ recovered) and $`V25`$ ($`89\%`$ recovered). Below these limiting magnitudes, almost none of the stars were recovered with light curves consisting of more than 6 data points (see section 2.2), let alone with the correct input period.

## 7 Searches for Non-Periodic Variable Stars

In addition to periodic variables, there may exist objects in the M4 field that vary in magnitude, but not periodically, such as cataclysmic variables, supernovae, and flare stars. There is also the possibility of the presence of objects whose period is too long compared with the time frame in which the observations were taken (or too short relative to the time sampling or exposure times), so that their periodicity may not be evident, as in the case of transits of stars by planets. Clearly, tools such as PDM are not useful in searches for these objects. The variability candidates chosen by some of the methods outlined thus far in this section (MS/WD outliers and variablity index, to be specific) *are* valid for these cases and those candidates chosen by these methods were searched for the existence of non-periodic variability in the M4 field.

### 7.1 Planetary Transits and Median-Smoothed Light Curves

As an estimate of the typical period lengths for planetary transits that can be expected to be seen in this dataset, we take an example of a Jupiter-sized planet. From the reliability tests described in section 6, we can expect to find periodic transits that have eclipse depths of about 20%. To get a 20% dip in the light curve of such a system, the radius of the parent star must be $`0.2R_{\mathrm{}}`$. This translates into a mass $`0.23M_{\mathrm{}}`$ (Montálban et al., 2000), and a magnitude $`V23.2`$ for this dataset. We set a minimum detectable transit time $`5600\text{ s}`$ based on the minimum separation in time between three consecutive exposures in the light curve. We estimate the orbital speed implied by such a transit by calculating the length of the transit to be the time it takes for the planet to pass from second contact to third contact. This implies a maximum orbital speed of $`31.6\text{ km\hspace{0.17em}s}^1`$. Assuming a circular orbit ($`v_{orb}=2\pi a/P`$), the orbital period is then given by:

$$P=\frac{2\pi G(M_{}+M_{pl})}{v_{orb}^3}71\text{ days.}$$
(1)
This is longer than the length in time over which the dataset was taken, so only one eclipse would be expected to be found. This is also a lower limit to the period, for the following reasons. From section 6, the lower magnitude limit at which we successfully recover periods over 90% of the time is $`V26`$. This corresponds to $`0.1M_{\mathrm{}}`$ for main sequence stars. A lower mass star has a smaller radius, and so the orbital velocity of the planet must be slower in order to have a transit time that can be detected in this dataset, which results in a longer orbital period. The same argument holds for a higher orbital inclination, since the separation between second and third contacts is smaller. Thus, periodic signatures of planetary transits would not likely be found for such systems in this dataset. Searches for potential one-time transits in this dataset are discussed below.

One method used to search for one-time planetary transits and variable star candidates with non-periodic variability involved smoothing the light curve of each star found in the dataset. This was done by replacing each data point in the light curve with the median of that point and its two nearest neighbouring points. This has the effect of reducing the scatter in the light curve, so that real variations in the light curve would be more pronounced. In addition, median smoothing was performed using four and six nearest neighboring points, to search for longer-term variability in fainter stars having noisier light curve data.

This method was executed as follows. The faintest data point in each median-smoothed light curve was selected, and the ratio of the difference in magnitude between that point and the mean magnitude of the light curve to the RMS scatter of the light curve was calculated. Outliers to the distribution of this ratio represent data points that deviate more than is expected from the scatter of the smoothed light curve, which would indicate the possibility of a real temporary variation in the star’s light curve. However, all the outliers found were simply due to saturated stellar images, diffraction spikes, or the star being too close to the edge of the image, causing it to disappear and reappear in the frame over time. Thus, no candidates for planetary transits, or other type of non-periodic variable star, were found using this method.

### 7.2 Searching for Supernovae

In addition to the median smoothing of light curves, additional techniques were employed to search for possible supernovae (SNe) occuring in the M4 field. Within the dataset, there is a 44.4-day gap in the data. For this section, and further discussions pertaining to searching for supernova with this dataset, *epoch 1* will refer to the data taken before this gap, and *epoch 2* to the data taken after the gap. The data taken during epoch 1 and epoch 2 were combined separately for each bandpass using MONTAGE2 , and PSF photometry was performed on each resulting image using DAOPHOT and ALLSTAR in exactly the same way as described in section 2.2. The combined epoch 1 images used 100 images in $`I`$ and 66 in $`V`$, taken over 14 days. The combined image for epoch 2 used 48 images in $`I`$ and 32 in $`V`$, taken over 9 days. Two methods employing these combined images were used to search for supernova candidates. The first involved comparing the magnitudes of stars between the two epochs. Outliers in the distribution of magnitude difference as a function of magnitude, shown in Figure 8, were chosen as possible candidates. However, none of the candidates resulting from this search passed the exclusion criteria described in section 4. For the second method, background galaxies were identified, and the second epoch image was subtracted from that of the first epoch (after performing sky subtraction and flux normalization for each combined image) in order to search for residual flux at the position of these identified galaxies. There was no evidence for any residual flux greater than about $`1\sigma `$ above (or below) the noise levels for the image, except where the nucleus of the galaxy was too bright, in which cases the quality of the subtraction was poor. However, some of these bright nuclei show stellar-like light profiles, perhaps indicating the presence of active galactic nuclei. This would make these galaxies interesting objects for possible further study.

## 8 Results and Discussion

After candidates were selected as described in section 4, their light curves were examined more closely and were analyzed using PDM to search for periodicity. Of all the variable star candidates listed in the tables in the previous sections, only one of these stars (PC, ID #97) showed real signs of variability in its light curve. Even for this remaining candidate star, indicators of periodicity were not significant enough to be conclusive. It is therefore possible that this remaining candidate is non-periodic, or else there is not enough data in the light curve of that star to decisively find a period for it through PDM analysis. This amounts to an observed photometric variability fraction of $`0.05\%`$ for this dataset.

### 8.1 Remaining Candidate

There is the still unresolved issue of the remaining candidate star found in the PC field of the observations. There is no clear reason that this star should not be considered as variable: it is relatively bright, it is not close to any image defects (such as cosmic rays or diffraction spikes from saturated stars) and is far from the edge of the image. Figure 9 shows a combined $`I`$-band image of the PC chip observations, in which the star in question is circled. In addition, it has been identified as a cluster member by its proper motion. The notable feature in the light curve of this star is the apparent rise in amplitude by $`0.1`$ magnitudes within the first few days of the onset of data collection, as seen in Figure 10. This peak is seen to occur in both bandpasses at the same time, further supporting the claim that this object is indeed variable. A clear second peak in brightness is not seen, which could mean that the star is not periodic. However, it may be that this object has a period such that a second peak occured during a gap in the observations, or else would probably have to be modulated in amplitude over time in order to explain the lack of additional peaks. Inspection of its light curve seems to indicate a possible period in the range of 5 - 9 days, but periodicity analysis performed on this star showed no significant result.

This star has a magnitude on the combined image of $`V=22.6`$ mag and color $`VI=0.71`$ mag, which places it at the bright end and just red of the white dwarf sequence in the CMD of M4, as shown in Figure 11. The location of this star in the color-magnitude diagram indicates that it may be a binary star, composed of a low-mass main sequence star and a white dwarf. In order to explore this possibility further, a simple decomposition of this candidate into two stars was performed. We investigated the magnitudes and colors of two separate stars, lying on either the main sequence or white dwarf sequence, that would be needed to produce a single observed star of magnitude and color similar to this candidate. This gives further insight as to the types of stars that might exist in a binary system found in the position of this candidate on the M4 CMD. We find that a binary system consisting of a white dwarf with magnitude and color $`V22.7`$, $`VI0.4`$ and a faint, late-type main sequence companion star, with upper limits to its magnitude and color of $`V26.8`$ and $`VI3.4`$, respectively, result in an object with a position on the CMD that is very near to that of the candidate star. This is illustrated in Figure 12. The positions of the two speculated binary star constituents are shown as the smaller open circles, and the resulting binary is shown as a filled circle. A large open circle is centered on the position of the actual candidate star. As one can see, the composite star is very close to the position of the candidate. Because of the proximity of the candidate variable to the white dwarf cooling sequence, it was expected that such an object would be dominated in luminosity by a white dwarf. The position of the hypothetical lower main sequence constituent would indicate that this object could be a late M-type star. However, the magnitude and color found for this companion are lower limits, since any companion that is fainter on the main sequence contributes little to the hypothetical binary system. It is then possible that this companion could be an even fainter, lower-mass main sequence star.

It may be possible that this object is a cataclysmic variable (CV). It is thought that such objects can be identified in the color-magnitude diagrams of globular clusters, between the main sequence and the white dwarf sequence (Townsley and Bildsten, 2002). In addition, those that are found close to the white dwarf sequence will be dominated in magnitude by the white dwarf component. The light curve of the object does not seem to show the characteristics of an eclipsing binary system. However, the possibility does exist that it is a face-on binary system that includes a BY Draconis variable. The period range of 0.4-10 days of the BY Dra stars found in the 47 Tuc study does not discount this as a possibility. In addition, that study reported amplitudes for these stars ranging from $`V0.001`$ \- $`0.061`$. This is not much smaller than the $`0.1`$ magnitude increase seen in the candidate object’s light curve. A further possibility is that the observed increase in brightness is due to flaring on the star’s surface, which would not necessarily be periodic, thus explaining the null result of the PDM analysis. However, this is unlikely, since the timescale typical of stellar flares rising to maximum brightness and diminishing is on the order of hours, whereas the rise in amplitude in the light curve of the candidate star seems to last for days. In any case, more data would be needed in order to determine the nature of this object’s variability. There is another star next to this candidate on the CMD. This star (ID #1), also found on the PC chip image, was very close to the edge of the frame. Its partial appearance and disappearance from the field over time, due to the dither pattern of the observations, was responsible for the observed variations in its light curve, discounting it as a possible variable star candidate.

### 8.2 A Lack of Variability

The lack of eclipsing binary stars in the field can be explained in a number of ways. Firstly, it may be that any binaries that are present in this field are face-on relative to the line of sight, although this possibility is an unlikely one. Another reason may be the fact that the observations were taken in a field that is considerably distant from the core of M4. It is expected that mass segregration in the cluster will, over time, preferentially cause binary systems to sink toward the cluster center, since binary stars are typically more massive on average than single stars. This may not have left many binary systems at the large cluster radius that was observed in this study. However, a considerable number of $`0.6M_{}`$ white dwarf stars are seen in the field (Hansen et al., 2002), as are main sequence stars up to $`0.65M_{}`$ (Richer et al., 2002), so it is not unreasonable to expect, for example, binary star systems comprised of $`2\times 0.3M_{}`$ stars in that same field. Another possibility is that due to the high-density environment of globular clusters, low-mass binaries are more prone to disruption before they have the opportunity to harden over time through interaction with third bodies. Conversely, the relatively large numbers of blue straggler stars in globular clusters (30 have been found in M4 (Rucinski, 2000)) may represent the coalescence of many low-mass binary systems which have acted as a heat source to the cluster in preventing core collapse. This could reduce the number of binaries expected in this magnitude range. However, the real explanation for the lack of eclipsing binaries in this sample may simply be small-number statistics. As speculated in section 3, based on the low fraction of stars in the study of 47 Tuc observed to be in binary systems, the predicted number of these objects recoverable in this dataset was low to begin with.

The small number of background galaxies observable in the M4 images is the most probable reason that no supernovae were found in the field observed. Previous studies of the Hubble Deep Field, whose images are comparable in depth to the M4 field, yielded the discovery of two supernovae (Gilliland et al., 1999). Due to the substantially larger background flux in the M4 observations compared to that of the HDF (see section 3), considerably fewer galaxies were observable, making it unlikely that any supernovae would be found in these observations.

While this research did not produce any confirmed discoveries of variability in the magnitude range $`V19\mathrm{\hspace{0.17em}26}`$ for M4, it is far from having been fruitless. Firstly, there is one interesting candidate variable, for which it seems that more observations will be needed to identify its nature. However, we can place an upper limit to the observed fraction of photometric variability for this study at $`0.05\%`$, within the magnitude ($`19V26\text{ and }17I24`$), amplitude ($`\mathrm{\Delta }V0.1`$ magnitude), and period ($`0.16\text{ days}P7.0\text{ days}`$) ranges to which we were sensitive.

The identification of the white dwarf stars in the M4 field suspected of being ZZ Ceti variables may prove to be of considerable value. The discovery of such stars in a globular cluster is of the utmost importance to understanding the chemical structure and evolution of population II white dwarfs, which currently suffers from a lack of observational data.

We are grateful to Robert E. Wilson for providing the Wilson-Devinney code, and to Puragra Guhathakurta for several very useful discussions on various aspects of this work. The authors also wish to thank the anonymous referee, who contributed many valuable suggestions. J. S. K. received financial support during this work through an NSERC PGS-B graduate student research grant. H. B. R. and G. G. F. are supported in part by the Natural Sciences and Engineering Research Council of Canada. H. B. R. extends his appreciation to the Killam Foundation and the Canada Council for the award of a Canada Council Killam Fellowship. R. M. R. and M. M. S. acknowledge support from proposal GO-8679 and B. M. S. H. from a Hubble Fellowship HF-01120.01 both of which were provided by NASA through a grant from the Space Telescope Science Institute which is operated by AURA under NASA contract NAS5-26555. B. K. G. acknowledges the support of the Australian Research Council through its Large Research Grant Program A00105171.


format: text

text:
# 1 Introduction

## 1 Introduction

Recently, large interest has been devoted to a variety of deformations of special relativity admitting two invariant fundamental scales, the speed of light and an energy scale, to be identified with Planck energy $`\kappa `$ \[1-3\]. These theories aim to describe the dynamics of particles up to the Planck region, where the structure of spacetime may change due to quantum gravity effects. They were first introduced algebraically, investigating quantum deformations of the Lorentz group , and later rederived starting from a set of suitable physical postulates, as for example the request that special relativity is recovered in the low-energy limit .

All these models have in common the assumption that the momentum of a particle transforms nonlinearly under the Lorentz group, but differ in several respects. In fact, first of all, it is possible to construct many different nonlinear representations of the Lorentz group satisfying the postulates mentioned above; moreover, the theory can be defined either in standard spacetime, or in a spacetime with noncommuting coordinates . In the latter case, when defining the classical dynamics through the Hamiltonian formalism, one is forced to use noncanonical Poisson brackets. Of course, depending on the choice of the previous basic assumptions, one obtains models that give rise to different physical predictions.

In order to compare the theory with experiment, one also needs a physical interpretation of the variables that enter in the theory in terms of measurable quantities. This is not always straightforward. We must consider in fact that the models under study, although classical, are supposed to be effective in the quantum domain, so that one cannot expect that all quantities (as for example velocity) can be defined in classical terms.

As we shall discuss in the following, for example, a difficulty in the interpretation of the theory arises from the fact that the transformation rules for the position of a particle depend on its momentum, and hence particles occupying the same position in a reference frame do not necessarily do so in another frame.

We believe that problems of this kind can be understod by recalling the quantum-mechanical origin of the models, which only in the limit of small momenta reduce to special relativity. This implies that the physical quantities should be interpreted accordingly. Clearly, we have no way to probe the sub-Planckian scales for which the theory is defined, and so we do not have to wonder if some of the usual assumption are not valid, provided that the theory is logically consistent.

In this paper we try to clarify some of these questions by discussing the issue of the transformation laws of position and momentum in a Hamiltonian setting, assuming either commuting spacetime coordinates and standard Poisson brackets or noncommuting coordinates and deformed symplectic structure. We treat explicitly the Magueijo-Smolin (MS) model , since this is the simplest algebraically, but our considerations are valid also for different models. In the appendix, we shortly discuss the Lukierski-Nowicki-Ruegg (LNR) model . For simplicity, we consider the case of 1+1 dimensions, but our results can be easily extended to four dimensions.

## 2 The transformation law for commuting coordinates

Let us consider the Hamiltonian formulation for a free particle in $`1+1`$ dimensions. We introduce commuting spacetime coordinates $`q^a`$ and momenta $`p_a`$ ($`a=0,1`$) for a particle, obeying canonical Poisson brackets, so that

$$\{q^a,q^b\}=0,\{p_a,p_b\}=0,\{q^a,p_b\}=\delta _b^a.$$
(1)
Deformed special relativity models \[1-3\] postulate a nonlinear Lorentz transformation law for the momenta, but do not specify the transformation law of the coordinates. More precisely, it is assumed that under an infinitesimal boost of generator $`J`$,

$$\delta p_a=\{J,p_a\}=w_a(p),$$
(2)
where $`w_a(p)`$ is a nonlinear function of the momentum $`p_a`$. For the MS model , the functions $`w_a`$ are

$$w_0=p_1\frac{p_0p_1}{\kappa },w_1=p_0\frac{p_1^2}{\kappa }.$$
(3)
Now, it is easy to see that the transformation law (2) implies that also the coordinates must transform in a nontrivial way under Lorentz transformations. This is a simple consequence of the Jacobi identities. In fact, under a Lorentz transformations $`\delta q^a=\{J,q^a\}`$, but

$$\{\{J,q^a\},p_b\}+\{\{p_b,J\},q^a\}+\{\{q^a,p_b\},J\}=0$$
(4)
and hence

$$\frac{\{J,q^a\}}{q^b}=\{\{J,q^a\},p_b\}=\frac{w_b}{p_a}$$
(5)
Integrating these relations, one finds that

$$\{J,q^a\}=\frac{w_b}{p_a}q^b+f^a(p)$$
(6)
where $`f^a(p)`$ are arbitrary functions of $`p^a`$, satisfying $`f^a/p_b=f^b/p_a`$, that can be set to zero.

From (3) and (6) one has

$$\delta q^0=q^1+\frac{p_1}{\kappa }q^0,\delta q^1=\left(1\frac{p_0}{\kappa }\right)q^0+\frac{2p_1}{\kappa }q^1,$$
(7)
Hence, the transformation properties of the position of a particle depend on its momentum, and two particles occupying the same position in a reference frame can occupy different positions in another. The momentum dependence of the transformations of the coordinates of a particle has been first remarked by Kowalski-Glikman  in a special case, but it seems to have been disregarded by other authors.

Actually, one can even derive the finite transformation properties of coordinates from those of momenta, by requiring covariance. In momentum space $`𝒫`$, the effect of a boost is given by a nonlinear transformation

$$p_ap_a^{}=W_a(p).$$
(8)
The space of coordinates can be identified with the tangent space to $`𝒫`$ and hence for covariance,

$$q^aq^a=Ł_b^aq^b,$$
(9)
where

$$Ł_b^a=\left(\frac{W_a}{p_b}\right)^1.$$
(10)
For the MS model, one has

$$W_0=\frac{p_0\mathrm{cosh}\xi +p_1\mathrm{sinh}\xi }{\mathrm{\Delta }},W_1=\frac{p_1\mathrm{cosh}\xi +p_0\mathrm{sinh}\xi }{\mathrm{\Delta }},$$
(11)
where

$$\mathrm{\Delta }=1+\frac{p_0}{\kappa }(\mathrm{cosh}\xi 1)+\frac{p_1}{\kappa }\mathrm{sinh}\xi ,$$

and $`\xi `$ is the rapidity parameter. It follows that

$$Ł_b^a=\mathrm{\Delta }\left(\begin{array}{cc}\mathrm{cosh}\xi +\frac{p_0}{\kappa }(1\mathrm{cosh}\xi )& \mathrm{sinh}\xi +\frac{p_1}{\kappa }(1\mathrm{cosh}\xi )\\ \mathrm{sinh}\xi +\frac{p_0}{\kappa }\mathrm{sinh}\xi & \mathrm{cosh}\xi +\frac{p_1}{\kappa }\mathrm{sinh}\xi \end{array}\right)$$
(12)
It is easy to check that the infinitesimal version of these relations is (7).

Given the transformation law (9), one can also define a covariant lagrangian,

$$L=p_a\dot{q}^aH(p).$$
(13)
In fact, $`p_a\dot{q}^a`$ is clearly invariant up to a total derivative under the Lorentz transformations (8-9). For a free particle, the Hamiltonian $`H`$ is given by the Casimir invariant of the algebra ,

$$H=\frac{1}{2}\frac{p_0^2p_1^2}{\left(1\frac{p_0}{\kappa }\right)^2}.$$
(14)
It may be interesting to write down the Hamilton equations for (14):

$$\dot{q}^0=\frac{H}{p_0}=\frac{p_0\frac{p_1^2}{\kappa }}{\left(1\frac{p_0}{\kappa }\right)^3},\dot{q}^1=\frac{H}{p_1}=\frac{p_1}{\left(1\frac{p_0}{\kappa }\right)^2},$$

$$\dot{p}_0=\frac{H}{q^0}=0,\dot{p}_1=\frac{H}{q^1}=0.$$

The velocity of a particle obeying these equations of motion is given by

$$v=\frac{\dot{q}_1}{\dot{q}_0}=\frac{p_1p_0p_1/\kappa }{p_0p_1^2/\kappa }.$$
(15)
As discussed in the literature , this definition of velocity is not satisfactory, since it implies that the velocity of a particle depends on its mass, and it is difficult to reconcile this fact with the role of velocity as the parameter of the Lorentz transformations.

A solution to this problem is to use deformed Poisson brackets . This is very natural from the point of view of $`\kappa `$-Poincaré models, since deformed Poisson brackets can be considered as the classical limit of noncommuting spacetime coordinates. In the following section we reconsider the MS model from this point of view.

## 3 The transformation law for noncommuting coordinates

In the MS model, a suitable definition of the velocity of a particle is given by $`v=p_1/p_0`$ . This coincides with $`\dot{q}_1/\dot{q}_0`$ if one introduces the following symplectic structure :

$`\{q^0,q^1\}={\displaystyle \frac{q^1}{\kappa }},\{p_0,p_1\}=0,\{q^0,p_0\}=1{\displaystyle \frac{p_0}{\kappa }},`$ (16)
$`\{q^1,p_1\}=1,\{q^0,p_1\}={\displaystyle \frac{p_1}{\kappa }},\{q^1,p_0\}=0.`$ (17)
The infinitesimal transformations of coordinates can then be deduced from the Jacobi identities as above. They read

$$\delta q^0=q^1+\frac{p_1}{\kappa }q^0,\delta q^1=q^0+\frac{p_1}{\kappa }q^1.$$
(18)
For a free particle, with Hamiltonian (14), the Hamilton equations derived from (16) are

$`\dot{q}^0`$ $`=`$ $`\left(1{\displaystyle \frac{p_0}{\kappa }}\right){\displaystyle \frac{H}{p_0}}{\displaystyle \frac{p_1}{\kappa }}{\displaystyle \frac{H}{p_1}}={\displaystyle \frac{p_0}{\left(1\frac{p_0}{\kappa }\right)^2}},`$ (19)
$`\dot{q}^1`$ $`=`$ $`{\displaystyle \frac{H}{p_1}}={\displaystyle \frac{p_1}{\left(1\frac{p_0}{\kappa }\right)^2}},`$ (20)
$`\dot{p}_0`$ $`=`$ $`\left(1{\displaystyle \frac{p_0}{\kappa }}\right){\displaystyle \frac{H}{q^0}}=0,`$ (21)
$`\dot{p}_1`$ $`=`$ $`{\displaystyle \frac{H}{q^1}}+{\displaystyle \frac{p_1}{\kappa }}{\displaystyle \frac{H}{q^0}}=0,`$
and hence $`v=p_1/p_0`$, as expected.

It is known that the Hamilton equations for systems with nonstandard symplectic structure can be derived from an action principle . Given a phase space with symplectic structure $`\{Q^A,Q^B\}=ø^{AB}`$, where $`Q^A`$ denote either the coordinates or the momenta, one defines the functions $`R_A(Q^A)`$ such that

$$\frac{R_A}{Q^B}\frac{R_B}{Q^A}=ø_{AB},$$
(22)
where $`ø_{AB}`$ is the inverse of $`ø^{AB}`$. The Hamilton equations can then be obtained varying with respect to $`Q^A`$ the action

$$I=𝑑\tau (R_A\dot{Q}^AH).$$
(23)
Note that in general the action so defined contains derivatives of the momenta.

In our case, this procedure yields

$$I=𝑑\tau \left[\kappa \mathrm{log}\left(1\frac{p_0}{\kappa }\right)\dot{q}^0+p_1\dot{q}^1\frac{p_1q^1}{\kappa \left(1\frac{p_0}{\kappa }\right)}\dot{p}_0H\right].$$
(24)
After integration by parts, (24) can also be written as

$$I=𝑑\tau \left[\frac{q^0+\frac{p_1}{\kappa }q^1}{\left(1\frac{p_0}{\kappa }\right)}\dot{p}_0+q^1\dot{p}_1+H\right].$$
(25)
The transformation laws for the coordinates can now be obtained by requiring that the variables $`r^a`$ conjugate to the momenta $`p_a`$ in (25) transform covariantly, i.e., according to $`r^ar^a=Ł_b^ar^b`$. From (25) one has

$$q^0=\left(1\frac{p_0}{\kappa }\right)r^0\frac{p_1}{\kappa }r^1,q^1=r^1,$$
(26)
and hence

$$q^0=\left(1\frac{W_0}{\kappa }\right)Ł_a^0r^a\frac{W_1}{\kappa }Ł_a^1r^a,q^1=Ł_a^1r^a.$$

Substituting (11) and (12), after tedious but elementary algebraic manipulations, one gets the surprisingly simple result

$`q^0=\mathrm{\Delta }(q^0\mathrm{cosh}\xi q^1\mathrm{sinh}\xi ),`$ (27)
$`q^1=\mathrm{\Delta }(q^0\mathrm{sinh}\xi +q^1\mathrm{cosh}\xi ).`$ (28)
Thus, the transformation laws of coordinates are identical to the usual Lorentz transformations, except for the momentum-dependent factor $`\mathrm{\Delta }`$. It is then easy to build a simple “invariant length”

$$(p_ap^a)q^aq_a\left(1\frac{p_0}{\kappa }\right)^2q^aq_a,$$
(29)
from which an invariant line element can be defined. One may also define new (commuting) coordinates $`\overline{q}^a=(1p_0/\kappa )q^a`$, which transform according to standard Lorentz transformations, so that $`\overline{q}^a\overline{q}_a`$ is invariant. These coordinates however do not lead to a correct definition of the velocity of a particle.

## 4 Conclusions

We have discussed the coordinate transformation laws that allow a covariant definition of the Hamiltonian formalism in models of deformed special relativity, both in the case of commuting coordinates with canonical symplectic structure, and of noncommuting coordinates with deformed Poisson brackets. The transformation laws for coordinates are momentum dependent and in general rather complicated, except in the case of the MS model with noncommuting coordinates. In this case the transformation laws assume a natural form, that allows the definition of a simple (momentum-dependent) line element.

Of course, the interpretation of the momentum dependence of the transformation laws for coordinates is highly nontrivial. The deformed Lorentz transformations act on the full phase space and not separately on coordinate and momentum space. This implies that the coincidence of two events becomes observer-dependent. At first, this prediction may seem unphysical, but it must be considered that the theory is assumed to be effective at sub-Planckian scales, of which we do not have any direct experience.

A more conservative interpretation of the results of this paper would be to consider the $`q^a`$ simply as labels of the position of a particle, which should be connected to the spacetime coordinates $`x^\mu `$ by a sort of vierbein field: $`q^a=e_\mu ^ax^\mu `$. The consistency of this approach is currently being investigated.

These considerations are also interesting in view of the inclusion of gravity in the theory . A consistent approach should presumably lead to some kind of phase space extension of general relativity.

Note. While completing this work, I became aware of a paper  where the same transformations of coordinates (27) are proposed, starting from a different point of view.

## 5 Appendix

We report here the calculations analogous to those presented above, for the case of the LNR model .

The infinitesimal Lorentz transformations for LNR are given by (2), with

$$w_0=p_1,w_1=\frac{\kappa }{2}\left(1\mathrm{e}^{2p_0/\kappa }\frac{p_1^2}{\kappa ^2}\right).$$
(30)
It follows from (6) that

$$\delta q^0=\mathrm{e}^{2p_0/\kappa }q^1,\delta q^1=q^0\frac{p_1}{\kappa }q^1.$$
(31)
The finite transformations are given by (8), with

$$W_0=p_0+\kappa \mathrm{log}\mathrm{\Gamma },W_1=\frac{p_1\mathrm{cosh}\xi +\frac{\kappa }{2}\left(1\mathrm{e}^{2p_0/\kappa }+\frac{p_1^2}{\kappa ^2}\right)\mathrm{sinh}\xi }{\mathrm{\Gamma }},$$
(32)
where

$$\mathrm{\Gamma }=\frac{1}{2}\left(1+\mathrm{e}^{2p_0/\kappa }\frac{p_1^2}{\kappa ^2}\right)+\frac{1}{2}\left(1\mathrm{e}^{2p_0/\kappa }+\frac{p_1^2}{\kappa ^2}\right)\mathrm{cosh}\xi +\frac{p_1}{\kappa }\mathrm{sinh}\xi .$$

From (10) follows that $`q^a=Ł_b^aq^b`$, with

$`Ł_{1}^1`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1\mathrm{e}^{2p_0/\kappa }{\displaystyle \frac{p_1^2}{\kappa ^2}}\right)+{\displaystyle \frac{1}{2}}\left(1+\mathrm{e}^{2p_0/\kappa }+{\displaystyle \frac{p_1^2}{\kappa ^2}}\right)\mathrm{cosh}\xi +{\displaystyle \frac{p_1}{\kappa }}\mathrm{sinh}\xi ,`$ (33)
$`Ł_{0}^1`$ $`=`$ $`\mathrm{e}^{2p_0/\kappa }\left[{\displaystyle \frac{p_1}{\kappa }}\left(1\mathrm{cosh}\xi \right)\mathrm{sinh}\xi \right],`$ (34)
and $`Ł_{0}^0=Ł_{1}^1/\mathrm{\Gamma }`$, $`Ł_{0}^1=Ł_{1}^0/\mathrm{\Gamma }`$.

Also in the LNR case the standard hamiltonian formalism does not give a consistent definition of the particle velocity. A suitable definition of velocity is given instead by the right-invariant velocity ,

$$v=\frac{\mathrm{e}^{p_0/\kappa }p_1}{\kappa \mathrm{sinh}\frac{p_0}{\kappa }+\mathrm{e}^{p_0/\kappa }\frac{p_1^2}{\kappa }},$$
(35)
and can be obtained introducing the following nonstandard symplectic structure :

$`\{q^0,q^1\}={\displaystyle \frac{q^1}{\kappa }},\{p_0,p_1\}=0,\{q^0,p_0\}=1,`$ (36)
$`\{q^1,p_1\}=1,\{q^0,p_1\}={\displaystyle \frac{p_1}{\kappa }},\{q^1,p_0\}=0.`$ (37)
The infinitesimal transformations of coordinates can then be obtained as before from the Jacobi identities. They read

$$\delta q^0=\frac{\kappa }{2}\left(1\mathrm{e}^{2p_0/\kappa }\frac{p_1^2}{\kappa ^2}\right)q^1+\frac{p_1}{\kappa }q^0,\delta q^1=q^0.$$
(38)
The Hamiltonian for a free particle is given by the Casimir invariant

$$H=\frac{1}{2}\left[\left(2\kappa \mathrm{sinh}\frac{p_0}{2\kappa }\right)^2\mathrm{e}^{p_0/\kappa }\frac{p_1^2}{\kappa ^2}\right],$$
(39)
and the Hamilton equations are

$`\dot{q}^0`$ $`=`$ $`{\displaystyle \frac{H}{p_0}}{\displaystyle \frac{p_1}{\kappa }}{\displaystyle \frac{H}{p_1}}=\kappa \mathrm{sinh}{\displaystyle \frac{p_0}{\kappa }}+\mathrm{e}^{p_0/\kappa }{\displaystyle \frac{p_1^2}{\kappa }},`$ (40)
$`\dot{q}^1`$ $`=`$ $`{\displaystyle \frac{H}{p_1}}=\mathrm{e}^{p_0/\kappa }p_1,`$ (41)
$`\dot{p}_0`$ $`=`$ $`{\displaystyle \frac{H}{q^0}}=0,`$ (42)
$`\dot{p}_1`$ $`=`$ $`{\displaystyle \frac{H}{q^1}}+{\displaystyle \frac{p_1}{\kappa }}{\displaystyle \frac{H}{q^0}}=0,`$ (43)
which lead to the definition (35) of velocity.

As explained above, eqs. (40) can be obtained from an action principle. The action reads in this case,

$$I=𝑑\tau \left[p_0\dot{q}^0+p_1\dot{q}^1\frac{q^1p_1}{\kappa }\dot{p}_0H\right],$$
(44)
or, after integration by parts,

$$I=𝑑\tau \left[\left(q^0+\frac{p_1}{\kappa }q^1\right)\dot{p}_0+q^1\dot{p}_1+H\right].$$
(45)
The transformation laws for the coordinates can now be obtained by requiring that the variables $`r^a`$ conjugate to the momenta $`p_a`$ in (45) transform covariantly. One has

$$q^0=r^0\frac{p_1}{\kappa }r^1,q^1=r^1,$$
(46)
and hence

$$q^0=Ł_a^0r^a\frac{1}{\kappa }W_1Ł_a^1r^a,q^1=Ł_a^1r^a.$$

After some algebra, one gets the result

$`q^0`$ $`=`$ $`q^0(\mathrm{cosh}\xi +{\displaystyle \frac{p_1}{\kappa }}\mathrm{sinh}\xi ){\displaystyle \frac{q^1}{2}}\left(1+\mathrm{e}^{2p_0/\kappa }{\displaystyle \frac{p_1^2}{\kappa ^2}}\right)\mathrm{sinh}\xi ,`$ (47)
$`q^1`$ $`=`$ $`q^0\left[\mathrm{sinh}\xi +{\displaystyle \frac{p_1}{\kappa }}(1\mathrm{cosh}\xi )\right]`$
$`+{\displaystyle \frac{q^1}{2}}\left[\left(1\mathrm{e}^{2p_0/\kappa }+{\displaystyle \frac{p_1^2}{\kappa ^2}}\right)+\left(1+\mathrm{e}^{2p_0/\kappa }{\displaystyle \frac{p_1^2}{\kappa ^2}}\right)\mathrm{cosh}\xi \right].`$
In this case, the transformation laws for coordinates are not especially simple and do not seem to lead to the same interesting developments as in the MS case.


text:
# Planck-Kerr Turbulence

## 1 Introduction

Observations show that the statistical distribution of galaxies and their clusters are chaotic and homogeneous, suggesting that something about the big bang mechanism was capable of randomly producing and randomly distributing the energy density seeds required to trigger nucleosynthesis and gravitational structure formation in patterns characteristic of high Reynolds number turbulence. Was the big bang turbulent? The big bang is presently the no-man’s-land of modern physics. Quantum mechanics (QM) and general relativity (GR) theories develop singularities under the big bang Planck scale conditions that have not yet been reconciled by multidimensional superstring theory or any other theory of the quantum gravitational dynamics (QGD) epoch. The string tension in string theory $`c^4G^1=1.1\times 10^{44}`$ kg m s-2 (Greene, 1999) appears in Einstein’s equation and matches the Planck force $`F_P`$ (see Table 1).

Turbulence theory is not in much better condition than QM and GR. No definition of turbulence is generally accepted, and it is widely believed that turbulence always cascades from large scales to small, despite clear evidence to the contrary in jets, wakes and boundary layers where turbulence length scales monotonically increase. How can turbulence exist at the time of the big bang when there are no large scales of anything to supply the energy? Furthermore, observations of cosmic microwave background (CMB) temperature anisotropies $`\delta T/T10^5`$ prove the plasma before transition to gas was not turbulent, since turbulence would produce much larger $`\delta T/T10^2`$ values (Gibson, 2000). Whatever turbulence emerged from the big bang was damped by processes in the plasma epoch.

In the following we apply the Planck-Kerr instability model to the QGD epoch before inflation, using dimensional analysis. It appears that not only did QGD turbulence exist, but that turbulence was the key mechanism responsible for the big bang. Evidence of high Reynolds number turbulence preserved by the CMB is discussed, and conclusions are presented.

## 2 Theory

General relativity theory requires the length scale $`L_{SS}`$ of black holes of mass $`m`$ to approximate the Schwarzschild radius $`L_{SS}Gm/c^2`$, where $`G=6.67\times 10^{11}`$ m3 kg-1 s-2 is Newton’s gravitational constant and $`c=3\times 10^8`$ m1 s-1 is the speed of light. Quantum mechanical wave-particle duality assigns a de Broglie wavepacket size $`L_{dB}h/mc`$ to particles of mass $`m`$, where $`h1.05\times 10^{34}`$ kg m2 s-1 is Planck’s constant. The GR and QM length scales match when the two theories break down, where $`L_{SS}L_{dB}`$ and the black hole mass and particle mass equal the Planck mass $`m_P[chG^1]^{1/2}=2.12\times 10^8`$ kg. Substituting $`m_P`$ gives the Planck length $`L_P[c^3hG]^{1/2}=1.62\times 10^{35}`$ m. The Hawking evaporation time of a black hole and its de Broglie wave propagation time $`c^1L_P`$ match at the Planck time $`t_P[c^5hG]^{1/2}=5.41\times 10^{44}`$ s. The Hawking temperature $`T_H[c^3hG^1k^1m^1]`$ of a Planck scale black hole (Planck particle) is $`T_H=[c^5hG^1k^2]^{1/2}T_P`$, where $`T_P=1.40\times 10^{32}`$ K is defined as the Planck temperature and $`k=1.38\times 10^{23}`$ kg m2 s-2 K-1 is the Boltzmann constant. The energy of Planck particles $`E_P[c^5hG^1]^{1/2}=1.94\times 10^9`$ kg m2 s-2. The Planck entropy $`S_PE_P/T_P=k`$.

How could the quantum gravitational dynamics epoch be turbulent while we know the plasma epoch at the CMB time of photon decoupling was not? Turbulence is an eddy-like state of fluid motion where the inertial vortex forces of the eddies are larger than any other forces that might tend to damp the eddies out (Gibson, 1991). Large Planck inertial vortex forces per unit mass $`[\stackrel{}{v}\times \stackrel{}{\omega }]_P[c^7h^1G^1]^{1/2}g_P=5.7\times 10^{51}`$ m s-2 must match and sometimes overwhelm the Planck gravitational acceleration $`g_P`$ as well as the Planck viscous forces for big bang turbulence to be triggered and persist. The size of the universe at the strong force freeze out time $`t_{SF}10^{35}`$ s is $`L_{SF}ct_{SF}=3\times 10^{27}`$ m at the strong force freeze out temperature $`T_{SF}10^{28}`$ K (Peacock, 2000). The strong force velocity of the Planck particles $`V_{SF}c[T_{SF}/T_P]^{1/2}=3\times 10^6`$ m s-1. The Planck viscosity $`\nu _P=cL_P=[c^1hG]^{1/2}=4.8\times 10^{27}`$ m2 s-1, assuming viscous stresses are transported by Planck particles at light speeds with Planck length scale collision distances. Thus, Reynolds numbers $`VL/\nu `$ in the QGD epoch increase to $`Re_{SF}V_{SF}\times L_{SF}/\nu _{SF}10^6`$ values that are well above critical. Small, weakly collisional particles like gluons, neutrinos and photons that inhibit turbulence by viscous stresses after the strong force freeze out cannot exist at QGD temperatures.

The Planck-Kerr mechanism for big bang turbulence formation is shown in Figure 1. Vacuum oscillations at Planck length scales allow a small possiblity for Planck particles and Planck anti-particles to appear spontaneously by quantum tunneling. Once they appear they increase the possibility of more Planck particles to be triggered by their enormous temperatures. Stable spinning quantum states analogous to positronium are possible (extreme Kerr black holes). Prograde accretions of Planck particles on these objects may trigger a turbulent Planck gas with vorticity matching the Plank-Kerr spin, taking advantage of the increased prograde energy release available from small marginally stable prograde orbits (Peacock 2000 p61) to increase the probability of further Planck-Kerr and Planck particle formations. Energies up to 42% $`E_P`$ are available.

The Planck-Kerr instability of Figure 1 has Planck power $`P_Pc^5G^1=3.64\times 10^{52}`$ kg m2 s-2, which is $`10^4`$ larger than the power of all galaxies in our horizon at the present time. The Planck turbulence stress per unit mass $`[\stackrel{}{v}\times \stackrel{}{\omega }]_P/L_{P}^{}{}_{}{}^{2}=\sigma _P[c^{13}h^3G^3]^{1/2}=1.26\times 10^{121}`$ m-1 s-2. From general relativity theory, space-time is created by sufficiently large negative stresses. Planck-Kerr turbulence supplies the power, stresses and entropy required to initiate space-time creation and reach the strong force phase transition state where negative stresses of the false vacuum can trigger inflation (Guth, 1997).

## 3 Observations

Fingerprints of high Reynolds number turbulence have been observed in the CMB temperature anisotropies using extended self similarity (ESS) coefficients (Bershadskii and Sreenivasan, 2002). Comparisons with high Reynolds number turbulence ESS coefficients and with ESS coefficients of other flows (Benzi et al., 1996) are shown in Figure 2 (Gibson, 2001). The CMB and turbulence ESS coefficients are in nearly perfect agreement with each other, but do not match ESS coefficients computed from non-turbulent flows.

## 4 Conclusions

A variety of indicators suggest that the big bang was not only turbulent but strongly turbulent before inflation, even though that time period was short and the length scales small. No strong sources of irreversibility besides turbulence and turbulent mixing exist at Planck scales where the temperatures are high and highly reversible vacuum oscillations produce little or no entropy. The Planck specific entropy $`s_P[c^1h^1Gk^2]^{1/2}`$ was only $`6.35\times 10^{16}`$ m2 s-2 K-1. The turbulent viscous dissipation rate $`\epsilon _P[c^9h^1G^1]^{1/2}`$ was $`1.72\times 10^{60}`$ m2 s-3.

The Planck-Kerr instability is identified between extreme Schwarzschild black holes and extreme Kerr black holes that can produce small scale vorticity and an expanding universe driven by Planck inertial vortex forces and turbulence stresses that overcome Planck gravity, stetch space, and create more turbulent Planck gas. Strong big bang turbulence mixes the resulting temperature fluctuations as the universe cools to the strong force phase transition temperature. Inflation fossilizes the turbulent temperature fluctuations by stretching the length scales beyond the scale of causal connection.

Fossil big bang temperature turbulence fluctuations retain their turbulence signatures by seeding nucleosynthesis, and are detected by observations of the homogenous chaotic pattern of large scale cosmic structures, and the extended self similarity coefficients and dissipation multiscaling parameters of the cosmic microwave background temperature anisotropies that identically match those of high Reynolds number turbulence and mixing.


text:
# Calculation of muon transfer from muonic hydrogen to atomic oxygen

## I Introduction

Negative muon transfer between exotic atoms (muonic hydrogen, for instance) and other atoms or molecules, has been extensively studied in the framework of muon catalysed nuclear fusion (see Ref. Ponomarev (2001) and literature cited therein). Also, the structural and spectroscopic properties of these species is of interest for metrology as well as a test of quantum electrodynamic theories Bakalov et al. (1993).

Very recently, a new method to measure the hyperfine structure ($`F=0,1`$) of muonic hydrogen ($`p\mu `$) based on the collisional energy dependence of muon transfer from the muonic hydrogen to an oxygen molecule, has been proposed Adamczak et al. (2001). When muonic hydrogen in the $`F=0`$ hyperfine ground state is laser excited to the $`F=1`$ level, collisions with $`H_2`$ convert this excess energy into kinetic energy, giving an additional 0.12 eV translational energy to the muonic hydrogen. When the muon is transfered to an oxygen atom, it is captured in high ($`n=5,6`$) states which promptly de-excite and emit X-rays. If the muon transfer rate to oxygen changes significantly from thermal (0.04 eV) to epithermal (0.16 eV) energies, the measurement of the X-ray emission intensity with and without laser excitation can be used for the determination of the hyperfine splitting in muonic hydrogen Bakalov et al. (1993); Adamczak et al. (2001). This proposal was based on the work of Werthmüller et al Werthmüller et al. (1998) which indicates that muon transfer from hydrogen to oxygen increases by a factor of 3-4 going from thermal to epithermal collision energies.

While the experimental transfer rate at thermal energies have been borne out by the calculations of Sultanov and Adhikari Sultanov and Adhikari (2000), the rates at epithermal energies have never been confirmed theoretically. Even more, such an increase of the transfer rate with energy does not have a clear theoretical explanation. Werthmüller et al Werthmüller et al. (1998) suggest the existence of a resonance at low energies. It is worth to stress that such an increase has not been observed in the case of muonic hydrogen colliding with sulfur Mulhauser and Schneuwly (1993) while for the case of CH4, experimental data suggest that the transfer rate actually decreases going from thermal to epithermal energies Kirch et al. (1998). In order to improve our understanding of this process and to assess the validity of the proposal of Ref. Adamczak et al. (2001), it is important to perform model calculations to study in detail this reaction as a function of energy. It is the purpose of this paper to present such a study.

Since the muonic hydrogen has to approach one of the oxygen nuclei very close in order for the muon to transfer Gershtein (1963), the process can be described as

$$(p\mu )+\text{O}^{8+}p+(\mu \text{O})^{7+}$$
(1)
Although there have been several full three-dimensional calculations of muon transfer rates at low energies between muonic-hydrogen and low-$`Z`$ atoms (see literature cited in Ref. Sultanov and Adhikari (2002)), there is none when the transfer involves nuclei with $`Z>3`$ . There are several good reasons for that. As $`Z`$ increases ($`Z=8`$ in our case) there is a larger initial-channel polarization and a stronger final-channels Coulomb interaction. In addition the number of open channels even at zero relative kinetic energy increases with $`Z`$ making the full three-dimensional calculation very heavy. Thus up-to-now only approximate calculations have been performed for the muon transfer rate between muonic-hydrogen and oxygen Haff et al. (1977); Gershtein (1963); Sultanov and Adhikari (2000). Among them, the most recent and accurate is the one done by Sultanov and Adhikari Sultanov and Adhikari (2000), who used a two-channel approximation to the integro-differential Fadeev-Hahn formalism to calculate the muon transfer rates from muonic-hydrogen to the $`n=5`$ states of oxygen.

We have chosen to perform exact calculations but in a restricted configuration with the muon moving along the line joining the proton to the oxygen atom. Indeed, for this system the colinear configuration is the most favorable for the transfer process. Furthermore, since the colinear configuration provides the minimum energy path for the reactance channel, even for an initial non colinear configuration there will be efficient orientational effects in particular for low translational energies.

Two model interaction potentials have been used in the calculations:

1. A pure coulombic potential with the bare oxygen nucleus $`\text{O}^{8+}`$.
2. A shorter range model potential with a distance dependent effective charge on the oxygen nucleus (calculated via a Thomas-Fermi electronic density model) to simulate the effect of the outer electrons.

The Hamiltonian was written in terms of hyperspherical coordinates. A piecewise diabatic basis set on the hyperspherical angle was used to expand the wave function. The basis functions have in turn been expanded in terms of first order Legendre functions to efficiently handle the two coulombic singularities at fixed hyperradii. The resulting close-coupling time-independent Schrödinger equations in the hyperradius were solved using a de Vogelaere algorithm and the partial and total muon transfer probabilities were determined by the standard S-matrix analysis at large distances. Since for energies below $`10^1`$ eV the muon transfer process studied here is equivalent to an ultra-cold collision (de Broglie wavelength ($`>1`$ Å) much larger than the characteristic distance ($`0.1`$ Å) of the potential interaction), special care had to be taken to the asymptotic analysis in the reactant channel.

The paper is organized as follows. Section 2 introduces the model and the methodology used in the calculations. Section 3 presents the calculated muon transfer probabilities together with their interpretation in terms of simple Landau-Zener and threshold models. A procedure to estimate their energy behavior in the 3-dimensional case, is also presented. In section 4, the muon transfer rates are discussed in the light of previous model calculations and available experimental data for this process. Finally, section 5 is devoted to the conclusions.

## II Methodology

The mass-scaled Jacobi coordinates adapted to the entrance channel of reaction (1) are :

$`R`$ $`=`$ $`\sqrt{{\displaystyle \frac{m_{\text{O},p\mu }}{m}}}\left(x_\text{O}{\displaystyle \frac{m_\mu x_\mu +m_px_p}{m_p+m_\mu }}\right);`$ (2)
$`r`$ $`=`$ $`\sqrt{{\displaystyle \frac{m_{p,\mu }}{m}}}(x_\mu x_p)`$ (3)
where

$`m_{\text{O},p\mu }`$ $`=`$ $`{\displaystyle \frac{m_\text{O}(m_p+m_\mu )}{m_\text{O}+m_p+m_\mu }};m_{p,\mu }={\displaystyle \frac{m_pm_\mu }{m_p+m_\mu }};`$ (4)
$`m`$ $`=`$ $`\left({\displaystyle \frac{m_\text{O}m_pm_\mu }{m_\text{O}+m_p+m_\mu }}\right)^{1/2}`$ (5)
The hyperspherical coordinates are then defined by

$$\rho =\sqrt{R^2+r^2};\mathrm{tan}\theta =r/R$$
(6)
where $`0\theta \theta _\mu `$ , with $`\theta _\mu =\text{arctan}(m_\mu /m)`$. In terms of these coordinates and after regularisation of the wave function by the factor $`\sqrt{\rho }`$ (the volume element is then given by $`d\rho d\theta `$), the time-independent Schrödinger equation at total energy $`E`$ is:

$`{\displaystyle \frac{\mathrm{}^2}{2m}}\left({\displaystyle \frac{^2}{\rho ^2}}+{\displaystyle \frac{1}{\rho ^2}}{\displaystyle \frac{^2}{\theta ^2}}\right)\psi (\rho ,\theta )`$
$`+\left(V(\rho ,\theta ){\displaystyle \frac{\mathrm{}^2}{8m\rho ^2}}E\right)\psi (\rho ,\theta )=0`$ (7)
where $`V(\rho ,\theta )`$ is the interaction potential given by

$`V={\displaystyle \frac{e^2Z_\text{O}^{}(|x_\mu x_\text{O}|)}{|x_\mu x_\text{O}|}}`$
$`+{\displaystyle \frac{e^2Z_\text{O}^{}(|x_px_\text{O}|)}{|x_px_\text{O}|}}{\displaystyle \frac{e^2}{|x_px_\mu |}}`$ (8)
Two sets of calculations have been performed. One with the bare oxygen nucleus, for which $`Z_\text{O}^{}=8`$ in (8), and another with $`Z_\text{O}^{}(d)`$ being an effective charge given by the Thomas-Fermi model Landau and Lifshitz (1977) and represented in Fig. 1.

The total wave function $`\psi (\rho ,\theta )`$ is expanded in terms of basis set wave functions depending on the hyperspherical angle $`\theta `$. We use a diabatic-by-sector representation. In each sector $`\rho _n\delta \rho _n\rho <\rho _n+\delta \rho _n;n=1,\mathrm{},N_\rho `$ we write:

$$\psi (\rho ,\theta )=\underset{i}{}F_i(\rho )\varphi _i(\theta ;\rho _n)$$
(9)
where $`\varphi _i(\theta ;\rho _n)`$ are eigenstates of the Hamiltonian at fixed $`\rho _n`$ distances. Their calculation requires the solution of a bound state problem for a potential showing Coulomb singularities at $`\theta =0`$ and $`\theta =\theta _\mu `$. Treatment of Coulomb singularities requires the use of specific methods such as increasing the density of grid points, or alternatively the oscillation frequency of the basis functions, near the singularities. Mapped Fourier Fattal et al. (1996); Lemoine (2000), Lagrange Baye and Vincke (1999) or Schwartz Schwartz (1985); Dunseath et al. (2002) interpolations are examples of such methods which have been used for cases with one singularity.

In our case we have two singularities to deal with at fixed hyperradius. Therefore we have used the auxilliary coordinate $`x=2\theta /\theta _\mu 1`$. This new variable is equivalent to one of the two hyperspherical elliptic angles Tolstikhin et al. (1995), but constrained to be in the range \[-1,1\] with the singularities being at the boundaries. If the $`\varphi _i(\theta ;\rho _n)`$ functions are renormalized as: $`\varphi _i(\theta ;\rho )=(1x^2)^{\frac{1}{2}}\overline{\varphi _i}(x;\rho )`$, both singularities of the potential are regularized. The $`\overline{\varphi _i}(x;\rho )`$ functions obey the differential equation:

$`\left({\displaystyle \frac{2\mathrm{}^2}{m\theta _{\mu }^{}{}_{}{}^{2}\rho ^2}}\widehat{D}+(1x^2)(V{\displaystyle \frac{\mathrm{}^2}{8m\rho ^2}})\right)\overline{\varphi _i}(x;\rho )`$
$`=(1x^2)ϵ_i(\rho )\overline{\varphi _i}(x;\rho )`$ (10)
where

$$\widehat{D}=(1x^2)\frac{^2}{x^2}2x\frac{}{x}\frac{1}{1x^2}$$
(11)
Eq. 10 is solved by expanding $`\overline{\varphi _i}(x;\rho )`$ functions on the basis set of the eigenvectors of $`\widehat{D}`$, which are the associated Legendre functions $`P_n^1(x)`$. This yields the generalized eigenvalue problem :

$`\left({\displaystyle \frac{2\mathrm{}^2}{m\theta _\mu ^2\rho ^2}}𝑫+𝑾\right)\overline{\mathit{\varphi }}_i(\rho )`$
$`=ϵ_i(\rho )𝑶\overline{\mathit{\varphi }}_i(\rho )`$ (12)
where $`\overline{\mathit{\varphi }}_i(\rho )`$ is the vector of the unknown coefficients of the function $`\overline{\varphi _i}(x;\rho )`$ in the Legendre basis set. The matrix representation of the differential operator $`𝑫`$ is diagonal, with diagonal elements $`n(n+1)`$. $`𝑾`$ is the potential coupling matrix and $`𝑶`$ is the representation of $`(1x^2)`$ in the Legendre basis. These matrices are calculated from a transformation of their diagonal representation in a grid of Gauss-Legendre quadrature points Corey and Lemoine (1992). The system (12) is transformed to a standard eigenvalue problem after left multiplication by $`𝑶^{\frac{1}{2}}`$. The method is expected to converge fast with the size of the basis set (or equivalently with the number of grid points) since the use of a Gauss-Legendre scheme provides an adequately high density of points near $`x=\pm 1`$. Since the wavefunctions concentrate in the vicinity of $`x=\pm 1`$ as energy decreases or as $`\rho `$ increases, we expect convergence to be most difficult for lower levels and large $`\rho `$. At large $`\rho `$, we also expect convergence to be more difficult to achieve for $`O\mu `$ states which are more compact than $`p\mu `$ states. We have computed the relative error for a given number $`n_L`$ of Legendre basis functions, using $`n_L=450`$ as a reference. For $`n_L=150`$, we obtained 86 states converged better than $`10^8`$ of relative error up to $`\rho a_0`$. For larger $`\rho `$ values, the convergence starts to deteriorate, first for the ground state $`O\mu (n=1)`$, then for the first excited state $`n=2`$ near $`\rho =2a_0`$…Using larger values of $`n_L`$ delays the limit of convergence deterioration. For instance, with $`n_L=250`$, the ground state is converged with $`10^8`$ relative precision up to $`\rho 5a_0`$, up to $`\rho 10a_0`$ with $`n_L=350`$. This latter value of $`n_L`$ was used to generate the results shown below, and provided a basis converged better than $`10^4`$ up to $`\rho 40a_0`$ for all states except the first four ones which do not play any significant role in the charge transfer process considered here. As expected, among all these states, the ones corresponding to bound states of $`p\mu `$ are the easiest to converge : the relative error at $`\rho =40a_0`$ for $`p\mu `$(n=1) is $`10^{10}`$ with $`n_L=250`$ and $`10^{12}`$ with $`n_L=350`$. High accuracy on the $`p\mu (n=1)`$ channel is especially important for our present study since we consider initial relative kinetic energies as low as $`10^6`$ eV.

The close-coupling equations are integrated along the hyperradius $`\rho `$ using the de Vogelaere algorithm Lepetit et al. (1986). In order to have convergence of the transfer probabilities to better than 1 % in the energy range considered here, we included 29 channels: ($`(p\mu )_{n=14}+\text{O}`$ and $`p+(\mu \text{O})_{n=125}`$). The integration of the coupled equations was performed from the origin to $`\rho _{\text{end}}30a_0`$.

The asymptotic analysis has been performed using the appropriate Jacobi coordinates for the entrance and for the product channels. The asymptotic states were expressed as products of Coulombic bound states and translational functions. For the product channels, the translational functions were approximated by simple plane waves. However, since the initial relative translational energy can be very small, the entrance channel has to be described with care to account for residual long range interactions. For the pure coulombic potential, this residual interaction has a $`1/R^2`$ dependence (charge-dipole interaction). In this case, Bessel functions of imaginary order have been used as a basis set in $`R`$. For the screened Thomas-Fermi potential, wave functions depending on $`R`$ have been numerically computed by inward integration from very large distances. These asymptotic wave functions and their derivatives are projected onto the hyperspherical basis functions $`\varphi _i(x;\rho _{\text{end}})`$ Lepetit et al. (1986). Transfer probabilities are then obtained by a standard S-matrix analysis.

## III Results of the calculations

We have performed calculations for collision energies in the range $`10^610^3`$ eV. In Fig. 2 we present the adiabatic energies $`ϵ_i`$ (see Eq. 10) as a function of the hyperradius $`\rho `$. The origin of energies has been chosen to coincide with the asymptotic limit of the entrance channel $`(p\mu )_{n=1}+\text{O}`$. Thus the calculations cover the energy range between this limit and the $`p+(\mu \text{O})_{n=10}`$ threshold.

The dynamics of muon transfer can be qualitatively understood by inspection of this figure. Starting in channel $`(p\mu )_{n=1}+\text{O}`$, the system crosses diabatically the channel $`p+(\mu \text{O})_{n=7}`$. Muon transfer is completely negligible as the coupling is very small compared with the collision energy. The couplings to channels $`p+(\mu \text{O})_{n=6}`$ and $`p+(\mu \text{O})_{n=5}`$ are larger as evidenced by avoided crossings. The other channels $`p+(\mu \text{O})_{n<5}`$ are weakly coupled to the initial one and they are not expected to be populated significantly.

In Figs. 3 and 4 we present the muon transfer probabilities into the different product channels $`p+(\mu \text{O})_{n=5,6}`$ together with the total transfer probability, for the pure coulombic (C) and the Thomas-Fermi (TF) potentials, respectively. The two channels are about equally populated in both cases and the total transfer probability is high (80 %) for energies in the intermediate range (between $`10^2`$ and $`10^2`$ eV). This can be understood in terms of non-adiabatic transitions in a simple 3 channels model ( $`(p\mu )_{n=1}+\text{O}`$ and $`p+(\mu \text{O})_{n=6,5}`$). Using the Landau-Zener expressions for the two crossings, a total transfer probability of $`0.6`$ between 0 and $`10^2`$ eV and a slow decrease at higher energies is obtained. This is in qualitative agreement with the fully converged results.

The main difference between the results for the C and TF potentials is at low energies. Whereas for the C potential the probabilities are constant at low energies, for the TF potential there is a typical threshold behavior. This is what is expected for short range potentials Landau and Lifshitz (1977). Indeed, in the colinear configuration the C potential behaves asymptotically as $`1/R^2`$ (charge-dipole interaction) while the TF potential goes to zero much faster. At low energies, the amplitude of the wave function for $`1/R^2`$ potentials near the origin is constant, while for short range potentials it decreases as energy goes to zero Landau and Lifshitz (1977). Thus the low energy dependence of the transfer probabilities is enterily determined by the asymptotic behavior of the potentials. This suggests that the dynamics can be described by a two step mechanism: a transmission in the incident channel through the long range part of the potential, followed by non adiabatic transitions in the interaction region. The charge transfer probability for low energies will be then given by

$$P(E)=|T(E)|^2P_{max}$$
(13)
where $`T(E)`$ is the transmittance in the incoming channel and $`P_{max}`$ the plateau transfer probability.

The transmittance $`T(E)`$ can be estimated by a simple one channel calculation with an effective potential of the form

$$V_{\text{eff}}(R)=\frac{\alpha Z_\text{O}^{}(R)}{R^2};\alpha =\frac{3\mathrm{}^2}{2m_{p,\mu }}$$
(14)
for $`R>R_0`$ and $`V_{\text{eff}}(R)=V_{\text{eff}}(R_0)`$ for $`R<R_0`$. The transfer radius has been taken at the avoided crossing with the $`p+(\mu \text{O})_{n=6}`$ channel. For $`R<R_0`$ the solution of this one-channel problem is given by

$$F(R)=T(E)\frac{\mathrm{sin}(KR)}{K^{1/2}}$$
(15)
with $`K=(2m_{\text{O},p\mu }(EV_{\text{eff}}(R_0)))^{1/2}/\mathrm{}`$, while for $`R\mathrm{}`$ we have $`F(R)=\mathrm{sin}(kR+\delta )/k^{1/2}`$ with $`k=(2m_{\text{O},p\mu }E)^{1/2}/\mathrm{}`$. Numerical solution of the one channel problem thus provides directly $`T(E)`$. For the pure Coulombic potential ($`Z_\text{O}^{}=8`$) we found $`T(E)=1`$, consistent with the result of our multichannel calculations which gives a constant transfer probability at threshold (see Fig. 3). This unusual behavior is due to the particular $`1/R^2`$ dependence of the effective potential Landau and Lifshitz (1977). The shorter range TF potential gives a more standard threshold behavior, with a transmittance going to 0 as energy decreases. This is apparent in Fig. 4 where the total transfer probabilities for the TF potential as a function of the energy are presented for both, the multichannel calculation and the simple estimation using Eq. (13). The two calculations agree within 10%.

This suggests a simple procedure for estimating the energy dependence in 3 dimensions (3D). In 3D, the ion-dipole interaction perturbes the $`p\mu _{n=1}`$ initial state at second order through a $`1/R^4`$ term Gershtein (1963). So even for the pure coulombic case there will be a threshold behavior at low energies. We have thus performed a one channel calculation of the transmittance $`T(E)`$ with the 3D functional form of the potential. This $`T(E)`$ was used in Eq. (13) together with $`P_{max}`$ provided by the colinear Coulombic potential to yield an approximate 3-dimensional transfer probability. Since the colinear configuration is the most favorable configuration, this procedure provides an upper limit of the total transfer probability. The results are presented in figure 5 together with those of the colinear calculations. Clearly, the threshold is displaced towards higher energies ($`10^1`$ eV). Thus we expect that for thermal energy collisions ($`\mathrm{4\hspace{0.17em}10}^2`$ eV) the muon transfer probability will be of the order of 0.6 or less and will increase smoothly with energy in the range $`10^210^1`$ eV.

## IV Discussion

In order to make contact with the experiments of Werthmüller et al Werthmüller et al. (1998), we have calculated the muon transfer rate according to:

$$\lambda (E)=Nv\sigma (E)$$
(16)
with $`N`$ being the number density of liquid hydrogen ($`\mathrm{4.25\hspace{0.17em}10}^{22}`$ cm-3), $`v=\sqrt{2E/m_{\text{O},p\mu }}`$ the relative velocity and

$$\sigma (E)=\frac{\pi }{k^2}P(E)$$
(17)
the muon transfer cross section with $`k=\sqrt{2m_{\text{O},p\mu }E}/\mathrm{}`$. In writing Eq. (17), we have assumed only one partial wave with zero angular momentum. This is valid for collision energies up to 0.1–0.2 eV, corresponding to the height of the barrier in the entrance channel for the $`p`$ partial wave.

In Fig. 6 we present our 3D estimates together with the results of the calculations of Sultanov and Adhikari Sultanov and Adhikari (2000) based on a two-states approximation of the Fadeev-Hahn equations. The agreement is good, our results being slightly higher than theirs. This was expected since our 3D estimation is based on the colinear calculations which provides an upper limit to the transfer probabilities. In addition, they have consider only the O$`\mu _{n=5}`$ final channel From Eqs. (16) and (17) it follows that for a constant transfer probability $`P(E)=1`$ the rate should decrease as $`1/\sqrt{E}`$. This is also plotted in Fig. 6. We note that our calculated rates decrease with energy slower than those for a constant transfer probability. The reason is that our calculated probabilities actually increase (although slower than $`1/\sqrt{E}`$) in the energy range considerde here as can be seen in Fig. 5.

In Fig. 6 we have also represented the rates of the two-component model used in Ref. Werthmüller et al. (1998) to fit the experimentally observed X-ray emission from excited muonic oxygen. While for thermal energies there is an order of magnitude agreement between our estimated values and the two-component model, there is a clear disagreement for epithermal energies. The two-component model predicts an increase of the rates going from thermal to epithermal energies in contradiction with our predictions and those of Sultanov and Adhikari Sultanov and Adhikari (2000). Moreover, for energies above 0.1 eV the rate predicted by the two-component model is three times larger than the maximum values of the rate obtained by assuming a muon transfer probability per collision of unity (see Fig. 6).

## V Conclusions

We have presented colinear calculations of muon transfer rates between muonic hydrogen and oxygen for relative translational energies between $`10^6`$ and $`10^3`$ eV. For the lower energies (below $`10^1`$ eV) the de Broglie wavelength is much larger than the characteristic distance of the potential interaction and the problem corresponds to an ultra-cold collision. A simple procedure to estimate the energy dependence of the muon transfer rate in three dimensions was proposed.

Our results show that the muon transfer rate decrease going from thermal to epithermal energies, in agreement with previous theoretical calculations Sultanov and Adhikari (2000) but in contradiction with the assumptions of a two-component model used in Ref. Werthmüller et al. (1998) to interpret the experimental data on this system.

In the experiment the X-ray emission shows a bi-exponential behavior. In the two-component model of Werthmüller et al Werthmüller et al. (1998), it was assumed that the prompt emission is due to muon transfer from epithermal muonic hydrogen while the delayed emission, similar to what is observed when oxygen is replaced by other atoms, is associated to thermalized muonic hydrogen. Assuming that the relative amounts of epithermal and thermal species are equal, Werthmüller et al concluded that the rate of muon transfer from epithermal muonic hydrogen to oxygen should be larger than the one from thermalized ones by a factor of almost 4. They suggest that a resonance could be responsible for this increase. We have shown that for s-waves it is not possible to obtain such an increase in the transfer rate, even if assuming a transfer probability per collision of one. Thus, if resonance effects exist they have to be due to the contribution of higher order partial waves. In an elastic cross section calculation, Kravtsov et al Kravtsov et al. (1996) have shown that for muonic hydrogen colliding with an oxygen nuclei, the contributions of $`p`$ and $`d`$ waves are small as compared to the $`s`$ wave in the range 0–0.2 eV and that only at $``$1.5 eV there is narrow resonance involving the $`d`$ wave. This will thus imply that the epithermal muonic hydrogen has much higher translational energy than what has been assumed in Ref. Werthmüller et al. (1998). Anyway, this conclusion is based on a calculation of the elastic scattering. It could be of interest to check whether it is valid also for the muon transfer process as well.

In this work the calculations were performed in a restricted colinear configuration. This approximation is known to give poor quantitative results in many cases. However for the system considered here, the colinear configuration is the most favorable for the transfer process. Also, since the colinear configuration provides the minimum energy path for the reactance channel, even for an initial non colinear configuration there will be efficient orientational effects in particular for low translational energies. In addition we have introduce a correction to the low energy dependence of the rate in order to take into account the correct asymptotic behavior of the potential in 3-dimensions. Therefore, we believe that although quantitatively the results may change if a full 3-dimensional calculation is performed, the conclusions, in particular concerning the energy dependence of the muon transfer rates, will remain valid. A full 3-dimensional treatment is not impossible however, although quite computationnally demanding. It will definitetly be very important to further improve our detailed understanding of this reaction. It would also be most interesting to study the energy dependence of the muon-transfer rates for other elements, such as C or S, for which the rates seem to behave differently as a function of the energy. This work is in progress.

## VI Acknowledgments

We thank J.M. Launay for helpful discussions, the IDRIS for an allocation of CPU time on a NEC SX5 vector processor, and NATO for a collaborative linkage grant PST.CLG.978454 between Bulgaria and France.