text:
# Deep HST-WFPC2 photometry of NGC 288. II. The Main Sequence Luminosity Function Based on observations made with the NASA/ESA Hubble Space Telescope at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. These observations are associated with proposal # GO-6804.

## 1 Introduction

After a decade of HST observations, the Luminosity Function (LF) of unevolved Main Sequence (MS) stars has been derived for several Galactic globular clusters (GGCs). In many cases the LF extends down to the faintest stellar objects ($`M0.1M_{}`$; see, e.g., King et al., 1998; Bedin et al., 2001). For unevolved stars, the luminosity directly tracks the mass, the most fundamental parameter of a star. The determination of the Initial Mass Function (IMF) from the observed LF of Main Sequence stars, as well as the study of the of the variation of the LF within a given cluster and among different clusters is now a well established field of research which has provided important insights on the dynamical evolution of GCCs (see, e.g., Ferraro et al., 1997; King et al., 1998; Marconi et al., 1998; Piotto & Zoccali, 2000; Paresce & De Marchi, 2000, and references therein). In particular, the signature of mass segregation, due to the energy equipartition established by two body encounters among clusters stars, has been detected in many clusters by comparing the LFs obtained at different distances from the cluster center. Also, the effect of the pruning of the outer extremities of the clusters by the Galactic tidal field or by the violent interaction with the Galactic bulge and disc have been shown with the analysis of the LF of MS stars (see, e.g., Andreuzzi et al., 2001, and references therein).

Here we present the MS LF of the globular cluster NGC 288 down to $`V24.5`$. The LF was based on HST-WFPC2 observations that are described in detail in a companion paper (Bellazzini et al., 2001b, hereafter Pap I), that was devoted to the study of the binary population of this cluster. In Pap I we used a method very similar to the one adopted by Rubenstein & Bailyn (1997), based on extensive sets of artificial star experiments, to estimate the global binary fraction ($`f_b`$) of the cluster ($`7`$ % $`f_b30`$ % ) and to show that binary systems are spatially segregated within the cluster. These same artificial star experiments can be used to obtain LFs in different regions of the cluster, corrected for all observational effects. In §2 we will briefly recall the characteristics of the adopted sample and of the artificial star experiments, as well as the results of Pap I that are relevant in the present context. In §3 we introduce a new method to correct the observed LF. The effect of the presence of binary systems on the final LF is also discussed and quantified with the support of the analysis performed in Pap I. In §4 we present the final LF of NGC 288 and compare it with the LFs of other clusters and, finally, in §5 we summarize our results.

### 1.1 The globular cluster NGC 288

NGC 288 is a loose ( $`\mathrm{log}\rho _0=1.80L_Vpc^3`$; Djorgovski, 1993) cluster. It is a classical old GGC (Rosenberg et al., 1999) with a blue Horizontal Branch (HB) morphology (see Bellazzini et al., 2001a; Catelan et al., 2001, and references therein). Its stars are not exceedingly metal deficient ($`[\mathrm{Fe}/\mathrm{H}]=1.39`$) and have the abundance pattern typical of halo stars (Shetrone & Keane, 2000). The cluster lies on a very inclined and eccentric orbit that carries it into the inner, denser region of the Galaxy, where the effects of Galactic tides and disc/bulge shocks are more disruptive . According to Dinescu, Girard & van Altena (1999) the perigalactic distance is $`R_p1.8`$ kpc, the eccentricity is $`e0.7`$, the orbital period is $`P220\times 10^6`$ years, and the “destruction rate” is one of the highest in the GGC system. The observed distance from the Galactic Center (Harris, 1996, $`R_{GC}11`$ kpc;) indicates that the cluster is presently near the apogalactic point of its orbit (see also Pasquali, Brigas & De Marchi, 2000).

## 2 The dataset and the artificial star experiment

The observations, the data reduction process and the artificial star experiments are described in detail in Pap I. Here we summarize the points that are relevant in the present context.

### 2.1 The sample

Two HST-WFPC2 fields were observed in the F555W and F814W passbands: one with the PC camera nearly centered on the center of the cluster (Int field), and one with the WF4 camera partially overlapping the WF3-Int field (Ext field; see Fig. 1 of Pap I). The Int field is almost completely enclosed in the region within 1 half light radius ($`1r_h1.6r_c`$; where $`r_c=85^\mathrm{"}`$ is the core radius) Djorgovski, 1993) while the Ext field samples a region contained in the annulus $`1r_h<r<2r_h`$, where $`r`$ is the distance from the cluster center. The data were transformed to the standard Johnson-Kron-Cousins photometric system by calibration to the ground-based data of Rosenberg et al. (2000). The Int sample contains 5766 sources from $`V13`$ to $`V25`$ and the Ext one contains 2013 sources between $`V16`$ and $`V25.5`$. The large majority of the observed stars are fainter than the Turn Off (TO) point ($`V19`$), and hence are unevolved MS stars (and/or binary systems containing two unevolved MS stars; see Pap I). The F336W and F225W observations of the Int field that has been shortly discussed in Pap I are not used in the present paper.

### 2.2 The artificial star experiments

Simulations containing more than 80,000 artificial stars per chip have been made on the $`V,I`$ frames. In each run of the artificial stars experiments $`100`$ artificial stars per chip were added at the same position in the median image (that we used as a master-frame to obtain a list of positions of bona fide stars; see Pap I) and in each of the single frames (that were used to obtain accurate photometry; see Pap I). The set of experiments for a given chip (sub-field) was completed with $`800`$ independent runs. The final total number of artificial star experiments is $`>1,500,000`$. For each artificial star a random input $`V`$ magnitude is drawn from a distribution similar to the observed LF extrapolated to magnitudes slightly below the observed lower limit. This fainter limit was required so that the simulations would contain a large number of stars fainter than the observed limiting magnitude to correctly sample the magnitude bins where the incompleteness is expected to be highest (see Pap I).

We ultimately derive the true LF in a way that does not depend on the assumed artificial star LF (see §3.2). The input $`I`$ is obtained from the input $`V`$ using the cluster ridge line. Each simulated star thus had the appropriate color and was simultaneously added to the corresponding $`V`$ and $`I`$ frames at the same position. The simulated stars were randomly distributed on the frames with an additional constraint (see Pap I and also Tosi et al., 2001, for a detailed description of the method) preventing any interference between simulated stars. In practice, for each run of the artificial stars experiment, only one artificial star is added in each $`80\times 80`$ pixel portion of a chip. The whole process of data reduction has been repeated on all the frames with artificial stars added and the output magnitudes have been recorded (see Pap I, for details).

Fig. 1 shows the input and output CMDs for the stars simulated in the PC-Int, PC-Ext, WF2-int and WF2-Ext frames. The WF2 samples have been chosen as representative of the WF cameras for the Int and Ext fields. Within a given field (Int and Ext) the observations with the different WF cameras have very similar properties (see Pap I). Note the different response to the same input color-magnitude distribution in the four cases represented. This arises because the larger pixels of the WF cameras allow them to reach fainter limiting magnitudes than the PC camera and because the Int and Ext fields have different average degrees of crowding.

All of the effects of the whole process of observation and data reduction on a large sample of stars representing a pure Simple Stellar Population (SSP) are shown in Fig. 1. These effects can be imagined to be the results of the application of a multi-dimensional function (that we will call Observation/Measure Function; hereafter OMF) on the true stars of the clusters, here represented by the input artificial stars lying on the ridge-line. In Pap I the empirical knowledge of the effects of the OMF has been used to extract the signal provided by the presence of a significant population of binary systems in the color distribution around the ridge-line from the observational noise due to random scatter, blending, etc. Here we will use it to correct the observed LF from all the spurious deformations due to the same factors.

### 2.3 Completeness and division in subsamples

The artificial star experiments have been carried out independently for each camera of each observed field. The completeness functions for each camera have been presented in Fig. 5 of Pap I. The crowding conditions of the dataset are never particularly critical, and the completeness is rather similar everywhere. The completeness factor ($`C_f`$) is larger than 80 % for $`V23.5`$ in all the samples except the PC-Int, the innermost one, where $`C_f80`$ % for $`V22`$. It has been shown that the spatial variation of completeness within a single camera are negligible. For the present purpose, we take advantage of the great similarity of crowding conditions in the WF cameras of the Int and Ext samples to subdivide the total sample in the following subsamples: (1) PC-Int, (2) WF-Int, containing the WF2-Int, WF3-Int and WF4-Int samples, (3) PC-Ext and, (4) WF-Ext, containing the WF2-Ext, WF3-Ext and WF4-Ext samples. The four subsamples are internally very homogeneous and cover different radial regions of the cluster. In this sense, the difference between the PC-Ext and the WF-Ext samples is small but we prefered to keep them separate because of the different limiting magnitude.

## 3 The correction of the observed LF

The effects of the OMF on the true properties of stars that may affect the LF derived from an observed sample are:

1. Star loss. Stars may be lost in the observation/measure process. This occurs for two different reasons:

1. Near the limiting magnitude of the observations a star can be succesfully detected if it falls on a positive fluctuation of the background that brings it above the detection threshold and can be lost if it falls on a negative fluctuation. This effect is purely stochastic; it does not depend on the crowding. It just depends on the amplitude of the random fluctuations of the background. There is a level of magnitude beyond which no possible positive fluctuation can bring a star above the detection threshold. All the stars fainter than this limit are lost (unobserved) and just contribute to the diffuse luminosity of the unresolved background.
    2. Stars may be lost if they fall on a brighter source (star or galaxy) and/or if their images are badly corrupted by defects of the chip or the hit of a cosmic ray (this two latter effects are expected to have a negligible impact in the present case because of the adopted observation + data reduction strategy; see Pap I). This effect depends on the local degree of crowding as well as on the magnitude of the star.
2. Blending. A star may be blended with a fainter source and its measured magnitude will be brighter than the true one because of the flux contributed by the fainter blended source. Considered from the perspective of artificial stars experiments the maximum possible effect of a blending is to change the observed magnitude of a star by $`0.75`$ mag. In fact, if we recover an artificial star by its (known) position in a frame with an output magnitude more than $`0.75`$ mag brighter than its input magnitude, this means that it is fallen on a real star brighter than the simulated star. Thus, we are not recovering/measuring the simulated star but the real one instead, i.e., thus the artificial star has been lost.
3. Observational scatter. The magnitude of the true star is altered by all the sources of stocastic noise associated with the process of observation/data reduction (e.g., photon noise, imperfect PSF fit, etc.). The maximum possible effect of this factor is magnitude dependent and may be roughly evaluated on the basis of the typical error on a given measure. For instance, in the present case we included in the final samples only the stars with errors111The error associated to each magnitude entry is the error on the mean over repeated measures, see Pap I. lower than 0.2 mag in both $`I`$ and $`V`$. Hence a $`3\sigma `$ fluctuation can change the magnitude of an observed star by $`0.6`$ mag, at most, near the limiting magnitude level.

The derived LFs are in practice histograms of the number of stars as a function of magnitude. The factors above change the true LF histogram by depopulating the bins as a function of (increasing) magnitude (mainly factor 1) and by moving stars from one bin to another (migration) due to the magnitude changes induced by factor 2 and/or 3. However factors 2 and 3 are strictly interwoven and their effect cannot be separated. For example, a given star may be moved, say, 0.4 mag brighter because it is blended with a fainter unresolved source, but it might simultaneously be affected by a negative fluctuation of the random observational scatter, so its final magnitude may differ from the true one by just 0.2 mag. A blending between two unresolved sources can bring them above the detection threshold, injecting in the sample two stars that should have been lost if not blended222Note that this phenomenon may significantly alter the star counts in the bins near the limiting magnitude, at least in cases in which the crowding conditions are critical, see Tosi et al. (2001), for discussion. The possibilities are too numerous to list here. The real contribution of the three effects cannot be known for a the specific case. On the other hand, it is clear that we can obtain a statistical knowledge of the actual LF.

The traditional way to recover the true LF from an observed one consists in multiplying the numbers in each bin of the observed LF by the inverse of the Completeness factor ($`C_f=N_R/N_S`$; where $`N_R`$ is the number of recovered stars and $`N_S`$ is the number of the simulated ones) at the center of the bins. It is clear that such procedure corrects at first order only for the effects of star loss (factor 1) while the effects of the other factors are ignored. A popular misconception is that the observational scatter cannot produce bin migration if the width of the bin is larger than the typical error on the estimate of the magnitude. This is obviously not true since a star whose true magnitude is near the edge of a bin can migrate because of an arbitrarily small fluctuation. The effects of migration may be amplified by the completeness correction, especially in the faintest bins of the LF.

Furthermore the Completeness factor may be obtained from the artificial star experiments only as a function of input magnitude (say, $`V_{\mathrm{in}}`$) while the observed LF is (obviously) a function of the observed magnitude, whose corresponding quantity in the artificial star experiment is $`V_{\mathrm{out}}`$333The input magnitude ($`V_{\mathrm{in}}`$) of artificial stars corresponds to the true magnitude of real stars ($`V_{\mathrm{true}}`$), while the output magnitude ($`V_{\mathrm{out}}`$) corresponds to the observed one ($`V_{obs}`$) under the extant conditions. Thus when an observed LF is corrected by multiplying it by $`1/C_f(V_{\mathrm{in}})`$, it is implicitly assumed that $`V_{\mathrm{in}}`$ and $`V_{obs}`$ are the same quantity, that, in general, is not true.

There have been few attempts to correct the observed LFs for all the effects above (Drukier et al., 1988; Stetson & Harris, 1988). The proposed methods require a complex mathematical approach, and may rely on the a priori assumption of an analytical form of the true LF (Stetson & Harris, 1988) or on approximate treatments of the propagation of errors under matrix inversions (see Drukier et al., 1988; Mighell, 1990, for references and discussion).

Here we propose a new method that is fully non-parametric and extremely simple in nature, taking also advantage of the huge enhancement in computational power of workstaions and PCs that took place since the times of the above quoted attempts. As we will show below (a) the proposed method takes full account of the effects of bin migration and (b) the completeness correction is consistently applied to a LF that is function of $`V_{\mathrm{in}}`$ as the completeness factor. The method is particularly appropriate for the application to LFs of SSPs, but it can be probably generalized to other cases. The prerequisite is the availability of huge sets of realistic artificial star experiments well sampling the whole spatial, luminosity, and color range of the observations.

### 3.1 The Equivalent Sample

Consider a large set of artificial star experiments performed in a region of uniform crowding conditions, e.g., the WF2-Ext set presented in Fig. 1. Suppose that the input magnitudes of the artificial stars be distributed as the “true” stars of the parent population (in this case the stars of NGC 288; the reasons of this hypothesis will become clear at the end ofthis section). From this set it is possible to extract a subsample of recovered stars with exactely the same dimension and the same distribution as a function of $`V_{\mathrm{out}}`$ as the observed sample. Call this subsample an Equivalent Sample (ES)444The operational procedure to extract an ES is the following. Consider the histogram representing the observed LF: each bin i contains $`N_i`$ stars. From the artificial star set we extract at random $`N_i`$ stars whose output magnitudes lie in the i-th bin. Repeating this step for each bin of the observed LF we obtain a sample of recovered artificial star that have an output LF equal to the observed one. Note that the same procedure can be applied on a star-by-star basis, by randomly extracting an artificial sister (see §3.4) for each observed star.. It is important to note that the ES contains only stars that have been succesfully recovered in the artificial star experiments process, hence that have suffered only the effects of bin migration (i.e, blending and/or observational scatter).

By definition the ES has the same output (e.g., observed) properties of the observed sample, it is one of the many possible relizations that are mapped by the OMF from the input/true domain to the output/observed one. Thus the input (e.g., true) properties of the ES are one of the possible true set of properties that may have been mapped into the observed sample by the OMF. It is immediately apparent that the LF of the input magnitudes of the ES is equivalent to the observed LF once corrected for all the effects of bin migration555In other words the LF of the input magnitudes of the ES is the observed LF seen before any blending and/or observational error take place.. Thus the identification of an ES allows (a) the removal of all the effects of bin migration and (b) the expression of the observed LF as a function of input magnitude, consistently with the completeness factor (see §3.). If we apply the completeness correction to the input LF of the ES we find one possible realization of the final LF, corrected for all the observational effects. Since different ES may be mapped into the same observed sample by the OMF, repeating the process for a number of (randomly extracted) different ES and taking the average will allow to obtain a final LF more robust to random fluctuations as well as a direct estimate of the uncertainty associated with the whole process.

However, there is another important factor we have to take into account. While the completeness correction is completely independent from the distribution in magnitude of the artificial star set, the derived correction for bin migration is not. Consider a given bin of the observed LF. Such a bin will be primarily filled by stars with $`V_{\mathrm{out}}V_{\mathrm{in}}`$, but also by a certain fraction of stars which have migrated from the nearby bins beacause of blending and photometric errors. Increasing the number of artificial stars in the nearby bins, with respect to the considered bin, will enhance the probability of extracting stars that enter in the bin because of migration, hence increasing the final “migration correction”. The right correction would be obtained only if the artificial star set is distributed as the true LF of the population which is unknown, being the very target of our analysis. We circumvent this problem by iteratively adjusting (a posteriori) the distribution of the artificial star set until the best approximation is obtained (see below for details). The approach has proven to be successful since we obtain the same final LF independently of the initial distribution of the artificial star set. In the following section the adopted operational procedure is described in detail.

### 3.2 Operational procedure

For each of the four observed subsamples (PC-Int, WF-Int, PC-Ext, WF-Ext) we adopted the following procedure:

1. An initial distribution of input magnitudes is assumed for the artificial stars and the largest possible subset accordingly distributed is extracted from the whole set of artificial stars.
2. Twenty different ES are extracted at random from the above subset and their $`LF(V_{\mathrm{in}})`$ are corrected for completeness.
3. The (bin per bin) average of the twenty obtained LFs is adopted as the “current corrected LF” and the standard deviation is adopted as the corresponding uncertainty.
4. The “current corrected LF” is assumed as the distribution of the artificial stars and a new subset of artificial stars accordingly distributed is extracted from the whole set of artificial stars as in step 1.
5. The steps 2,3 and 4 are repeated twenty times. The twentieth “current corrected LF” is adopted as the final LF with the corresponding uncertainties.

Step 4 ensures that at each new iteration the adopted set of artificial stars will become more and more similar to the “true” one, giving increasingly appropriate rates of bin migration.

In all the cases considered here the process converged in $`6`$ iterations, e.g., the “current corrected LFs” become very stable after the sixth iteration. Furthermore, the result is independent of the initial distribution of input magnitudes that is assumed for the artificial star subset. We have tried many different distributions (some of which are shown in Fig. 2) and in all cases the process converged to the same final LF in few iterations.

Some examples of the procedure are shown in Fig. 2, applyied to the WF-Ext sample. In Fig. 2a the artificial star subset was initially distributed as the observed LF down to $`V25.5`$ and uniformly from this limit down to $`V=27`$. A scaled version of the adopted distribution is plotted as a thick dotted line. The thin dotted histogram shows the “current corrected LF” after the first iteration of the process, the thin continuous histogram shows the “current corrected LF” after the sixth iteration and the open dots are the final LF, corresponding to the twentieth iteration of the process, with the currrent standard deviation of the twenty realizations of the corrected LFs as errorbars. The encircled dot is the first bin of the LF for which $`C_f<0.5`$, i.e. the completeness correction is larger than a factor 2. The fast convergence and the stability of the final solution can be readily appreciated.

Fig. 2b and 2c have the same symbols of Fig. 2a but shows the convergence of the process starting from (very) different initial distributions of the artificial star subset. The open dots are the same final LF shown in Fig. 2a. Thus Fig. 2a and 2b shows that the process converges to the same final LF independently of the assumed initial distributions of the artificial stars subset.

Finally, Fig. 2d shows in detail the results of the twentieth (i.e. final) iteration of the process shown in Fig. 2a. The twenty histograms are the corrected LFs obtained by the twenty different ES extracted from the finally adopted subset of artificial stars. The full dots and errorbars are their means and standard deviations.

#### 3.2.1 The final LFs

Since the obtained LFs are independent from the initial distribution of the adopted subset of artificial stars, we derived our final (corrected) LFs for all the considered samples assuming the same initial distribution of the artificial stars shown in Fig. 2a. The final LFs are given in Table 1 (PC-Int and WF-Int samples) and in Table 2 (PC-Ext and WF-Ext samples). The table gives the $`V_{\mathrm{in}}`$ (e.g., $`V_{\mathrm{true}}`$) of the center of the bins and, for each sample, the observed LF, the completeness factor $`C_f`$, the final LF and the total error on the final LF (see below).

The final LF is given down to the first bin that has $`C_f<0.50`$. In the following plots this point is circled to serve as a reminder that it is has been derived applying a large correction. We retain these points as educated guesses of the behaviour of the LFs beyond the range in which moderate and safe corrections can be made. We consider as fully reliable only the points of the LF for which $`C_f0.50`$. We note that the last fully reliable points, according to the above criterium, have $`C_f=0.607,0.769,0.806,0.842`$ for the PC-Int, WF-Int, PC-Ext, and WF-Ext samples, respectively. Thus the completeness is quite high for the whole considered range.

The final total error reported in Tab. 1 and 2 has been obtained by summing in quadrature the error on the Completeness factor (estimated according to Eq. 18 by Mighell, 1990), to the standard deviation of the twenty different realizations of the LF that are obtained in the last step of the procedure (see §3.2 and Fig. 2d). In all the following plots the errorbars associated with the LFs of NGC 288 represent this final total error.

The final LFs (open dots) are shown with the corresponding observed LF (full dots) in Fig. 3.

As a further check we compared our final LFs with those obtained by dividing the observed $`LF(V_{obs})`$ by $`C_f(V_{\mathrm{in}})`$. We note that the differences are within $`1\sigma `$ in most bins and $`<2\sigma `$ in any bin. As expected, in this particular application the effects of bin migration are not large. This would not be the case if much more crowded fields were considered (see §3. and Tosi et al., 2001).

### 3.3 The effect of binary systems on the LF

There is an intrinsic kind of blending that cannot be directly tackled with artificial star experiments. Physical binary systems are strictly equivalent to chance superposition blendings. The LF for single stars cannot be obtained if an independent estimate of the incidence of binary systems is not available. Further, since the actual binary fraction ($`f_b`$) depends on the (unknown) distribution of mass ratios ($`F(q)`$; see Pap I), any correction would be model dependent. Indeed, a direct estimate of the impact of binary systems on the LF of a globular cluster has only been possible for NGC 6752 (Rubenstein & Bailyn, 1999).

Now we can add a second cluster using the results of Pap I. We measured $`f_b`$ in the Int and Ext regions of NGC 288. In the Int region we found that the observations are compatible with 8%$`f_b`$ 38%, and the most probable range is 10% $`f_b`$ 20%. In the Ext region $`f_b`$ 10 % and the most probable value of $`f_b`$ is zero, independently of the assumed $`F(q)`$. Thus, a sizeable fraction of binary systems is present in NGC 288, at least in the Int sample, and may affect the derived LF of single stars.

To quantify the effect we computed the mean bin-by-bin ratio between the LFs obtained by 5 ES with $`f_b=20`$ % and by 5 ES with $`f_b=0`$ %. This number gives (approximately) the fractional variation of the star content of the bins due to binary systems. Furthermore, dividing the final LF by the quoted ratio, the LF is corrected for the binary content and mapped to the “pure single stars” case. In Fig. 4 the final LF of the WF-Int sample not corrected for the binary content (filled dots with error bars) is compared with its binary-corrected versions in the assumption $`f_b=20`$ % and for three different $`F(q)`$ (see caption and Pap I for details). The comparison is performed in the magnitude range in which the binary estimate of Pap I has been performed ($`20V23.5`$).

There are several noteworthy features in Fig. 4: 1) The maximum correction is $``$ 3%. 2) All the corrected LFs lie within the uncertainties of the un-corrected one. 3) The PLMR $`F(q)`$ case is most similar to the uncorrected case. This was expected, since in this case the large majority of binary system has a very low mass secondary. This means that the migration of the primary is small, and that most of the secondary stars are so faint that they fall in an unobserved region of the LF. 4) The correction for the PHMR case depletes the intermediate bins ($`21.5V23`$; because it corrects for the migration of primary) and enhances the last bin because of the recovery of some secondary stars as single stars. 5) No bin is exclusively depleted or enhanced by the redistribution of the binary components as single stars.

From the above results and discussion we conclude that:

1. The correction for binary systems can be safely neglected for the Ext samples where the binary fraction is much lower than $`f_b=20`$ % and possibly null.
2. The correction can be neglected also for the Int samples, at least for $`V<23.5`$, since it is smaller than the uncertainties of the final LF. However, it cannot be excluded that the faintest bins of the single stars LF need a significant enhancement to take into account the effect of binaries, particularly if the actual $`F(q)`$ favors high mass ratios.

### 3.4 Further developements and other applications

The case under study is particularly simple. NGC 288 is a very loose and relatively nearby cluster, so our WFPC2 observations well resolve the large majority of the stars down to the limiting magnitude of each subsample with a very high degree of completeness. Further, the density gradient is a very fortuitous match to WFPC2 allowing a few large subsamples to nicely probe radial variations.

It is important to note that the Equivalent Sample concept does not require such simplifications and may be successfully applied under any condition. If a strictly rigorous approach is adopted the ES approach may have far reaching applications in the study of stellar populations. Here we provide two possible examples to show the potential of the method.

* Extremely crowded field with strong density gradient. The method can tackle with samples with arbitrarily large density gradient, for instance by subdividing the observed sample in subsamples whith nearly uniform stellar density (as done in the present application). If the dimension of the observed sample is not sufficiently large to allow this approach, the ES method can be implemented also on a star-by-star basis. It will suffice to assemble the ESs by associating an artificial star analogue for each observed star, randomly extracted from a set of artificial stars having (approximately) the same $`V_{\mathrm{out}}`$ and position in the frame of the considered observed star (artificial sisters). The only basic requirement to successfully follow this approach is to have a very large set of artificial stars, sampling with large multiplicity the whole range of magnitudes and positions covered by the observed sample.
* Mapping synthetic CMDs into the observational plane. The standard tool used to derive the Star Formation History (SFH) from the CMD of a resolved galaxy consist in the reproduction of observed CMD by a synthethic CMD in which the stars are extracted from theoretical evolutionary tracks (see Lejeune et al., 2001, for references and discussion of the various adopted techniques). To transform the synthetic CMD to the observational plane the effects of incompleteness, blending and observational scatter must be added to the synthetic stars. This task is not easy and it is usually done with simplified approaches. For instance, the observational scatter is usually introduced by randomly extracting ‘errors’ from gaussian distributions with standard deviations equal to the typical $`1\sigma `$ photometric error. However, the distribution of the observational scatter is not gaussian, in general, and it is not easy to parametrize. Furthermore, the effects of blending are not separable from the effects of observational scatter, and gaussian distributions do not provide a good reproduction of the real distribution, introducing undesired noise in the comparisons.

With the ES approach, one can easily map the input magnitudes (e.g., “true” magnitudes) extracted by the evolutionary track in to the observational plane with the correct OMF, by “asking” to the artificial stars experiments what is the probabilty that such star would be successfully recovered or not. If a star passes the “incompleteness barrier”, it can be correctely mapped into the observational plane, for instance by assigning to it the output properties of an artificial analog extracted at random from a set of artificial sisters.

In this way, the whole effects of the OMF would be correctly reproduced without dangerous approximations and/or parametrizations.

## 4 Comparison of LFs

Once the observed LFs have been corrected for all the observational effects and the impact of binary systems have been quantified, we can use the final LF to study the dynamical status of NGC 288 either comparing the LFs in the Int and Ext regions of the cluster, either comparing the LF of NGC 288 with those of other clusters. In the following we will normalize all the samples to be compared to the number of stars in the range $`18.5V20.5`$, following the approach of Piotto & Zoccali (1999, hereafter PZ99).

NGC 288 is not sufficiently nearby to allow a safe derivation of the LF down to the hydrogen burning limit even with HST. The comparison with other “state of the art” LFs of more nearby clusters (see Fig. 6, below) shows that such limit would occur a couple of magnitudes below our faintest valid point (the faintest bin of the WF-Ext LF). The Mass-Luminosity relation adopted in Pap I indicates that the lowest mass efficiently sampled by our deepest LF is $`M0.35M_{}`$. Usually MF are compared according to their slope in the plane $`\mathrm{log}N\mathrm{log}M`$ for $`M0.5M_{}`$ (see PZ99 and references therein). Our best LF has just three points below this limit, thus the slope of the MF would be poorly constrained. Because of this limitation we avoid any conversion of our LFs to Mass Functions (MF), since that would necessarily be somewhat model dependent and uncertain (see, e.g., Bedin et al., 2001). Instead we will make direct comparisons between LFs in the allowed range.

### 4.1 Mass segregation within single stars in NGC 288

The LFs of the PC fields are identical (to within the errors) to those of the respective WF samples, in the range where the comparison is possible. Since the WF samples are much larger and reach fainter magnitudes we will limit our analysis to these homogeneous datasets in the following.

In Fig. 5 the LFs from the WF-Int and the WF-Ext samples are compared. It is evident at a first glance that the LF of the WF-Int sample is flatter than that of WF-Ext. Given the adopted normalization, this is a clear indication of a depletion of low mass stars in the region within $`1r_h`$ (Int) with respect to the Ext field. If the brightest bin is included in the normalization the effect is significantly enhanced. We note also that any correction for binaries would be much smaller than the observed difference (see Fig. 4). On the other hand, if fainter stars ($`20.5<V24`$) are used for normalization, one would infer a substantial excess of the more massive stars in the Int sample.

Independently of the assumed normalization it is worth establishing if the difference between the two LFs is statistically significant. To check this hypothesis we performed a $`\chi ^2`$ test, adopting the appropriate definition of the statistic according to Eq. 14.3.3 by Press et al. (1992). The comparison has been performed in the safely comparable range $`18.5V24`$. The resulting reduced $`\chi ^2`$ is 2.31 with 12 degrees of freedom. The difference between the two LFs is found to be significant at the $`99.1`$ % confidence level, i.e. highly significant.

We have already demonstrated that mass segregation exists in NGC 288, since both binary systems and their by-product Blue Stragglers Stars are also more centrally concentrated than the the single stars. Although the evidence is not as strong (because of the possible ambiguity in the interpretation associated with the choice of the normalization range), we feel that we now find some evidence for further mass segregation in the difference between the LFs of unevolved stars in the Int and Ext regions of the cluster.

### 4.2 Low mass star depletion in NGC 288

In Fig. 6 the WF-Int and the WF-Ext LFs of NGC 288 are compared with the LFs of M10, M22 and M55 by PZ99. The transformation from apparent to absolute magnitude has been obtained by assuming $`(mM)_0=14.67`$ and $`E(BV)=0.03`$, according to Ferraro et al. (1999). The faintest point of the WF-Ext LF reaches $`M_V10`$ while the LF of the (more nearby) PZ99 clusters reaches $`M_V12`$–13. The faintest reported point appear to coincide with the turn-over point found in the LFs of the PZ99 clusters, so this possible feature in the LF of NGC 288 is not detectable in our data (but it have been apparently detected by Pasquali, Brigas & De Marchi, 2000, hereafter PBD00).

On the other hand, the comparison in the range $`6<M_V10`$ shows that the LF of NGC 288 is significantly flatter than those of M10, M22 and M55. Note that the LFs of these clusters were obtained near the half-mass radius (usually approximated with the half light radius $`r_h`$, see Djorgovski, 1993). This is the region where the mass distribution of the cluster population is expected to be fairly representative of the (present day) global population (see Meylan & Heggie, 1997; Vesperini & Heggie, 1997, and references therein). Since the Int and Ext sample are separated at 1 $`r_h`$, the LF at the half-mass radius must be intermediate between the WF-Int and the WF-Ext, and the derived LFs should be representative of the global LF of NGC 288.

The comparison presented in Fig. 6 strongly suggests that the global population of NGC 288 is significantly depleted of low mass stars with respect to M10, M22 and M55 (see PZ99 for an indirect comparison with other clusters).

NGC 288 has a mass ranging from $`\frac{1}{2}`$ to $`\frac{1}{6}`$ that of the comparison clusters. Its central density ranges from $`\frac{1}{2}`$ to $`\frac{1}{60}`$ (Pryor & Meylan, 1993) that of the comparison clusters. Thus NGC 288 is intrinsically less resistent to harrassment and heating by external forces. Furthemore, among the clusters considered here, it moves on the orbit having the smallest perigalactic distance, the highest eccentricity and the highest inclination to the Galactic plane. Hence, it is expected to suffer strong bulge and disk shocks, and indeed may be one of the clusters with the highest disruption rates (Dinescu, Girard & van Altena, 1999; Gnedin & Ostriker, 1997). The coupled effects of internal energy equipartition (mass segregation), pushing the lighter stars to outer, shallowly bound orbits, and of the tidal field of the Galaxy, removing the less bound stars venturing in the outer parts of the cluster halo, may well be at the origin of the depletion of lighter stars in NGC 288. We regard this result as a reassuring confirmation of the theoretical predictions that account for the orbital characteristics of actual GGCs (Dinescu, Girard & van Altena, 1999; Gnedin & Ostriker, 1997).

Recently, PBD00 presented a LF in the $`J`$ and $`H`$ near-infrared bands obtained with the HST NIC3-NICMOS camera (parallel observations). The sample is small but reachs fainter magnitudes (lighter masses) than our data. PBD00 conclude that the mass function of NGC 288 is very similar to that of all other studied clusters, and that it shows no sign of the strong harrassment predicted by Dinescu, Girard & van Altena (1999) and Gnedin & Ostriker (1997). We regard the results (and the conclusions) of the present paper as more robust than those by PBD00 since: (a) the LF by PBD00 is based on a sample of 75 stars, a very small sample that may be subject to strong statistical fluctuations, and is not supported by artificial star experiments, while our LF is based on a sample of more than $`6,000`$ stars corrected for all the observational biases with the artificial star experiments; and (b) the field observed by PBD00 is a small stamp ($`51\stackrel{}{\mathrm{.}}2\times 51\stackrel{}{\mathrm{.}}2`$, $`15`$ times smaller than the field observed here) located at $`r2.4r_h`$, far outside the region of our sample which is expected to be representative of the global cluster LF.

## 5 Summary and Conclusions

The LF of the globular cluster NGC 288 (down to $`M_V10`$) has been obtained from a field covering a region comprised between the center of the cluster and $`r2r_h`$.

A new method to correct the observed LFs for all the observational effects have been introduced and applied. The method is based on the Equivalent Sample concept, which may have many interesting applications in the study of stellar populations.

The effect of the presence of binary systems in the final LF has been quantified and it has been found negligible, in the considered magnitude range.

Comparison of the LFs obtained in different regions of the clusters, indicates the presence of mass segregation within NGC 288. Independently of the assumed normalization the LFs of the inner (WF-Int) and outer (WF-Ext) sample turns out to be different at the $`99.1`$ % confidence level.

The comparison with the LFs of other clusters strongly suggests that the global population of unevolved stars in NGC 288 is significantly depleted of low mass stars ($`M0.5`$$`0.6M_{}`$), in general agreement with the predictions of recent theoretical predictions taking into account also the orbital properties of the cluster.

The financial support of the Agenzia Spaziale Italiana (ASI) is kindly acknowledged. The financial support to M. Messineo has been provided by the Osservatorio Astronomico di Bologna. Part of this work has been the subject of the Thesis of Degree of L. Monaco (Dept. of Astronomy, Bologna University). Part of the data analysis has been performed using software developed by P. Montegriffo at the Osservatorio Astronomico di Bologna. RTR is partially supported by NASA LTSA Grant NAG 5-6403 and STScI grant GO-8709. This research has made use of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# Ehrenfest-time-dependent excitation gap in a chaotic Andreev billiard

(February 2002)

## Abstract

A semiclassical theory is developed for the appearance of an excitation gap in a ballistic chaotic cavity connected by a point contact to a superconductor. Diffraction at the point contact is a singular perturbation in the limit $`\mathrm{}0`$, which opens up a gap $`E_{\mathrm{gap}}`$ in the excitation spectrum. The time scale $`\mathrm{}/E_{\mathrm{gap}}\alpha ^1\mathrm{ln}\mathrm{}`$ (with $`\alpha `$ the Lyapunov exponent) is the Ehrenfest time, the characteristic time scale of quantum chaos.

The density of states in a normal metal is suppressed near the Fermi energy when it is brought into contact with a superconductor. The history of this proximity effect goes back to the 1960’s Gen63 . It was understood early on Mcm68  that the energy range of the suppression is the inverse of the typical life time $`\tau _\mathrm{c}`$ of an electron or hole quasiparticle in the normal metal. This life time is finite (even at zero temperature) because an electron is converted into a hole by Andreev reflection at the interface with the superconductor And64 . The energy scale $`E_\mathrm{c}=\mathrm{}/\tau _\mathrm{c}`$, known as the Thouless energy, is the product of the mean level spacing $`\delta `$ in the normal metal and the dimensionless conductance of the contact to the superconductor. (For example, $`E_\mathrm{c}=N\delta `$ for coupling via an $`N`$-channel ballistic point contact.) The appearance of an excitation gap of the order of the Thouless energy is the essence of the traditional proximity effect.

Some years ago it was realized Mel96 ; Lod98 ; Sch99 ; Pil00 ; Ihr01 ; Tar01  that the proximity effect is essentially different if the normal metal becomes so small and clean that scattering by impurities can be neglected. This applies to a quantum dot in a two-dimensional electron gas Wee97 , and because of the resemblance to a billiard (cf. Fig. 1) one speaks of an “Andreev billiard” Kos95 ; Wie02 . Depending on the shape of the billiard, the classical dynamics varies between integrable and chaotic. No excitation gap is induced by the proximity effect in an integrable billiard Mel96 ; Ihr01 . An excitation gap does appear in a chaotic billiard Mel96 ; Sch99 , but its magnitude is only given by the Thouless energy if the chaos sets in sufficiently rapidly Lod98 ; Tar01 .

The characteristic time scale of quantum chaos is the Ehrenfest time $`\tau _\mathrm{E}=\alpha ^1\mathrm{ln}(L/\lambda _\mathrm{F})`$, defined in terms of the Lyapunov exponent $`\alpha `$ (being the rate at which nearby trajectories diverge exponentially in time) and the relative magnitude of the Fermi wave length $`\lambda _\mathrm{F}=2\pi /k_\mathrm{F}`$ and a typical dimension $`L`$ of the billiard Ale96 . Chaotic dynamics requires $`\alpha ^1\tau _\mathrm{c}`$, but $`\tau _\mathrm{E}`$ could be either smaller or larger than $`\tau _\mathrm{c}`$. In the regime $`\tau _\mathrm{E}\tau _\mathrm{c}`$ the excitation gap is set as usual by the Thouless energy. Established techniques (random-matrix theory, non-linear $`\sigma `$-model) provide a complete description of this regime Mel96 ; Vav01 ; Ost01 ; Lam01 . The opposite regime $`\tau _\mathrm{E}\tau _\mathrm{c}`$ has no analog in the conventional proximity effect. Random-matrix theory is helpless and this regime has also shown a frustrating resilience to solution by means of the ballistic $`\sigma `$-model Tar01 . In particular, no mechanism has yet been demonstrated to produce the hard gap at $`\mathrm{}/\tau _\mathrm{E}`$ conjectured by Lodder and Nazarov Lod98 .

Here we report an attack on this problem by an alternative approach, starting from the semiclassical Andreev approximation to the Bogoliubov-De Gennes (BdG) equation And64 . The limit $`\tau _\mathrm{E}\mathrm{}`$ yields the Bohr-Sommerfeld approximation to the density of states Mel96 ; Lod98 ; Sch99 ,

$$\rho _{\mathrm{BS}}(E)=\frac{2}{\delta }\frac{(E_\mathrm{c}/4E)^2\mathrm{cosh}(E_\mathrm{c}/4E)}{\mathrm{sinh}^2(E_\mathrm{c}/4E)},$$
(1)
which is gapless (cf. Fig. 1). We have found that diffraction at the contact with the superconductor is a singular perturbation to $`\rho _{\mathrm{BS}}`$ that opens up a gap at the inverse Ehrenfest time, and provides an intuitively appealing mechanism for the gap phenomenon.

We recall the basic equations. The electron and hole components $`u(𝐫)`$ and $`v(𝐫)`$ of the spinor wave function satisfy the BdG equation

$$\left(\begin{array}{cc}H& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& H\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)=E\left(\begin{array}{c}u\\ v\end{array}\right),$$
(2)
which contains the single-particle Hamiltonian $`H=^2+V(𝐫)E_\mathrm{F}`$ (with confining potential $`V`$) and the pair potential $`\mathrm{\Delta }(𝐫)`$ (vanishing in the normal metal and equal to $`\mathrm{\Delta }_0`$ in the superconductor). The energy $`E`$ is measured relative to the Fermi energy $`E_\mathrm{F}=k_\mathrm{F}^2`$, in units such that $`\mathrm{}^2/2m1`$. (In these units the mean level spacing $`\delta `$ is related to the area $`𝒜`$ of the billiard by $`\delta =4\pi /𝒜`$.) We assume that the motion inside the billiard is ballistic ($`V=0`$) and that the interface with the superconductor is a ballistic point contact of width $`W\lambda _\mathrm{F}`$ (so that the number of channels $`N=2W/\lambda _\mathrm{F}1`$ and the Thouless energy $`E_\mathrm{c}=N\delta \delta `$). We work in the regime $`\mathrm{\Delta }_0\mathrm{}v_\mathrm{F}/W`$ (which also implies $`\mathrm{\Delta }_0E_\mathrm{c}`$), to ensure that the excitation spectrum is independent of the properties of the superconductor.

For a semiclassical description one substitutes $`(u,v)=(\overline{u},\overline{v})Ae^{iS}`$, with $`\mathrm{}S`$ the action along a classical trajectory at the Fermi energy. The wave amplitude $`A`$ is related to the classical action by the continuity equation $`(A^2S)=0`$, while $`S`$ itself satisfies the Hamilton-Jacobi equation $`|S|^2=E_\mathrm{F}V`$ (so that $`\mathrm{}S`$ is the momentum along the trajectory). The BdG equation takes the form

$$\left(\begin{array}{cc}2ik_\mathrm{F}_s+\delta H& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& 2ik_\mathrm{F}_s\delta H\end{array}\right)\left(\begin{array}{c}\overline{u}\\ \overline{v}\end{array}\right)=E\left(\begin{array}{c}\overline{u}\\ \overline{v}\end{array}\right),$$
(3)
with $`\delta H\overline{u}=A^1^2(A\overline{u})`$. The derivative $`_s=k_\mathrm{F}^1(S)`$ is taken along the classical trajectory. The Andreev approximation consists in neglecting the term $`\delta H`$ containing second derivatives of the slowly varying functions $`A,\overline{u},\overline{v}`$.

We consider a classical trajectory that starts as an electron at a point $`q(0,W)`$ along the interface with the superconductor, making an angle $`\varphi (\pi /2,\pi /2)`$ with the normal (cf. Fig. 1). The product $`b=q\mathrm{cos}\varphi `$ is the “impact parameter”. The trajectory returns to the interface after a path length $`\mathrm{}`$, and then it is retraced in the opposite direction as a hole . The coordinate $`s(0,\mathrm{})`$ runs along one repetition of this trajectory. We count trajectories with measure $`dqd\mathrm{sin}\varphi =dbd\varphi `$, corresponding to a uniform measure in phase space. Equivalently, we can sum over scattering channels $`n=1,2,\mathrm{}N`$, related to $`\varphi `$ by $`nN|\mathrm{sin}\varphi |`$.

If we ignore the term $`\delta H`$ in Eq. (3) we recover the Bohr-Sommerfeld density of states Mel96 ; Lod98 ; Sch99 . Indeed, without $`\delta H`$ the solution of the eigenvalue problem is

$$\left(\begin{array}{c}\overline{u}_m\\ \overline{v}_m\end{array}\right)=\frac{1}{\sqrt{2\mathrm{}}}\left(\begin{array}{c}e^{im\pi s/2\mathrm{}}\\ ie^{im\pi s/2\mathrm{}}\end{array}\right),E_m=m\pi k_\mathrm{F}/\mathrm{},$$
(4)
with $`m=\pm 1,\pm 3,\pm 5\mathrm{}`$ running over positive and negative odd integers. The path length $`\mathrm{}`$ in a chaotic billiard varies in a quasi-random way upon varying the initial conditions $`q`$ and $`\varphi `$, with an exponential distribution $`P(\mathrm{})=\overline{\mathrm{}}^1\mathrm{exp}(\mathrm{}/\overline{\mathrm{}})`$. (The mean path length is $`\overline{\mathrm{}}=4\pi k_\mathrm{F}/E_\mathrm{c}`$.) The density of states

$$\rho (E)=N_0^{\mathrm{}}𝑑\mathrm{}P(\mathrm{})\underset{m=1,3,5}{\overset{\mathrm{}}{}}\delta [EE_n(\mathrm{})]$$
(5)
then evaluates to the $`\rho _{\mathrm{BS}}`$ of Eq. (1).

The key assumption that will enable us to go beyond the Andreev and Bohr-Sommerfeld approximations is to assume that the amplitude $`A`$ varies more slowly in space than the spinor components $`\overline{u}`$ and $`\overline{v}`$, so that we can approximate $`\delta H`$ by $`^2`$ (neglecting derivatives of $`A`$). The characteristic length scale $`L_A`$ for the spatial dependence of $`A`$ is set by the smoothness of the confining potential $`V`$, while the characteristic length scale for $`\overline{u},\overline{v}`$ is the contact width $`W`$. By assuming $`L_AW`$ we consider the case that diffraction occurs predominantly at the interface with the superconductor, rather than inside the billiard. Since $`A`$ depends on the shape of the billiard, this is the regime in which we can hope to obtain a geometry-independent “universal” result.

To investigate the effect of $`\delta H`$ we restrict the dimensionality of the Hilbert space in two ways: Firstly, we neglect any mixing of the $`N`$ scattering channels. (This is known to be a good approximation of the diffraction that occurs when a narrow constriction opens abruptly into a wide region Sza89 ; it does not require smooth corners in the contact.) Secondly, since we are interested in excitation energies $`EE_\mathrm{c}`$, we include only the two lowest eigenstates $`m=\pm 1`$ of the zeroth-order solution (4). \[The contributions from higher levels are smaller by a factor $`\mathrm{exp}(E_\mathrm{c}/2E)`$.\] We need to include both $`E_1`$ and $`E_1`$, although the excitation spectrum contains only positive eigenvalues, because of the (virtual) transitions between these two levels induced by $`\delta H`$. With these restrictions we have for each scattering channel a one-dimensional eigenvalue problem. The effective Hamiltonian $`_{\mathrm{eff}}`$ is a $`2\times 2`$ matrix differential operator acting on functions of $`b`$.

We write $`_{\mathrm{eff}}=_0+_1`$, where $`_0`$ corresponds to the Andreev approximation and $`_1`$ contains the diffractive effects. The zeroth-order term is diagonal,

$$_0=\left(\begin{array}{cc}\pi k_\mathrm{F}/\mathrm{}(b)& 0\\ 0& \pi k_\mathrm{F}/\mathrm{}(b)\end{array}\right).$$
(6)
The relation between $`\mathrm{}`$ and $`b`$ is determined by the differential equation $`d\mathrm{}/db=g(b)\mathrm{exp}(\kappa \mathrm{})`$, which expresses the exponential divergence of nearby trajectories (in terms of a Lyapunov exponent $`\kappa =\alpha /v_\mathrm{F}`$ given as inverse length rather than inverse time). The pre-exponential $`g(b)`$ is of order unity, changing sign at extrema of $`\mathrm{}(b)`$. Upon integration one obtains

$$\kappa \mathrm{}(b)=\mathrm{ln}|\kappa b|+𝒪(1),$$
(7)
where we have shifted the origin of $`b`$ such that $`b=0`$ corresponds to a local maximum $`\mathrm{}_{\mathrm{max}}\overline{\mathrm{}}`$ of $`\mathrm{}(b)`$. \[The logarithmic singularity is cut off at $`|\kappa b|\mathrm{exp}(\kappa \mathrm{}_{\mathrm{max}})`$.\] There is an exponentially large number $`𝒩(\mathrm{})\mathrm{exp}(\kappa \mathrm{}\mathrm{}/\overline{\mathrm{}})`$ of peaks around which Eq. (7) applies.

To obtain the diffractive correction $`_1`$, in the regime that $`\delta H=^2`$, we express the Lapacian in the local displacements $`ds`$ and $`d\zeta =e^{\kappa s}db`$. Since these are approximately orthogonal for $`\kappa s1`$, we have

$$\delta H=_s^2e^{2\kappa s}_b^2.$$
(8)
The first term $`_s^2`$ is a small correction to the zeroth order density of states. The second term $`_b^2`$, in contrast, is a singular perturbation because it associates a kinetic energy with the variable $`b`$. The resulting zero-point motion implies a non-zero ground state energy, and hence it is responsible for the opening of an excitation gap. Projecting $`_1`$ onto the space spanned by the two lowest eigenfunctions $`n=\pm 1`$ of Eq. (4), and retaining only the leading order terms in $`1/\kappa \mathrm{}`$, we find

$$_1=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)\frac{\pi }{4}\frac{d}{db}[\kappa \mathrm{}(b)]^2\frac{d}{db}+𝒪(\kappa \mathrm{})^3.$$
(9)
The effective Hamiltonian can be brought into a more familiar form by the unitary transformation $`_{\mathrm{eff}}e^{i\sigma _1\pi /4}_{\mathrm{eff}}e^{i\sigma _1\pi /4}`$ (with $`\sigma _i`$ a Pauli matrix), followed by the change of variable $`x=\kappa b\kappa b\mathrm{ln}|\kappa b|`$ (in the range $`|x|<1`$). We work again to leading order in $`1/\kappa \mathrm{}`$, and find

$$_{\mathrm{eff}}=\pi k_\mathrm{F}\kappa \left(\begin{array}{cc}ϵ_x^2& i/\mathrm{ln}|x|\\ i/\mathrm{ln}|x|& ϵ_x^2\end{array}\right),ϵ\frac{\kappa }{4k_\mathrm{F}}.$$
(10)
This effective Hamiltonian has the same form as the BdG Hamiltonian (2), for a fictitious one-dimensional system having $`V=E_\mathrm{F}`$ and having a pair potential $`\mathrm{\Delta }(x)`$ that vanishes logarithmically $`1/\mathrm{ln}|x|`$ at the origin (cf. Fig. 2). The kinetic energy $`ϵ_x^2`$ gives a finite excitation gap, even though $`ϵ1`$. Let us now compute this gap.

Since $`e^{i\sigma _2\pi /4}_{\mathrm{eff}}^2e^{i\sigma _2\pi /4}`$ is a diagonal matrix, the spectrum of $`_{\mathrm{eff}}`$ is given by the scalar eigenvalue problem

$$\left|ϵ\frac{d^2}{dx^2}+\frac{i}{\mathrm{ln}|x|}\right|^2\mathrm{\Psi }(x)=\left(\frac{E}{\pi k_\mathrm{F}\kappa }\right)^2\mathrm{\Psi }(x).$$
(11)
The ground state energy is the excitation gap $`E_{\mathrm{gap}}`$. To generate an asymptotic expansion of $`E_{\mathrm{gap}}`$ for small $`ϵ`$, we first multiply both sides of Eq. (11) by a factor $`Z^2`$ and then substitute $`x=X\sqrt{ϵZ}`$. This results in

$`\left|{\displaystyle \frac{d^2}{dX^2}}+iU\right|^2\mathrm{\Psi }=\left({\displaystyle \frac{ZE}{\pi k_\mathrm{F}\kappa }}\right)^2\mathrm{\Psi },`$ (12)
$`U(X)={\displaystyle \frac{2Z}{\mathrm{ln}ϵZ}}\left[1{\displaystyle \frac{2\mathrm{ln}|X|}{\mathrm{ln}ϵZ}}+𝒪(\mathrm{ln}ϵZ)^2\right].`$ (13)
We now choose $`Z`$ such that $`Z^2=\mathrm{ln}^3ϵZ`$ and obtain the biharmonic equation

$`(d^4/dX^4+16\mathrm{ln}|X|)\mathrm{\Psi }=\omega \mathrm{\Psi },`$ (14)
$`\omega =(ZE/\pi k_\mathrm{F}\kappa )^24Z^{2/3}+𝒪(Z^{1/3}).`$ (15)
The ground state of Eq. (14) is at $`\omega _0=14.5`$. Substituting in Eq. (15), and using $`Z^{2/3}=|\mathrm{ln}ϵ|\frac{3}{2}\mathrm{ln}|\mathrm{ln}ϵ|+𝒪(1/\mathrm{ln}ϵ)`$, we arrive at

$$E_{\mathrm{gap}}=\frac{2\pi k_\mathrm{F}\kappa }{|\mathrm{ln}ϵ|}\left(1+\frac{3\mathrm{ln}|\mathrm{ln}ϵ|}{2|\mathrm{ln}ϵ|}+\frac{\omega _0}{8|\mathrm{ln}ϵ|}+𝒪(\mathrm{ln}ϵ)^{3/2}\right).$$
(16)
Only the leading order term is significant in view of the approximations made in Eq. (10). Restoring the original variables we have

$$E_{\mathrm{gap}}=\frac{\pi \mathrm{}\alpha }{\mathrm{ln}(v_\mathrm{F}/\alpha \lambda _\mathrm{F})}.$$
(17)
The Ehrenfest time $`\tau _\mathrm{E}=\alpha ^1\mathrm{ln}(L/\lambda _\mathrm{F})`$ contains the classical length $`L=v_\mathrm{F}/\alpha `$, which is of the order of the linear dimension of the billiard.

The density of states $`\rho (\omega )`$ of the biharmonic equation (14) can be calculated numerically Gre97 . The density of states $`\rho (E)`$ near the gap is related to $`\rho (\omega )`$ by

$$\rho (E)=\frac{8𝒩|\mathrm{ln}ϵ|}{E_{\mathrm{gap}}}\rho [\omega =\omega _0+8|\mathrm{ln}ϵ|(E/E_{\mathrm{gap}}1)],$$
(18)
and is plotted in Fig. 2 for $`|\mathrm{ln}ϵ|=10`$. The factor $`𝒩\mathrm{exp}(\pi k_\mathrm{F}\kappa /E_{\mathrm{gap}}E_\mathrm{c}/E_{\mathrm{gap}})`$ counts the number of peaks in $`\mathrm{}(b)`$ around which $`_{\mathrm{eff}}`$ applies. The Bohr-Sommerfeld approximation (1) corresponds to the large-$`\omega `$ asymptote $`\rho (\omega )=\frac{1}{16}\mathrm{exp}(\omega /16)`$. Since $`\omega \omega _01`$ implies $`E/E_{\mathrm{gap}}11/|\mathrm{ln}ϵ|`$, the width $`\mathrm{\Delta }EE_{\mathrm{gap}}/|\mathrm{ln}ϵ|`$ of the energy range above the gap in which the Bohr-Sommerfeld approximation breaks down is small compared to the gap itself.

Because $`_{\mathrm{eff}}`$ has only a few levels in the range $`\mathrm{\Delta }E`$, the density of states $`\rho (E)`$ oscillates strongly in this range. These levels are highly degenerate (by a factor $`𝒩`$) in our approximation. Tunneling between the levels will remove the degeneracy and smooth the oscillations. (A small amount of smoothing has been inserted by hand in Fig. 2.) We surmise that some oscillatory energy dependence will remain, but this is an aspect of the problem that needs further study.

In conclusion, we have analyzed a mechanism for the “gap phenomenon” in the proximity effect of chaotic systems. Diffraction at the contact with the superconductor is described by an effective Hamiltonian $`_{\mathrm{eff}}`$ that contains (1) a kinetic energy which vanishes in the classical limit and (2) a pair potential with a logarithmic profile. The resulting excitation gap $`E_{\mathrm{gap}}`$ (being the ground state energy of $`_{\mathrm{eff}}`$) vanishes logarithmically as the ratio of the Fermi wavelength and a classical length scale (set by the Lyapunov exponent) goes to zero. The time scale $`\mathrm{}/E_{\mathrm{gap}}`$ is the Ehrenfest time, providing a manifestation of quantum chaos in the superconducting proximity effect.

This work was supported by the Dutch Science Foundation NWO/FOM. We thank A. N. Morozov and P. G. Silvestrov for helpful discussions.


text:
# How Cumbersome is a Tenth Order Polynomial?: The Case of Gravitational Triple Lens Equation

## 1 Gravitational Triple Lens Equation

The lens equation of a set of three gravitationally bound point masses is written with three (real) relative mass parameters, $`ϵ_1`$, $`ϵ_2`$, and $`ϵ_3`$, and three complex parameters for their 2-dimensional positions. If $`z`$ is the position of an image and $`\omega `$ is the position of its source, the lens equation reads as follows.

$$\omega =z\frac{ϵ_1}{\overline{z}\overline{x}_1}\frac{ϵ_2}{\overline{z}\overline{x}_2}\frac{ϵ_3}{\overline{z}\overline{x}_3}zf(\overline{z};\overline{x}_j)$$
(1)
The mass parameters are subject to a constraint $`ϵ_1+ϵ_2+ϵ_3=1`$ where $`1=M`$ is the total mass, and the position variables include coordinate degrees of freedom: two degrees of freedom for translation and one degree of freedom for roation in the two-dimensional lens plane. We choose a coordinate system such that the position of a mass, $`x_1`$, and the center of mass of the other two elements, $`x_4`$, define the lens axis along the real axis of the complex plane.

$$(ϵ_2+ϵ_3)x_4=ϵ_2x_2+ϵ_3x_3$$
(2)
Then, $`x_1`$ and $`x_4`$ are real, and the triple lens system is completely specified by a set of five parameters.

$$\mathrm{}=|x_1x_4|;\mathrm{}^{}=|x_2x_3|;\gamma =[0,\pi ];ϵ_2=(0,1);ϵ_3=(0,1)$$
(3)
The angle $`\gamma `$ is the angle between $`(x_2x_3)`$ and the lens axis. In order to find the images of a given source, we need to solve the lens equation.

When the lensing system and a source star are far apart, there are four images of the source star, which we can verify easily for a source at $`\omega =\mathrm{}`$. From the lens equation (1), we find that the images of the source at $`\omega =\mathrm{}`$ are at the position of the source, $`z=\mathrm{}`$, and at the three lens positions, $`z=x_1,x_2`$, and $`x_3`$. The magnifications of the images at the lens positions are $`0`$, and they are hardly images in any practical sense because there are no photon fluxes related to the image positions. In other words, there is only one physical image for a source at $`\mathrm{}`$ – namely, the image of the unmagnified source star. However, if we move the source toward the lens system, the images at the lens positions move away from the lens positions, and they do have non-vanishing photon fluxes. Thus, we freely speak of image positions with zero magnification as the continuity limit so that we don’t have to cut out the lens positions from the image space. In practice, $`\mathrm{}`$ is a large distance limit, which may be no more than a few hundred au in any physical relevance of Galactic lensing.

If we consider a source trajectory away from the caustic regions, the images form four smooth curves anchored at the base points of $`z=x_1,x_2,x_3`$, and $`\mathrm{}`$ which are the image positions of the source at $`\omega =\mathrm{}`$. The complex lens plane can be considered a large two sphere with one point at $`\mathrm{}`$, and then, the four smooth disjoint image curves are loops with fixed points at the base points. The image curves can be calculated esaily from the base points using Newtonian transportation method (the algorithm can be found in Numerical Recipes). However, the physical interest of lensing lies in caustic regions where the lensing signals are most obvious, and this makes the Newtonian method useless in microlensing business. The large lensing signals arise where the inverse Jacobian determinant of the lens equation is large, and the large Jacobian determinant makes Newtonian transportation method unstable and unusable. The method of infinitesimal corrections based on derivative values is valid where the Jacobian determinant is practically finite and the mapping is non-degenerate. Laguerre’s method (see Numerical Recipes) is based on complex analysis and uses first and second derivatives of the polynomial to help converge to the solutions from initial values which one is free to guess.

Another complication of the caustic region is that the size (which may be $`110\mu `$as) of the source star can not be ignored when the source star crosses a caustic curve. In fact, the luminosity profile of the source star can manifestly affect the shape of the light curve, and the luminosity profile dependence of the light curve shape during a line caustic crossing lasts $`3`$ stellar radius crossing times (Rhie and Bennett, 1999). In order to incorporate the finite size effects, one can locally pixelize the image plane ($`z`$-plane) in the neighborhood of the critical curve near the images and count up the image pixels that are mapped into the source disk (Bennett and Rhie, 1996). The luminosity profile of the source can be easily incorporated by weighting each image pixel by the luminosity shape function value at the corresponding source position. The pixel size can be adjusted for a desired resolution. In the ray shooting method which has been designed mainly to handle a system of a large number of lensing elements solving whose lens equation is impractical if not impossible, the entire lens plane is pre-pixelized and each pixel in the image plane is tested whether the center of the pixel is mapped into the area defined by the source disk in the source plane. The ray shooting method is inefficient for fitting light curves of low-multiplicity point mass lenses as microlensing planet systems.

The astro-ph version of Gaudi, Naber, and Sackett (1998) conveys an impression that solving the lens equation by finding the roots of the 10-th order polynomial analytic equation is slower than using the ray shooting method for triple lens light curve fitting. However, Gaudi (private communication) recently informed us that the comment on the calculation speed was a comparison between binary lens and triple lens but not a comparison between root finding method and ray shooting method for the triple lens equation.

## 2 The Tenth Order Polynomial Equation

The lens equation is an explicit function from an image position $`z`$ to its source position $`\omega `$, and $`\omega (z,\overline{z})`$ is a genuine real function, namely a function of both $`z`$ and $`\overline{z}`$. However, the $`z`$-depdence is linear, and this simplicity is behind the analyticity of the differential behavior which can be completely described by one analytic function $`\kappa _z\overline{\omega }`$. The linearity in $`z`$ also makes it easy to find an analytic equation where the lens equation is embedded. Using $`z\omega =f(\overline{z};\overline{x}_j)`$ and $`\overline{z}\overline{\omega }=f(z;x_j)`$,

$$z\omega =f(f(z;x_j)+\overline{\omega };\overline{x}_j)$$
(4)
If we let $`Hz_1z_2z_3`$ and $`Gϵ_1z_2z_3+ϵ_2z_3z_1+ϵ_3z_1z_2=_{\mathrm{cyc}}ϵ_iz_jz_k`$, then $`f=G/H`$, and it is simple to see that equation (4) is a tenth order polynomial equation. If we let $`\overline{\omega }_j\overline{\omega }\overline{x}_j`$,

$$0=(z\omega )(G+\overline{\omega }_1H)(G+\overline{\omega }_2H)(G+\overline{\omega }_3H)\underset{\mathrm{cyc}}{}Hϵ_i(G+\overline{\omega }_jH)(G+\overline{\omega }_kH)$$
(5)
An analytic polynomial equation has the same number of solutions as the order (see any textbook on complex variable or mathematical physics), there can be up to ten images in a triple lensing. There are only four images for $`\omega =\mathrm{}`$, the number of images changes by two at a caustic crossing, and the caustic curves form heierachical structures of domains for high multiplicity images. There are triple lenses with domain $`𝒟^3`$ the sources therein produce ten images, the maximum possible number of images (Rhie, 1997). As we repeatedly emphasized in Rhie (1997), the triple lens equation is equivalent to the tenth order polynomial analytic equation only in the domain $`𝒟^3`$, and the statement on the equivalence in section 3 of Gaudi, Naber, and Sackett (1998) seems to be a misquote which can mislead the readers to think that the two equations are equivalent everywhere.

In order to fit a triple microlensing light curve, we need to solve the tenth order analytic equation (and select the image solutions that satisfy the lens equation). This can be done numerically using root finders available in the literature, and we only need to type in the coefficients. The coefficients may appear to be cumbersome as declared in Gaudi, Naber, and Sackett (1998), and it is indeed the case if we, for example, calculate the coeffients using Mathematica. The output of an algebraic computing package is (unnecessarily) messy even for the binary lens equation. Thus, it is useful to group (or not to unfold) the coefficients, which is a natural intermediate process in hand calculations.

In the center of mass system, $`ϵ_jx_j=0`$, $`H`$ and $`G`$ are written with four coefficient functions, $`a,b,c`$, and $`d`$: $`H=z^3+az^2+bz+c`$ and $`G=z^2+az+d`$, where $`a(x_1+x_2+x_3),bx_1x_2+x_1x_3+x_2x_3,cx_1x_2x_3`$, and $`d_{\mathrm{cyc}}ϵ_ix_jx_k`$. If we let $`a_\omega \overline{\omega }_1+\overline{\omega }_2+\overline{\omega }_3`$, $`b_\omega \overline{\omega }_1\overline{\omega }_2+\overline{\omega }_2\overline{\omega }_3+\overline{\omega }_3\overline{\omega }_1`$, $`c_\omega \overline{\omega }_1\overline{\omega }_2\overline{\omega }_3`$, and $`d_\omega _{\mathrm{cyc}}ϵ_i\overline{\omega }_j\overline{\omega }_k`$, then equation (5) becomes

$$0=G^3(z\omega )+G^2H((z\omega )a_\omega +1)+GH^2((z\omega )b_\omega +a_\omega \overline{\omega })+H^3((z\omega )c_\omega +d_\omega )$$
(6)
If we let $`G^3H_{0k}z^k(k6)`$, $`G^2HH_{1k}z^k(k7)`$, $`GH^2H_{2k}z^k(k8)`$, and $`H^3H_{3k}z^k(k9)`$, the equation becomes

$$0=\underset{k=1}{\overset{10}{}}\mathrm{cff}(k)z^k$$
(7)
where the polynimial coefficients are

* $`\mathrm{cff}(k)=\left(H_{0k1}+H_{1k1}a_\omega +H_{2k1}b_\omega +H_{3k1}c_\omega \right)\left(H_{0k}\omega +H_{1k}(\omega a_\omega 1)+H_{2k}(\omega b_\omega +a_\omega \overline{\omega })+H_{3k}(\omega c_\omega +b_\omega )\right)`$

The coefficients $`H_{ij}`$ are polynomials of $`a,b,c`$, and $`d`$ where the polynomial coefficients are simple combinatoric integers.

* $`H_{39}=1;H_{38}=3a;H_{37}=3b+3a^2;H_{36}=3c+6ab+a^3;H_{35}=6ac+3b^2+3a^2b;H_{34}=6bc+3a^2c+3ab^2;H_{33}=3c^2+6abc+b^3;H_{32}=3ac^2+3b^2c;H_{31}=3bc^2;H_{30}=c^3`$.
* $`H_{28}=1;H_{27}=3a;H_{26}=d+2b+3a^2;H_{25}=2ad+4ab+a^3+2c;H_{24}=2db+da^2+4ac+2a^2b+b^2;H_{23}=2dc+2dab+2a^2c+ab^2+2bc;H_{22}=2cad+db^2+2abc+c^2;H_{21}=2bcd+ac^2;H_{20}=c^2d`$
* $`H_{17}=1;H_{16}=3a;H_{15}=2d+3a^2+b;H_{14}=4ad+a^3+2ab+c;H_{13}=d^2+2a^2d+2bd+ba^2+2ac;H_{12}=ad^2+2abd+2cd+ca^2;H_{11}=bd^2+2acd;H_{10}=cd^2`$
* $`H_{06}=1;H_{05}=3a;H_{04}=3d+3a^2;H_{03}=6ad+a^3;H_{02}=3d^2+3a^2d^2;H_{01}=3ad^2;H_{00}=d^3`$

### 2.1 Comments

An interested party may download the source file of this manuscript to avoid tying the coefficients. We encourage to check the coefficients, however. It is a quick exercise once one adopts the poor person’s calculation with pencil and paper as shown above; also with the free biocomputer which is harder to hack either internally or externally barring the long term process of brainwashing. We also found it useful to test the symmetric cases in rh97 whose image solutions behavior (for example, the number of images) is known. The critical curve is obtained by solving $`\kappa =e^{i2\phi }`$ which is an eighth order polynomial equation. The caustic curve is obtained by applying the lens equation to the critical curve solution. It is fine to use $`\phi =[0,2\pi )`$ as the parameter for equal interval sampling. Let $`\delta =2\pi /N`$ for $`N`$ not too large and observe the intervals (or speeds) of the solutions on the critical curve and the caustic curve. Note especially the density of the solutions around the cusps of the caustic curve. It is worth pausing for a moment counting the relative numbers of the solutions on the stellar caustic and planetary caustics for planet systems lenses.

For the measurements and theoretical calculations of the anomalous magnetic moment of the electron, we have consulted Peskin and Schroeder (1999) and Kinoshita (1990). We have considered drawing all the Feynman diagrams to lay out the degree of lengthy squiggly messiness but given up, and our misadventure may be considered an indirect testimony of the degree of mud wrestling the exquisite anomalous magnetic moment requires.

Barring the notion that lensing community may have been chosen to be dealt with laxed scrutiny of nature, we have no doubt that we have no luxury to complain about a bit of algebra we encounter in lensing. In fact, we find it a cherished treasure that low multiplicity point lenses are exactly solvable and their light curves can be reconstructed and interpreted without ambiguities. Microlensing events do share the transiency with the scattering events in accelerator particle physics. The both need high resolution data for minute rare prized signals and complete data to interpret them. The both rely on methodical analyses building from the simpler and dominant events to more rare events. Thus, it is important to have a homogeneous and comprehensive data set of microlensing events where consistencies can be tested within in order to find microlensing planets. The best bet for such data set is microlensing from space. Accelerator particle physics of the last century is at the foundation of the Standard Model or the Theory of Matter. We expect that space microlensing planet search will lay a foundation of extrasolar planet physics within a few years of operation of a small space telescope (Bennett and Rhie, 2000). In the comprehensive data set, the so-called high magnification events with only stellar caustic signals will form an independent data set that can be used for a consistency check of the interpretation of the planetary light curves as a subset.

This note is based on the work with D. Bennett for Bennett et al. (1999).