text:
# Sub-Subgiants in the Old Open Cluster M67?

## 1 Introduction

The open cluster M 67 is one of the most comprehensively studied of all star clusters, and has long been the prototype for old (4 Gyr) open clusters in the Galaxy. Indeed, one of the first photoelectric color-magnitude diagrams was derived for M 67 by Johnson & Sandage (1955). Since that time the precision of stellar photometry has steadily improved, and with it has the definition of the M 67 giant branch (Janes & Smith 1984, Montgomery et al. 1993). Indeed the remarkable narrowness of the M 67 giant branch has served as a precise touchstone against which innumerable single-star evolution models have been tested (e.g., Dinescu et al. 1995).

However, even after application of strict proper-motion membership criteria, the color-magnitude diagram of M 67 remains littered with stars that do not fall on a single-star isochrone (Fig. 1). Some of these seeming anomalies can be accounted for. For example, one luminous star to the blue of the giant branch is a spectroscopic binary whose composite light can be explained by a giant–main-sequence pairing; another is a giant–white dwarf pair with a complicated history (Mathieu et al. 1990, Verbunt & Phinney 1995, Landsman et al. 1997). Many of the stars immediately above the turnoff and subgiant branch are spectroscopic binaries and consequently overluminous. And the parallel sequence of stars to the red of the main sequence is assuredly comprised of binaries as well (Montgomery et al. 1993).

Still, there remain stars that are not so easily explained. Most famous are the blue stragglers, first noted in M 67 by Johnson & Sandage. The origin of blue stragglers in an open cluster environment is still not securely understood today (Bailyn 1995, Hurley et al. 2001). Similarly, the yellow giant S 1072, lying 1.2 magnitudes in V above the main-sequence turnoff, has defied explanation. A kinematic member in all three dimensions and lying in projection in the cluster core, its explanation as a non-member remains possible but not satisfying (Mathieu & Latham 1986, Nissen et al. 1987). In this paper we consider two more stars – S 1063 and S 1113 – which by every indication are cluster members, yet whose location in the M 67 color-magnitude diagram is dramatically inconsistent with single-star evolutionary theory. Specifically, as shown in Fig. 1, these stars lie below the subgiant branch.

Our attention has converged on these stars from two directions. We have underway a multi-decade survey of the spectroscopic binary population in M 67 (Mathieu et al. 1990, Latham et al. 1997). Both S 1063 and S 1113 were found to be spectroscopic binaries, with S 1113 in particular being notable for the rapid rotation of its primary star. Belloni et al. (1998) have also undertaken a comprehensive study of stellar X-ray sources in M 67 based on ROSAT observations. S 1063 and S 1113 were independently discovered to be X-ray sources, of which there are only 25 known among cluster members. Both stars are also photometric variables (van den Berg et al. 2002), which are infrequent in a cluster of this age (Stassun et al. 2002).

In this paper we present a comprehensive discussion of the orbital, spectroscopic, and photometric properties of S 1063 and S 1113. Regrettably, like the blue stragglers and S 1072, their interpretation remains a puzzle.

## 2 The stars

The stars S 1063 and S 1113 had attracted some attention prior to the radial-velocity and X-ray studies. Both stars are included in the highly precise proper-motion study of Girard et al. (1989), in which S 1063 is given a proper-motion membership probability of 97% and S 1113 of 98%. Being a bit less than one core radius from the cluster center, S 1063 has been included in most photometric studies of M 67; S 1113 is located at three core radii, and thus only has been observed in the wider-field photometric surveys. Racine (1971) noted that S 1063 was photometrically unusual in that it lies roughly a magnitude in $`V`$ below the subgiant branch. S 1113 has photometric properties very similar to those of S 1063, but seems to have escaped notice until mentioned by Kaluzny & Radczynska (1991). In Fig. 1 we show the color-magnitude diagram of M 67 derived from the photometry of Montgomery et al. (1993); S 1063 and S 1113 are indicated with boxes.

Both stars have been found to be photometric variables. Racine (1971) noted a large range (0.18 mag) in photometric observations of S1063, and its variability was subsequently confirmed by other observers (Rajamohan et al. 1988, Kaluzny & Radczynska 1991). S 1113 has a variable star designation – AG Cnc – and Kaluzny & Radczynska (1991) found S 1113 to be photometrically variable at the 0.05 mag level. Based on the variability (but not light curves) of the stars and their location in the color-magnitude diagram, Kaluzny & Radczynska (1991) suggested that both stars are highly evolved W UMa binaries with extremely small mass ratios and consequently small photometric amplitudes.

Both S 1063 and S 1113 were included in a large survey for spectroscopic binaries among M 67 proper-motion cluster members begun in the mid-1980’s (last summarized in Latham et al. 1997). S 1063 was found to be a single-lined spectroscopic binary, with no evident distinction otherwise. S 1113 was immediately recognized as an unusual double-lined spectroscopic binary in that the primary star was a rapid rotator.

Independently, ROSAT PSPC observations of M 67 identified these two stars as among 25 cluster members detected to have X-ray emission (Belloni et al. 1993, 1998). As members of this select group, S 1063 and S 1113 received optical spectroscopic attention by Pasquini & Belloni (1998) and van den Berg et al. (1999). Both stars show strong emission cores in the Ca II H and K lines, indicative of chromospherically active stars. Both stars also show H$`\alpha `$ emission. In S1063 the H$`\alpha `$ line is asymmetric, showing emission which is blue-shifted with respect to an absorption feature. The position of the absorption line agrees with the velocity of the primary, perhaps suggesting that the absorption can be attributed to the primary star and that the H$`\alpha `$ emission originates elsewhere in the system. Alternatively, the H$`\alpha `$ line profile is very similar to that seen from HK Lac by Catalano & Frasca (1994), who attribute the emission to a large flare lasting 6 days. In S 1113 the H$`\alpha `$ emission is broad and clearly centered on the velocity of the primary star, with an equivalent width of 15 Å.

## 3 Observations and data analyses

### 3.1 Speedometry and orbital solutions

Radial-velocity observations of S 1063 were begun in 1987, and 28 observations were obtained through 1989. 18 radial-velocity observations were obtained for S 1113 between 1989 and 1998. All observations were obtained with the Center for Astrophysics (CfA) Digital Speedometers (Latham 1992). Two nearly identical instruments were used on the Multiple Mirror Telescope 111Some of the observations reported here were obtained with the Multiple Mirror Telescope, a joint facility of the Smithsonian Institution and the University of Arizona. and the 1.5-m Tillinghast Reflector atop Mt. Hopkins, Arizona. Echelle spectrographs were used with intensified photon-counting Reticon detectors to record about 45 Å of spectrum in a single order near 5187 Å with a resolution of 8.3 km s-1 and signal-to-noise ratios ranging from 8 to 15 per resolution element.

The 28 radial velocities for S 1063 were measured using the one-dimensional correlation package rvsao (Kurtz & Mink 1998) running inside the IRAF 222 IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc., under contract with the National Science Foundation. environment. The template spectrum was drawn from a new grid of synthetic spectra (Morse & Kurucz, in preparation) calculated using Kurucz’s code ATLAS9. We correlated our observed spectra against an array of solar-metallicity templates spanning effective temperature, surface gravity, and projected rotation velocity. The template giving the highest peak correlation (averaged over all observations) had $`\mathrm{log}g=3.5`$, $`T_{\mathrm{eff}}=5000`$ K, and $`v\mathrm{sin}i`$ = 6 km s-1. The radial velocities derived with this template are presented in Table 1. The average internal error of these observations is 0.86 km s-1.

A single-lined orbital solution with period $`P=18.396`$ days is easily derived from these data; the elements are given in Table 2 and the theoretical velocity curves are superimposed on the data in Fig. 2. The center-of-mass velocity of $`\gamma `$ = 34.3 $`\pm `$ 0.2 km s-1 is consistent with the cluster mean radial velocity of 33.5 km s-1 and observed velocity dispersion of 0.5 km s-1 (Mathieu 1983), so the star is a kinematic cluster member in all three dimensions. Otherwise, the orbital elements are unremarkable.

The analyses of the S 1113 spectra were done using the two-dimensional correlation technique TODCOR (Zucker & Mazeh (1994), as implemented at CfA). A two-dimensional array of templates in effective temperature and projected rotation velocity, again derived from Kurucz spectra, were correlated with all of the S 1113 spectra. The effective temperatures and rotation velocities of the templates were selected so as to maximize the correlation peak heights, averaged over all 18 spectra and weighted according to signal-to-noise ratio. We find projected rotational velocities $`v\mathrm{sin}i`$ of 53 km s-1 for the primary and 11 km s-1 for the secondary, respectively, with a 1-$`\sigma `$ uncertainty of 2 km s-1 on each. Van den Berg et al. (1999) previously measured $`v\mathrm{sin}i`$ for both stars, finding 45$`\pm `$6 km s-1 for the primary and 12$`\pm `$1 km s-1 for the secondary. Our derived effective temperatures are 4800 K and 5500 K for the primary and secondary, respectively, with a 1-$`\sigma `$ uncertainty of 150 K. Finally, we derive a luminosity ratio at 5187Å of 0.35$`\pm `$0.02 (secondary/primary). With color corrections we find a luminosity ratio at V of 0.32 or 1.23 mag.

The primary and secondary radial velocities derived with these template parameters are given in Table 3. We obtained a double-lined orbital solution using equal weighting for all velocities. If the eccentricity is allowed to vary, the best-fit value is $`e=0.022\pm 0.010`$. As this measure is not significantly different from zero, we present elements for a circular (e=0) orbital solution in Table 4 and show the orbit curves superimposed on the phased data in Fig. 3. (In these elements, $`T_0`$ represents a time of maximum radial velocity.) Again, the center-of-mass velocity of $`\gamma `$ = 33.4$`\pm `$0.4 km s-1 is consistent with the cluster mean radial velocity. The mass ratio M2/M1 of the system is $`q=0.703\pm 0.012`$.

### 3.2 Photometry and variability

Van den Berg et al. (2002) have photometrically monitored both stars, with 1-m telescopes at ESO (La Silla), Kitt Peak and La Palma, in order to look for light variations that may provide insight into the nature of the systems. We here summarize their main results.

Van den Berg et al. (2002; Figure 5) also find S 1063 to be variable, but as yet no period has been identified in the photometric data. The variability is clearly of a long-term nature. If it is periodic, the 18-day-long ESO light curve establishes that the period is longer than the pseudo-synchronous period of 14.64 days (Hut 1981, eq.42; see also Claret & Giménez 1993). Furthermore, the ESO and Kitt Peak data do not phase up on the orbital period of 18.4 days. Longer intervals of continuous observation are required to search for periodicity at longer periods.

Color variations of S 1063 appear to be wavelength dependent. The $`BV`$ color shows little or no variation as the system brightens. Both $`UV`$ and $`VI`$ become bluer by as much as 0.05 mag with increased system brightness.

In contrast, a periodicity comparable to the orbital period was quickly evident in the S 1113 data. We show in Fig. 4 the $`V`$ light curves derived from the 18-day ESO observations and phased to the radial-velocity orbit solution, which is also shown in Fig. 4. The morphology of the light curve is sinusoidal with a high degree of symmetry. Maximum brightness occurs at orbital phase 0.5, when the stars are aligned perpendicular to our line of sight and the primary is approaching. The system becomes bluer in all bands as it brightens, although again it is most evident in $`VI`$.

A time-series analysis on all of the photometric data indicates a possible period of 2.822$`\pm `$ 0.001 days. This period is very similar to the orbital period of 2.823094$`\pm `$0.000014 days, indicating a close association. There is some indication that the light curve shape is not stable. As is evident in Figure 4, three sets of photometric observations obtained at intervals of 4-6 weeks and 2 years do not repeat to within the formal photometric errors, suggestive that the source of the photometric modulation has undergone small changes on these time scales.

## 4 Cluster membership

We use kinematic evidence to assess the membership of these binaries in the M67 cluster. As noted above, Girard et al. (1989) find the proper-motion membership probability of S 1063 to be 97% and of S 1113 to be 98%. The location of S 1113 at 3 core radii in projection makes the argument for its membership somewhat less secure. However, even weighting by the cluster surface-density distribution Girard et al. (1989) give it a membership probability of 94%, reflecting the fact that S 1113 lies within the half-mass radius of the cluster (Mathieu 1983). At the same time, both binaries have radial velocities within 1 km s-1 of the cluster radial velocity.

The proper motions of Girard et al. (1989) are of very high precision, and thus their ability to distinguish cluster members and field stars is excellent (see their Figure 5). Nonetheless, even the highest probability members have some chance of being field stars. Our radial-velocity survey provides independent information which permits us to evaluate the percentage of field stars among the proper-motion high-probability members. Specifically, we consider the sample of stars with proper-motion membership probability greater than 95%. Within this sample there are 140 stars with multiple radial-velocity measurements obtained in the CfA survey and no indication of velocity variability. Of these 140 stars, only one has a velocity (25.6 km s-1) indicating non-membership. The remainder show a gaussian velocity distribution with a mean velocity of 33.55 km s-1 and a velocity dispersion of 1.06 km s-1. This result remains true for the subset of 39 stars at projected radii greater than 2 core radii, which show a mean velocity of 33.45 km s-1 and a velocity dispersion of 1.07 km s-1. Based on the single field-star interloper, we conclude that the field contamination among those high-probability proper-motion members is of order 1%.

A fraction of these field-star contaminants may have radial velocities similar to that of the cluster. The field radial-velocity dispersion is expected to be much larger than the observed cluster velocity dispersion of 0.5 km s-1. As yet the true field radial-velocity distribution in the direction and magnitude range of the proper-motion study has not been measured. In its place we model the field radial-velocity distribution with a gaussian having a typical galactic disk velocity dispersion of 20 km s-1 (Sparke & Gallagher 2000) and conservatively center the field velocity distribution on the M67 radial velocity.333As a plausibility check, we have examined the distribution of all radial-velocity non-members in our survey. These stars have a mean radial velocity of 22 km s-1 with an r.m.s dispersion of 29 km s-1. Since the survey selection was based on proper-motion membership, this distribution is not independent of the proper-motion distribution and thus may not be appropriate for use in this analysis. Nonetheless these results give plausibility to the adopted radial-velocity distribution. With these assumptions, we find that 4% of field stars have velocities within 1 km s-1 of the cluster velocity, corresponding to two times the observed velocity dispersion of cluster members. Thus we estimate that the percentage of field stars with both high proper-motion membership probability and radial velocities consistent with cluster membership is no more than 0.04%.

There are 246 stars in the proper-motion study with membership probabilities in excess of 95%. Using the binomial theorem and our estimated probability of 0.04%, there is a 0.4% chance that two of these stars are field stars with radial velocities consistent with cluster membership, and a 9% chance that one such star is a field star.

Given that the chances of one or even two field stars being three-dimensional kinematic members of M67 are not entirely negligible, the possibility that S1063 and S1113 are field stars rather than anomalous members must be explored. However, as we shall see in the next section, neither of these systems can be completely explained as field stars. Furthermore, S1063 and S1113 are not the only three-dimensional kinematic members of M67 without evident explanation; for example, S1072 already has been noted in the Introduction. The probability that all of these seemingly anomalous cluster members are simply field stars is indeed negligible. Finally, several of the authors would argue that S1063 and S1113 have been drawn from the much smaller sample of ROSAT-detected stars, and therefore would conclude that the chances that they are field stars are yet smaller than presented here.

## 5 Properties of the stars

S 1063 The most notable property of the orbit of S 1063 may be its ordinariness. If S 1063 fell on the cluster main sequence or subgiant branch, its orbit would merit no further attention. For example, the eccentricity of $`e=0.206`$ is quite typical for main-sequence and subgiant binaries in M 67 with periods $`20`$ days (Mathieu, Latham & Griffin 1990).

If the primary star is more massive than the unseen secondary, then the orbit places an absolute lower limit on the primary and secondary masses of 0.06 $`M_{}`$, which is too low to be useful. Given the anomalous location of S1063 in the M67 CMD, we have no estimate of the primary mass based on stellar evolution theory. However, we will show below that the primary is an evolved star, with the surface gravity of a subgiant. Thus we consider primary masses between 0.7 $`M_{}`$ and 1.3 $`M_{}`$, delimited by the main-sequence lifetime equal to the age of the Galaxy and the upper mass limit on the formation of a subgiant branch. For this range on the primary mass the orbit places a lower limit on the secondary mass between 0.23 $`M_{}`$ and 0.34 $`M_{}`$.

Evidently the primary star dominates the light at 5187 Å. As noted above, the CfA spectra indicate an effective temperature of $`T_{\mathrm{eff}}=5000`$ K and a gravity of $`\mathrm{log}g=3.5`$. New analyses of the higher-signal-to-noise spectra of van den Berg et al. (1999) yield similar measures. Comparison of the V I 6251.83 Å/Fe I 6252.57 Å line ratio with the spectral analyses of Gray (Gray 1989, Gray & Johanson 1991) yield a $`T_{\mathrm{eff}}`$ of 5150 K for a dwarf and 4900 K for a giant. Alternatively, in Fig. 5 we use the spectral diagnostic $`I_s`$ as defined by Malyuto & Schmidt-Kaler (1997) to derive a spectral type of G8–K0. Such a spectral classification indicates effective temperatures similar to the spectral line results (Schmidt-Kaler 1982, Bessell & Brett 1988, Bessell et al. 1998).

Detailed analyses of the optical spectra obtained by van den Berg et al. (1999) indicate that the primary of S 1063 is neither a giant nor a dwarf. We have taken two approaches to determining the surface gravity. First, we have compared the gravity-sensitive Mg b lines with Kurucz spectra. In Fig. 6 we show the spectrum of S 1063 and Kurucz spectra at $`T_{\mathrm{eff}}=5000`$ K over a range of $`\mathrm{log}g`$. Based on the Mg b triplet line shapes, S 1063 is neither a giant ($`\mathrm{log}g=2.53.0`$) nor a dwarf ($`\mathrm{log}g=4.55.0`$). Spectral fitting to the Mg I b lines (with $`T_{\mathrm{eff}}`$ a free parameter) results in $`T_{\mathrm{eff}}=5000`$ K, $`\mathrm{log}g`$ = 3.5, and $`v\mathrm{sin}i`$ = 8 km s-1, corroborating the results obtained from the lower-signal-to-noise CfA spectra. However, the referee has correctly noted several disparities of relative line depths in other wavelength regions between the observed spectrum and the models, for example comparing the Fe I 5187.9 and Ti II/Ca I 5188.8 lines. More detailed spectroscopic analyses and modeling may prove fruitful.

Second, in Fig. 5 we compare the values of $`I_s`$ and $`I_2`$ for S 1063 with those of stars of all luminosity classes. Again, for stars of similar $`I_s`$ S 1063 falls between the dwarfs and giants. Thus all of the spectroscopic evidence indicates that the primary of S1063 is a subgiant with an effective temperature of $`T_{\mathrm{eff}}5000`$ K.

Montgomery et al. (1993) provide UBVI photometry (Table 5) which allows an analysis of the extinction toward S1063 under the assumption that the luminosity of the system derives solely from the primary star. We initially consider the extinction of E(B-V) = 0.032 found for M67 by Nissen et al. (1987). Adopting the extinction curve of Cardelli et al. (1989), the dereddened colors of S1063 are given in Table 5. For comparison we also give in Table 5 the colors of a solar-metallicity 4-Gyr 1.37 $`M_{}`$ subgiant with the same value of (B-V)o, using the Yale isochrones with the GDK color transformation (Yi et al. 2001). The intrinsic colors of S1063 are both 0.15 mag bluer in (U-B) and 0.14 mag redder in (V-I). Furthermore, the effective temperature implied by the (B-V)o = 1.02 mag is 4577 K, significantly cooler than the effective temperature of 5000 K derived from spectroscopy.

These color differences can be resolved if the extinction toward S1063 (or equivalently the intrinsic color (B-V)o of S1063) is taken as a free parameter. For example, we consider the colors of a $``$ 5000 K subgiant, specifically a 4 Gyr 1.35 $`M_{}`$ subgiant as defined by the Yale models and presented in Table 5. Such a star implies a reddening of 0.18 mag, and the consequent intrinsic colors of S1063 are also given in Table 5. The dereddened (U-B)o, (B-V)o, and (V-I)o colors for S1063 agree very well with the 4989 K subgiant model.

We note that given the observed photometric variability of S1063 (Section 3.2), detailed modeling of colors as in Table 5 must be considered with care. S1063 shows little variation in (B-V), and so combined with the spectroscopic derivation of the effective temperature the argument for enhanced extinction toward S1063 remains unchanged. The argument is also not sensitive to the surface gravity of the primary star, with which the colors vary only slightly. On the other hand, the observed (U-B) and (V-I) of S1063 show variations of up to 0.05 mag peak-to-peak, and so there is uncertainty in the use of the single-epoch Mongtomery et al. photometry to specify the nature of the system. (Perhaps surprisingly, the (U-B) measurements of Mongtomery et al., Racine (1971) and Sanders (1989) agree to within 0.01 mag; the (B-V) measurements agree to within 0.02 mag.)

The amount of inferred extinction also depends on the choice of stellar model. The Yale models have also been transformed to observables using the LCD transformations (Yi et al. 2001). In Table 5 we repeat the subgiant analysis with these colors. While the derived reddening of E (B-V) = 0.14 mag is somewhat smaller than with the GDK transformations, it nonetheless remains larger than the extinction toward M67 itself.

Interestingly, correcting for a reddening of E (B-V) = 0.18 mag moves S1063 to a position 0.15 mag in V below the base of the M67 giant branch. Given an observed variation in V of 0.15 mag, S1063 remains a candidate for a normal cluster subgiant if an explanation for the higher extinction can be found. Adopting this reddening, the extinction curve of Cardelli et al. (1989), the bolometric correction of the Yale model, and a true distance modulus to M67 of 9.6, we can derive a radius of 2.4 $`R_{}`$ and a luminosity of 3.3 $`L_{}`$ for the primary star (Table 5). Given a log g = 3.5 this radius implies a mass of 0.7 $`M_{}`$; however we have not determined the gravity to better than a factor 2-3. We stress that all of these values are only meaningful given cluster membership.

Alternatively, the colors of S1063 can be modeled with a background giant star having an effective temperature somewhat in excess of 4900 K. For specificity, we show in Table 5 the intrinsic colors of a 0.8 Gyr 2.3 $`M_{}`$ giant, using the Yale models with the LCD transformation. This model agrees well with the intrinsic colors of S1063 given E (B-V) = 0.14 mag. With $`M_v=1.1`$ mag, the distance of such a giant would be 2.8 kpc, roughly 1.5 kpc above the galactic plane. Placing S1063 at such a distance does not encounter significant difficulties with other observations. For example, the resulting luminosity in the Ca H and K lines is typical for RS CVn’s with giant primaries (van den Berg et al. 1999). Similarly, the X-ray luminosity would become 8 x 1031 ergs/sec (0.1 - 2.4 keV). Examination of Figure 4 in Belloni et al. (1998) shows such a flux to be at the high end of observed X-ray fluxes from similar RS CVn primary stars.

However, the large distance does not explain the required extinction. The extinction maps of Schlegel et al. (1998) indicate a Galactic extinction in the direction of M67 of E(B-V) = 0.03 mag, essentially placing all of the Galactic extinction in front of the cluster. This giant model for the primary also has the modest difficulty of a lower surface gravity (log g = 2.9) compared with the log g = 3.5 derived from the spectral analyses.

In summary, these color analyses suggest a reddening of S1063 of order 0.15 mag more than can be attributed to Galactic extinction in the direction of M67. The system might be interpreted as a cluster subgiant or as a distant ($``$ 3 kpc) halo giant, but neither interpretation naturally explains the enhanced reddening. A conjecture is that the reddening might be due to circumstellar or circumbinary material local to the system. Indeed, the slow but substantial V brightness variation observed in this system might be attributed to variations in local extinction; however such an explanation is confounded by the lack of associated variation in (B-V). We also note that the system shows no evidence of a near-infrared excess. Infrared colors have been obtained from the 2MASS survey and are given in Table 5. Theoretical colors are also available from the Yale models with the LCD transformation, also given in Table 5 for the 1.35 $`M_{}`$ model at 4 Gyr. The dereddened observed colors and theoretical colors agree to within the 4% measurement errors on the colors. Reducing these errors and extending measurements to longer wavelengths would be worthwhile.

Of course, these color analyses presume that the colors of the primary star are not abnormal for a 5000 K star, and that the composite light of the system derives only from the primary photosphere. We encourage more careful spectroscopic analyses of S1063, among other things to precisely compare abundances with cluster members, to better determine the surface gravity, and to search for evidence of extincting matter. We also recommend photometric observations over a broader wavelength range to search for other contributors to the composite light.

Finally, we briefly discuss the rotation of the primary star in the context of an M67 1.3 M subgiant with a system E (B-V) = 0.18 mag. The non-detection of the secondary in our spectra requires that the secondary star have $`V>16`$, roughly two magnitudes fainter than the primary (here taken as equivalent to the composite light). If the secondary is a main-sequence star in the cluster, then based on this magnitude limit and extinction the Yale isochrones give an upper limit on the mass of the secondary of 0.9 $`M_{}`$ and an associated lower limit on the inclination angle of 28 deg from the orbital solution. For aligned rotational and orbital axes and a measured vsin i of 6 km s-1, the true equatorial rotation velocity could be as large as 13 km s-1. For a primary radius of 2.4 $`R_{}`$ the rotational period could be as short as 9.5 days. Such a short rotation period would be consistent with the observed X-ray emission, chromospheric activity, and possible flares. Perhaps most plausibly the pseudo-synchronous rotation period of 14.64 days (cf. Claret & Giminez 1993), requiring an inclination angle of $`i=46^{}`$, is permitted by the primary mass range and the upper limit on the secondary mass. However, we remind the reader that the photometric data show no evidence of periodicity at this rotation period, nor is there strong H$`\alpha `$ emission from the primary (van den Berg et al. 1999).

S 1113 The spectroscopic detection of the secondary permits substantially more information to be derived for this binary. Here we analyze the system first from a distance-independent geometric perspective, and then we take a distance-dependent photometric approach.

For the geometric argument, we presume that the stellar rotational periods are synchronized with the orbital period, and that the rotational and orbital axes are aligned. Both are expected for a circularized orbit, since alignment and synchronization happen on shorter timescales than circularization (Hut 1981), and synchronization is further supported by the observed photometric variation on a period very similar to the orbital period.

Given these two assumptions, the mean density of either star can be written as

$$\rho \frac{\left(M\mathrm{sin}^3i\right)}{\left(Pv\mathrm{sin}i\right)^3}.$$
(1)
The mean density is thus independent of the unknown inclination angle and completely determined by the observables Msin$`{}_{}{}^{3}i`$, P, and vsin i. Using the values for these observables and Teff given in Section 3.1, both stars can be placed on the density-effective temperature diagram, as shown in Figure 7. Also shown in Figure 7 are solar-metallicity Yale isochrones for a set of ages. The 4 Gyr isochrone, comparable to the age of M67, is shown as a dashed line.

The mean density of the secondary indicates a main-sequence star with an age of between 1-10 Gyr. While the uncertainties are too large to precisely determine the secondary age, the inferred secondary mass varies little as a function of age. For a 4 Gyr isochrone the mass of the secondary star is 0.93 $`M_{}`$, and the range on mass implied by the uncertainties extends from 0.89 $`M_{}`$ to 0.97 $`M_{}`$. For a mass of 0.93 $`M_{}`$, the dynamical mass ratio implies a mass for the primary star of 1.33 $`M_{}`$ with a range of 1.27 $`M_{}`$ to 1.38 $`M_{}`$. The inclination of the system becomes 48, and the semi-major axis is 11.0 $`R_{}`$. The primary and secondary radii, presuming synchronous rotation for both stars, are 4.0 $`R_{}`$ and 0.90 $`R_{}`$, respectively.

This derived mass for the primary star is notable in that the mass of a 4 Gyr giant is 1.37 $`M_{}`$. That is, the derived mass of the primary star is very similar to that expected for an evolved M67 member. Similarly, the location of the primary star in Figure 7 is closely consistent with the 4 Gyr isochrone of M67, although the uncertainties in effective temperature permit a wide range of ages for the primary star.

A second independent analysis of the primary and secondary stars can be done from the observed photometry and the spectroscopically determined magnitude difference of $`\mathrm{\Delta }V=1.23`$ mag, assuming no other contributors to the composite light. Adopting again the photometry of Montgomery et al. (1993) for the composite system, the primary and secondary $`V`$ magnitudes are derived immediately from the measured flux ratio (Table 6).

To obtain colors for the primary and secondary stars, we begin by presuming membership in M67 and use the 5 Gyr Yale models (GDK transformation). The derived secondary $`V`$ magnitude, a reddening of E (B-V) = 0.032 mag, and an apparent distance modulus of 9.6 mag (Montgomery et al. 1993) gives a secondary absolute magnitude $`M_v`$ = 5.6, indicating a 4 Gyr 0.90 $`M_{}`$ main-sequence model for the secondary star (Table 6). Notably, this mass is similar to the secondary mass of 0.93 $`M_{}`$ derived from independent geometric arguments. It is also encouraging that the effective temperature of the 0.90 $`M_{}`$ cluster member is 5311 K, within the measurement uncertainties of the spectroscopically determined effective temperature of 5500 K (Section 3.1). 444Turning this argument around, we adopt a 4 Gyr, 5500 K star as a model for the secondary. The implied mass is 0.95 $`M_{}`$ with an absolute V magnitude of Mv = 5.3. Unreddened, the distance modulus to such a star is 10.0 mag; given reddening this is an upper limit. Thus, if the secondary is a main-sequence star, S1113 cannot be located far beyond M67.

If we adopt the model colors for a 0.90 $`M_{}`$ 5 Gyr main-sequence star, then we can determine colors for the primary star as presented in Table 6. These results for the primary and secondary are also plotted as open boxes in Fig. 1. The primary lies somewhat further below the subgiant branch and further to the red than indicated by the composite light.

This color deconvolution leaves an interesting loose end, in that the primary shows a (U-B) blue excess of 0.3 mag relative to either main-sequence or giant colors. A giant model selected to match the dereddened primary (B-V) color is shown in Table 6; the (U-B) color of this model is significantly redder than derived for the primary star. The (U-B) color of a main-sequence star of the same (B-V) is 0.89. Both the U-band measurements of Sanders (1989) and Montgomery et al. (1993) support such an excess, although their measurements differ by 0.07 mag and possibly indicate variability. This excess ultraviolet light may in fact be associated with the primary star, or it may derive from elsewhere in the system.

Since the gravity of the primary star is weakly constrained by the spectroscopy, we use both main-sequence and giant models of the same B-V color as the primary to provide bolometric corrections and luminosities at the distance of M67. The giant model implies a luminosity for the primary of 2.0 $`L_{}`$ and an effective temperature of 4534 K. The main-sequence model provides a luminosity of 2.1 $`L_{}`$ and effective temperature of 4726 K. Both of these temperatures are somewhat lower than the spectroscopically determined effective temperature of 4800 K (Section 3.1).

Unlike the secondary star, the model for the primary star inferred from these photometric arguments is markedly different from the model for the primary developed from the geometrical arguments. The photometric radius for the primary of 2.0 $`R_{}`$ is significantly smaller than the geometric radius of 4.0 $`R_{}`$, or equivalently the photometric luminosity of 2.0 $`L_{}`$ is substantially smaller than the geometric luminosity of 7.7 $`L_{}`$.

This discrepancy between the primary luminosities derived from geometric and photometric arguments can be cast in a distance-independent manner. The geometric argument provides a distance-independent prediction for the ratio of primary to secondary V flux. Continuing to assume synchronization and alignment of axes, the ratio of primary to secondary radii is the same as the ratio of vsin i. Using also the ratio of our spectroscopically determined effective temperatures (Section 3.1), the ratio of bolometric luminosities then becomes $`L_{2,bol}/L_{1,bol}=0.085`$. Using the Yale bolometric corrections, we find that $`L_{2,V}/L_{1,V}=0.11`$. This ratio is significantly smaller than the measured flux ratio at V of 0.352 $`\pm `$ 0.016 (Section 3.1). The sense is consistent with the primary luminosity being a factor 3 less luminous than predicted by the geometric argument, as in the last paragraph.

## 6 Discussion

S 1063 The binary S1063 is remarkable in three ways: low apparent brightness compared to M67 subgiants of similar color, enhanced reddening compared to Galactic extinction along the line of sight, and weak H$`\alpha `$ emission observed at a velocity different from the primary star velocity. In addition, the long timescale light variations are not easily explained if they are aperiodic. If they are periodic on timescales longer than 18 days, then they likely are due to star spots on the primary star and indicate that the primary is not rotating pseudosynchronously.

Given both its close kinematic association with M67 and projection upon the core of the cluster, we have explored explanations for S1063 in the context of cluster membership. In this context, we have previously inferred a radius of 2.4 $`R_{}`$ for the primary star. As such, the primary star does not approach its critical surface at periastron passage. Thus the radius of the star is not confined by the presence of the secondary, nor is there motivation for present mass transfer. Similarly, the eccentricity of the orbit does not suggest either a large evolved star or mass transfer in the past, both of which would likely have circularized the orbit.

We have explored the possibility that rather than being an underluminous subgiant the primary star is instead an overluminous main-sequence star. For example, a recent merger of two main-sequence stars would deposit kinetic energy into the merging stars. Here we consider two types of mergers: coalescence and collision. In the coalescence scenario, the S 1063 system was originally a triple system in which the present secondary was the tertiary star. If we assume that the orbital period that we observe now was equal to the period of the original tertiary star, we can use stability arguments to place an upper limit on the period of the original inner binary star. Using the coplanar formulation of Mardling & Aarseth (1998) for the stability of a triple system, and adopting a mass ratio of $`q=0.5`$ (outer star to inner binary), we find that an original inner binary would have had an orbital period of less than 2.22 days. Such a period is physically permitted, and so that S 1063 was formerly a triple system is possible. Coalescence via multiple mechanisms is plausible on time scales of a few Gyr (Stȩpień 1995). Four contact binaries are known in M67 (see light curves in van den Berg et al. (2002) and Stassun et al. (2002)), and the old (6 Gyr) cluster NGC 188 has many contact binaries (Baliunas & Guinan 1985, Rucinski 1998).

A collision scenario has been suggested by Leonard & Linnell (1992) as a possible mechanism for the formation of some blue stragglers in M 67. In the low-density environment of an open cluster, collisions of single stars are improbable. However, binary-binary encounters have substantially higher probabilities, and in the course of such encounters the probability of a stellar collision is not negligible. In this scenario, the present S 1063 would consist of three of the stars involved in the binary-binary encounter, two having merged into the primary and a third being the present secondary. An eccentric orbit is a natural result of such a resonant encounter. A prediction of either of these merger scenarios is that the resulting star will rotate rapidly. In the coalescence scenario rapid rotation is expected as a result of synchronous rotation prior to the merger. With respect to the collision scenario, the simulations of Sills et al. (2001) show that a collision product is born with a high rotation rate. We argued in Sect. 5 that for aligned rotation and orbital axes, the true equatorial rotation velocity could be as large as 13 km s-1, and for a primary radius of 2.4 $`R_{}`$ the rotational period could be as short as 10 days. However such a rotation period is still long compared to those expected from mergers or collisions. In the collision scenario, of course, the primary rotation axis and the orbital axis need not be aligned, particularly since tidal circularization has not come to completion. As such, the surface rotation velocity of the primary could be significantly higher. A high inclination angle could also possibly explain a lack of periodic photometric modulation.

Most problematic, these impulsive origins must face the challenge that – even if the merger product would briefly have the properties of the S 1063 primary star – the thermal time scale of the primary is very short compared to the cluster lifetime. If, for the sake of example, we adopt a primary mass of 0.7 $`M_{}`$, the present thermal time scale $`E_{\mathrm{pot}}`$/$`L`$ is only 3.4 Myr. A merger product would be expected to readjust to its new mass on such time scales (and perhaps become a present or future blue straggler depending on the combined mass). Given such short time scales, the probability of observing the primary prior to its new equilibrium state is very low. On the other hand, should a star be found shortly after a merger, the presence of residual circumstellar material resulting from the merger might also be expected, and represents a possible explanation for the enhanced reddening of this system. Lastly, if formed recently such a system may not yet be pseudo-synchronized.

We stress that the nature of the secondary is unknown; we know only that it is substantially fainter than the primary at $`V`$. Thus the secondary is permitted to be a white dwarf. An isolated white dwarf of the cluster age would have a mass of order 0.6–0.7 $`M_{}`$ (Wood 1992, Richer et al. 1998); a white dwarf deriving from a prior mass transfer scenario could be much different. A hot source of radiation might help explain the weak H$`\alpha `$ emission observed at a velocity different from the primary star velocity. How the secondary could become a white dwarf without circularizing the orbit during its prior post–main-sequence evolution would need to be explained; again, a dynamical exchange might play a role.

S 1113 In Section 5 we derived the properties of the stars in S1113 from two independent lines of reasoning, geometric and photometric. The two arguments produced quite different conclusions regarding the luminosity ratio of the stars, or similarly regarding the nature of the primary star. To bring these two lines of reasoning into agreement, we must give up a basic premise of at least one of the arguments. Here we explore two possibilities, and then discuss the measurement uncertainties.

a) Both the geometric and photometric arguments converge on a 0.9 $`M_{}`$ main-sequence secondary. This suggests that the resolution of the contradiction in luminosity ratios might be found in the primary star. Under the assumption of membership in M67 the geometric argument must be found wrong, for its implied 4800 K primary star of radius 4.0 $`R_{}`$ would be 1.2 mag more luminous in V than observed. The problem cannot be solved with extinction, since the reddening vector for the S1113 primary does not pass through an M67 4.0 $`R_{}`$ giant.

Examining the premises of the geometric argument, we first note that relaxing the assumption of alignment of orbital and rotational axes cannot in itself resolve the problem. Even if the inclination angle of the primary rotation axis is taken to be 90, its radius is reduced only to 3.0 $`R_{}`$. (Recall that the primary radius derived from the photometric argument is 2.0 $`R_{}`$.)

Supersynchronous rotation of the primary leads to more success. Specifically, a rotation period of 1.4 days - or supersynchronous rotation by a factor of 2 - would bring the primary radii derived from both the geometric and photometric arguments into agreement on a primary radius of 2.0 $`R_{}`$.

Such supersynchronism is not expected theoretically given the close proximity of the stars, the primary’s large radius, and the surface convection zones on both stars. Unless continuously driven, the duration of a supersynchronous state would be short. Furthermore, supersynchronous rotation of the primary would also require that the observed variability of the composite light with the orbital period must derive from elsewhere in the system. A cool spot on the secondary seems an unlikely origin. Given a secondary/primary V luminosity ratio of 0.35, a 38% light modulation from the secondary would be required to reproduce the observed 10% variation in the composite light at V.

We note that making the secondary subsynchronous would imply a larger, more luminous secondary which would also resolve the flux ratio discrepancy. However the consequent change in mean density moves the secondary off the main sequence into a domain not occupied by standard stellar evolution (Figure 7). While perhaps this is the case, this path permits essentially unconstrained modeling of the binary.

Finally, the unusual CMD position of S1113 is not easily solved by enhanced extinction of a cluster member. Large extinctions (in excess of 1 mag) along a standard reddening vector would imply a primary star near the top of the M67 main sequence, with much higher effective temperatures than found spectroscopically.

b) Alternatively, we consider arbitrarily discounting the spectroscopically determined flux ratio at V of 0.32 (Sec. 3.1) and adopt the geometric model for the system. This model indicates a secondary/primary flux ratio at V of 0.11. Corresponding decomposition of the composite V magnitude leads to a primary V = 13.88 and a secondary V = 16.31. The geometric model indicates a 0.9 $`M_{}`$ main-sequence secondary, and so without considering additional reddening the secondary V magnitude places the binary 1 mag beyond M67. A similar conclusion is reached by comparing the luminosity of a 4800 K, 4.0 $`R_{}`$ primary star to the composite light. At this distance the 1.3 $`M_{}`$ primary becomes very similar to a star at the base of the M67 giant branch, and as such can be explained by standard stellar evolution theory (albeit not as a cluster member).

The larger primary radius of 4.0 $`R_{}`$ implied by the geometric model represents a large fraction of the primary Roche radius. The Roche radii are 4.5 $`R_{}`$ and 3.8 $`R_{}`$ for the primary and secondary, respectively. An equipotential surface about the primary whose volume equals that of a 4.0 $`R_{}`$ sphere extends 74% of the way to the $`L_1`$ point. Thus, within the geometric model, the role of mass transfer in the evolution of S1113 merits consideration. Evidently the 0.9 $`R_{}`$ main-sequence secondary lies well within its Roche radius.

Interestingly, both the phasing of the periodic photometric variation and the ultraviolet excess are suggestive of a hot spot near the secondary star powered by an accretion stream. Such spots are formed on the following side of the secondary (e.g., Flannery 1975), in qualitative agreement with the phasing of the light and velocity curves shown in Figure 4. The 10% system brightness variation at V requires a mass accretion rate of 2.5 x 10-8 $`M_{}`$/yr, presuming that the energy is released by a 10,000 K spot on the surface of the secondary star.

We note that at present no spectral signatures of such a hot spot have been seen. Van den Berg et al. (1999) did find broad H$`\alpha `$ from S1113, but this emission was kinematically associated with the primary star. They argued that this emission could be consistent with the higher chromospheric activity driven by the rapid rotation of the star, in analogy to V711 Tau. In this context, we note that the photometric variation can also be well modeled by a large cool spot on the primary star. For example, we have reproduced the light variations with the addition of a cool spot in the upper hemisphere of the primary and located 90 from the major axis of the orbit. Specifically the spot properties are: latitude 40, longitude 270, radius 27, temperature 0.82 of the stellar surface temperature. The location of the spot at 270 longitude is motivated by the observations rather than an independent physical argument.

The large extent to which the primary star fills its Roche lobe in this model implies that its shape is significantly asymmetric. Using the Wilson-Devinney formalism we have investigated the expected magnitude of photometric variations due to the anticipated ellipsoidal shape of the primary. Adopting the spot parameters above, a temperature of 4800 K for the primary star, a limb darkening coefficient (linear law) of 0.6 for each star, a gravity-brightening coefficient of 0.32, and a reflection coefficient of 0.5, we find the peak-to-peak ellipsoidal variations to be 0.06 mag in $`V`$. Most importantly, these ellipsoidal variations produce a double-peaked light curve over the orbital period, in marked contrast to the observed single-peaked light curve. It is possible that these small ellipsoidal variations have gone undetected in contrast to the larger photometric variations of the system. For example, in Fig. 8 we show the light curve combining both the cool spot described above and the ellipsoidal variations. Evidently the ellipsoidal variations are lost to inspection in this synthetic light curve, which in fact looks very similar to the observed light curve.

To summarize, the geometric line of reasoning suggests that S1113 is a field binary located behind M67. In this scenario it may be a detached RS CVn whose rapid rotation is producing the several emission diagnostics of enhanced chromospheric activity as well as a large spot. Alternatively, mass transfer may be underway, producing a hot spot in the vicinity of the secondary star. However, to adopt this interpretation of the system, both the kinematic association with M67 and the agreement of the primary mass with the M67 turnoff mass must be taken as merely coincidental. In addition, the measured V flux ratio must be ignored. As we discuss below, we find this last requirement to be a serious counterargument.

c) Measurement uncertainties We have explored whether the discrepancy between the conclusions of the geometric and photometric analyses are the result of our measurement uncertainties. On the photometric side, broadband photometric measurements of S 1113 are several and corroborative. The effective temperature for the primary derived from photometric colors is similar to the effective temperature derived from our high-resolution spectra, indicating that both determinations of the effective temperature are reasonable. Finally, we find it unlikely that the bolometric correction for the primary star overestimates the luminosity by a factor 3, particularly given the observed excess flux in the U-band. However, given that only UBV photometry is available, photometry over a wider range of wavelength is much needed.

The geometric luminosity ratio can be brought into accord with the observed V ratio if the effective temperatures and projected rotation velocities of both stars are all adjusted by 1.5 $`\sigma `$ in the appropriate senses. While ad hoc and improbable, the required adjustments are small enough relative to their uncertainties that further precise measurements of these quantities are in order.

The measured flux ratio at V is pivotal in this discussion. External checks of this quantity on other eclipsing binary systems, using similar spectroscopic material and with the same techniques used here, have not shown significant systematic errors. For example, a comparison with independent determinations available from the light-curve solutions of double-lined eclipsing binaries typically agree within 10% or better with our spectroscopic determinations (see, e.g., Lacy et al. 1997; 2000). Thus we have no reason to believe that our flux-ratio measurement for S1113 would be in error by as much as a factor 3. Nonetheless, given the importance of this measurement we encourage additional observational study.

In closing, we note that Hurley et al. (2001) have suggested that S1113 is subluminous due to an evolutionary response to mass transfer. Albrow et al. (2001) identify stars below the subgiant branch of the globular cluster 47 Tuc which they note are similar to S1063 and S1113. They suggest that such systems may result from mass transfer, with the subsequent contraction of the primary Roche radius deflating the primary star and making it less luminous. Our analyses indicate that, considered as a cluster member (i.e., the photometric model), both stars in S1113 presently fall well within their Roche radii. Put another way, if the primary of S1113 fills its Roche lobe (e.g., the geometric model) and has an effective temperature of 4800K, the observed brightness places it behind M67 (and as noted above this conclusion of non-membership cannot be removed via extinction arguments). Thus at least presently S1113 would not seem to be a candidate for a Roche-filling cluster binary.

## 7 Summary

The stars S 1063 and S 1113 in the field of the old open cluster M67 have attracted attention for their location in the cluster color-magnitude diagram roughly 1 magnitude below the cluster subgiant branch (Fig. 1). Comprehensive photometric (optical and X-ray), spectroscopic, and astrometric studies have shown them to be RS CVn-like systems with observed three-dimensional motions that are the same as that of the cluster to high precision. Their location in the color-magnitude diagram has led us to describe these stars as candidate “sub-subgiant” cluster members.

Specifically, S 1063 is a single-lined spectroscopic binary with a period of 18.4 days and an eccentric orbit ($`e=0.2`$). The primary effective temperature is 5000 K. Spectroscopic surface gravity measurements indicate a subgiant with a surface gravity log g $``$ 3.5. UBVIJHK colors indicate a reddening E(B-V) = 0.14 - 0.18, higher than the cluster reddening of E(B-V)=0.032. As the entire Galactic reddening in the direction of M67 is only 0.03 mag, the enhanced reddening of S1063 is likely not interstellar. There is no evidence in the near-infrared colors for warm circumstellar dust in this system. The strong X-ray emission suggests rapid rotation of the primary, although the measured $`v\mathrm{sin}i`$ is only 6 km s-1. Pseudosynchronous rotation at a period of 14.64 days is permitted by existing observations, but there is no evidence for photometric variability at this period. Photometric variability is seen on longer timescales, with no evident periodicity for periods less than 18 days.

S 1113 is a double-lined spectroscopic binary with a period of 2.8 days, a circular orbit, and a mass ratio of 0.7. Effective temperatures of 4800 K and 5500 K are found for the primary and secondary stars, respectively. The photometric properties of S1113 are well modeled by a binary with a 0.9 $`M_{}`$ main-sequence secondary star and 1.3 $`M_{}`$ a subgiant primary star (2.0 $`R_{}`$) at a distance similar to M67. However, this model does not easily explain the large projected rotation velocity (53 kms-1) of the primary star. Supersynchronous rotation of the primary is required, in contrast to the short synchronization timescales, the circular orbit, and the observed variability of the composite light with the orbital period.

Alternatively, mean densities for both stars are derived geometrically assuming a tidally relaxed system. These densities indicate a 0.9 $`M_{}`$ main-sequence secondary star and a 1.3 $`M_{}`$, 4.0 $`R_{}`$ primary star that very nearly fills its Roche lobe. In this model S1113 would be an RS CVn (possibly transferring mass) located behind the cluster. However, the V flux ratio derived from this model differs by a factor 3 from the spectroscopically measured flux ratio. In addition, the observed light curve shows no evidence for asymmetry of the nearly-Roche-filling primary star. Thus neither the photometric nor the geometric analyses are able to provide a fully consistent model for the S1113 system.

The issue of cluster membership is pivotal for understanding both S1063 and S1113. Given their low apparent luminosity, the interpretation of both binaries as background field RS CVn’s is a natural solution. This interpretation runs counter to the very close kinematic three-dimensional association of both binaries to the cluster. Still, our statistical analyses indicate that the probability of one of the two being a background binary is not entirely negligible, although both being background stars is improbable. Equally important, we have not been able to create model background binaries which reproduce all of the observed properties of either S1063 or S1113.

Star clusters have long been recognized as dynamically active environments, particularly in the form of binary encounters with both single stars and other binaries. It would be no surprise if we were to find products of such encounters that run counter to single-star evolutionary theory. Nor would it be a surprise if the products of such encounters were binary stars. S1063 and S1113 may be recent products of such encounters, and as such may define alternative evolutionary tracks in the cluster environment. Regrettably we have not yet been able to explain their low luminosities, although we conjecture that their evolutionary histories include mass-transfer episodes, mergers, and/or dynamical stellar exchanges.

Albrow et al. (2001) have noted stars in the globular cluster 47 Tuc which fall in a similar domain of the color-magnitude diagram as S1063 and S1113. Their finding of more examples of stars like S1063 and S1113 emphasizes the need for searches for sub-subgiants in other clusters.

Finally, the primaries of both S1063 and S1113 are candidates for future mass transfer. It is intriguing to consider these binaries as progenitors of the blue stragglers, and particularly of the short-period binary blue straggler F190 in M67.

We thank R. Davis, E. Horine, A. Milone and J. Peters for acquiring many of the CfA spectra. We are grateful to J. Orosz for assistance in computing the ellipsoidal variation light curves. We thank Pierre Demarque and Sukyoung Yi for providing their stellar evolution models. We also appreciate a very helpful review by the referee Roger Griffin. RDM and KS are supported by National Science Foundation research grant AST 9731302. MvdB was supported by the Netherlands Organization for Scientific Research.


format: text

text:
# Lack of self-averaging in neutral evolution of proteins

(March 21, 2002)

## Abstract

We simulate neutral evolution of proteins imposing conservation of the thermodynamic stability of the native state in the framework of an effective model of folding thermodynamics. This procedure generates evolutionary trajectories in sequence space which share two universal features for all of the examined proteins. First, the number of neutral mutations fluctuates broadly from one sequence to another, leading to a non-Poissonian substitution process. Second, the number of neutral mutations displays strong correlations along the trajectory, thus causing the breakdown of self-averaging of the resulting evolutionary substitution process.

Introduction. – The neutral theory of molecular evolution is the simplest and most elegant theory of protein evolution. It was formulated in the late 1960’s Kimura:1968+King/Jukes:1969  to explain the high substitution rate of amino acids observed in proteins of many vertebrates and the large intra-specific genetic variations between most species. The theory assumes that most of the amino acid substitutions occurring in an evolving population do not bring any selective advantage but are selectively neutral, maintaining the biological activity of the protein at the original level. Such neutral mutations are assumed to occur at random at a rate $`\mu x`$, where $`\mu `$ represents the genomic mutation rate and $`x`$ is the fraction of mutations which happen to be neutral. The theory predicts that (i) the rate of substitutions (mutations which become fixed in the population) equals the neutral mutation rate $`\mu x`$ and depends only on the protein considered, independent of the size of the population and its ecology, and that (ii) the number of amino acid substitutions taking place in a time $`t`$ follows a Poisson distribution with mean value $`\mu xt`$, thus giving an explanation to the ‘molecular clock’ observed in the early 1960’s Zuckerkandl/Pauling:1962 . Later studies in the 1970’s and 1980’s have revealed, however, that the variance of the substitution process is larger than its mean value Ohta/Kimura:1971+Langley/Fitch:1973 ; Gillespie:1989 , pointing to an underlying non-Poissonian process. Since then, different alternatives have been suggested to explain this feature. Some authors have extended the neutral theory by including into it slightly deleterious mutations Ohta:1976 , others have rejected the neutral theory completely and have suggested that most mutations are fixed by positive selection Gillespie:1991 . An interesting proposal within the realm of neutral theory is its modification in terms of a fluctuating neutral space model Takahata:1987 , which can account for the non-Poissonian statistics.

Progress in the understanding of the folding and thermodynamics of biomolecules has opened the way to assess the thermodynamical stability of biomolecules involved in evolution through computational methods, thus providing powerful tools to complement the traditional population genetic approach. This structural approach has been introduced in the study of neutral networks (i.e., the set of sequences connected by structure conserving mutations) of RNA secondary structures Schuster/etal:1994-98 , and has found fruitful applications in the study of protein evolution OTHERSTUDIES . Based on this structural approach, we have studied a model of neutral evolution denoted structurally constrained neutral model (SCN model), in which conservation of the thermodynamic stability of the native structure is imposed Bastolla/etal:2002 . The SCN model yields evolutionary trajectories consisting of sequences connected through neutral mutations. In this Letter, we show that the fraction of neutral mutations obtained along the trajectories fluctuates strongly, and consequently the SCN model produces a non-Poissonian substitution process, consistent with a fluctuating neutral network scenario Takahata:1987  and with the statistics of protein evolution Gillespie:1989 . Moreover, we find that the number of neutral neighbors displays strong auto-correlations along the trajectories, which deeply influence the statistics of protein evolution as we discuss below.

We have studied seven protein folds: myoglobin (Protein Database (PDB) code 1a6g), cytochrome c (PDB code 451c), lysozyme (PDB code 3lzt), ribonuclease (PDB code 7rsa), rubredoxin (both from a mesophilic and a thermophilic species, PDB codes 1iro and 1brf\_A), ubiquitin (PDB code 1u9a\_A), and the TIM barrel (PDB code 7tim\_ABastolla/etal:2002 . Although these proteins cover a broad spectrum of different biological activities, according to the SCN model their neutral evolution occurs in a rather similar way, suggesting that our model captures ‘universal’ features of protein evolution. In what follows, we concentrate on myoglobin to illustrate this common scenario.

The SCN model Bastolla/etal:2002 ; Bastolla/etal:1999+2000 . – The simulated trajectories are started from a given “wild-type” sequence $`𝐀^{}`$ of $`N`$ amino acids (residues) that folds onto the native structure $`𝐂^{}`$, where both $`𝐀^{}`$ and $`𝐂^{}`$ are taken from the PDB. The target structure $`𝐂^{}`$ is kept fixed throughout the simulation. Following the neutral theory, a given mutated sequence $`𝐀^{}`$ is considered to be either neutral, if it folds onto $`𝐂^{}`$ preserving thermodynamic stability, or lethal otherwise. Neither advantageous nor slightly deleterious mutations are considered. At each iteration, we generate all possible sequences $`\{𝐀^{}\}`$ obtained through point mutations of the current sequence $`𝐀`$ Note , and determine the fraction belonging to the neutral network. One of these neutral sequences is chosen at random and becomes the new current sequence. The whole process is iterated, resulting in an evolutionary trajectory $`\{𝐀_0=𝐀^{},𝐀_1,𝐀_2,\mathrm{}\}`$.

To assess the conservation of the thermodynamic stability of a test sequence in the target structure, we rely on two well-established empirical parameters, the normalized energy gap $`\alpha `$ Bastolla/etal:1999+2000 , measuring the minimal value of the energy gap between an alternative conformation and the target one divided by their structural distance, and the $`Z`$-score Bowie/etal:1991+Goldstein/etal:1992 , measuring the difference between the native energy and the average energy of alternative compact conformations in units of the standard deviation of the energy. A positive and large value of the $`\alpha `$-parameter ensures both that the target conformation has lowest energy and that the energy landscape is well correlated, the latter in the sense that conformations very different from the native have much higher energy. Moreover, the native energy allows a rough estimate of the folding free energy. Alternative structures are generated by aligning in all possible ways the test sequence with all non-redundant protein structures in the PDB. Conformational energies are computed through the effective energy function described in Ref. Bastolla/etal:2000+Bastolla/etal:2001 . The resulting estimate of the thermodynamic stability provides a realistic genotype to phenotype mapping, although the energy function used is approximate and alternative structures are non-exhaustively sampled.

Neutral connectivity. – In the SCN model, the quantity of interest is the fraction of neutral neighbors $`x_i=X_i/X_{\mathrm{max}}(0,1]`$ associated to each sequence $`𝐀_i`$, as this quantity determines the evolutionary substitution process (see below). Here, $`X_i`$ is the number of neutral neighbors of sequence $`𝐀_i`$, and $`X_{\mathrm{max}}`$ is its maximal possible value Note . Thus, an evolutionary trajectory is an (ideally infinite) series $`𝐱=\{x_0,x_1,x_2,\mathrm{}\}`$. For each position $`k`$, we define analogously evolutionary trajectories $`𝐱^{(k)}`$, where $`x_i^{(k)}`$ counts the fraction of the $`18`$ possible changes at position $`k`$ which are neutral. The probability distribution $`P(x)`$, where $`P(x)\mathrm{d}x`$ is the probability to find a fraction of neutral neighbors between $`x`$ and $`x+\mathrm{d}x`$ in the trajectory $`𝐱`$, obtained for myoglobin is shown in Fig. 1(a). The fraction of neutral neighbors is broadly distributed, reminiscent of a fluctuating neutral network Takahata:1987 .

To get a deeper insight into the mechanism of neutral evolution, it is helpful to investigate the auto-correlation of the trajectory $`𝐱`$. For this purpose, we study the ‘profile’ $`p(n)`$ of walks of $`n`$ steps, defined as $`p(n)=_{i=1}^n\left[x_i\overline{x}\right]`$, where $`\overline{x}=_0^1xP(x)dx`$ is the average over all $`x_i`$ in the trajectory (cf. inset in Fig. 1(b)). The ‘roughness’ of the profile tells us about the presence or absence of correlations along the trajectory $`𝐱`$. They can be determined quantitatively by calculating the fluctuations of $`p(n)`$ within a window of width $`\mathrm{}`$ as $`F(\mathrm{})=<[p(n)p(n+\mathrm{})]^2>^{1/2}`$ (see Ref. Peng/etal:1994+Koscielny-Bunde/etal:1998 ). If $`𝐱`$ is a series of independent variables, $`F(\mathrm{})`$ scales as $`\mathrm{}^{1/2}`$. Instead, if the series $`𝐱`$ is correlated (anti-correlated) over a length $`\mathrm{}_0`$, the scaling is modified into $`F(\mathrm{})\mathrm{}^\nu `$ with $`\nu >1/2`$ ($`\nu <1/2`$) on the same length scale. The information provided by $`F(\mathrm{})`$ is equivalent to the one given by the auto-correlation function, but the former has the advantage that different regimes and intermediate crossovers are more easily detected.

The plot of the fluctuations of an evolutionary trajectory $`𝐱`$ obtained for myoglobin is shown in Fig. 1(b), together with the analogous analysis of the fraction of neutral neighbors for residues $`25`$ and $`115`$, respectively the residues with the largest and the smallest mean fraction of neutral neighbors. The fraction of neutral neighbors for the whole protein is strongly auto-correlated, with an apparent exponent $`\nu 0.82`$ for about $`20`$ evolutionary steps, after which the random walk exponent $`\nu =1/2`$ is attained. Single residues have a somewhat larger apparent exponent between $`\nu 0.85`$ (residue $`25`$) and $`\nu 0.9`$ (residue $`115`$) for about $`10^2`$ to $`10^3`$ evolutionary steps, which is roughly the product of number of steps it takes the correlations to vanish for the whole protein times the protein length. These auto-correlations can be qualitatively explained by the fact that more stable sequences have a larger number of neutral neighbors (they are more tolerant to mutations), and stability itself is auto-correlated along an evolutionary trajectory as long as the number of mutations is small with respect to sequence length. These correlations, although being short range, have an important influence on the evolutionary substitution process, as shown below.

Modeling the substitution process. – In order to construct the substitution process, we have to obtain the number of mutations accepted within a time interval $`t`$. The substitution process consists of three steps: (i) generation of an evolutionary trajectory using the SCN model; (ii) determination of the number $`k`$ of mutation events taking place within the time $`t`$, which is assumed to be a Poissonian variable of average $`\mu t`$, so that the probability of $`k`$ mutations within time $`t`$ is $`P_t(k)=\mathrm{exp}(\mu t)(\mu t)^k/k!`$; (iii) determination of the number $`n`$ of accepted mutation events out of $`k`$, where the corresponding conditional probability $`P_{\mathrm{acc}}(nk)=\left[_{i=1}^nx_i\right]_{\{m_j\}}_{j=1}^{n+1}(1x_j)^{m_j}`$. Here, the $`\{m_j\}`$ are all integer numbers between zero and $`kn`$ satisfying $`_{j=1}^{n+1}m_j=kn`$. In other words, the probability that a mutation is accepted is $`x_0=x(𝐀_0)`$ as long as the protein sequence is $`𝐀_0`$, $`x_1=x(𝐀_1)`$ as long as the sequence is $`𝐀_1`$, and so on. If all $`x_i=x`$ are equal, $`P_{\mathrm{acc}}(nk)`$ reduces to a binomial distribution $`P_{\mathrm{acc}}(nk)=\left(\genfrac{}{}{0pt}{}{k}{n}\right)x^n(1x)^{kn}`$. The probability $`\mathrm{\Pi }_t(n)`$ that $`n`$ mutations are accepted within time $`t`$ is the weighted sum over $`k`$ of the acceptance probability, $`\mathrm{\Pi }_t(n)=_{kn}P_t(k)P_{\mathrm{acc}}(nk)`$, which in the case of equal $`x_i`$ reduces to a simple Poisson distribution $`\mathrm{\Pi }_t(n)=\mathrm{exp}(\mu tx)(\mu tx)^n/n!`$ with average value $`\mu xt`$ and substitution rate $`\mu x`$, as in the original neutral model.

Statistics of the substitution process. – In our analysis, two kinds of averages must be distinguished. We indicate by angular brackets $`<>`$ the average over the mutation and acceptance process for a given realization of the evolutionary trajectory, and by an overline $`\overline{}`$ the average over evolutionary trajectories. We determine the mean number of substitutions $`\overline{<S_t>}`$ within a time $`t`$, and show this quantity in Fig. 2 together with the normalized mutation variance $`R_\mu =(\overline{<S_t^2>}\overline{<S_t>^2})/\overline{<S_t>}`$, the normalized trajectory variance $`R_x=(\overline{<S_t>^2}\overline{<S_t>}^2)/\overline{<S_t>}`$, and the normalized total variance $`R=R_\mu +R_x`$. The latter quantity is also known as dispersion index. Notice that if all $`x_i=x`$ are the same, one gets $`R_\mu =1`$ as for all Poissonian processes, and the normalized trajectory variance $`R_x=0`$, since $`x`$ is constant in all sequence space.

The results of the substitution process based on the evolutionary trajectories are shown in Fig. 2(a). We also show in Fig. 2(b) results based on annealed trajectories, obtained by extracting at random the values of $`x_i`$ at each substitution event according to the observed distribution $`P(x)`$. In this case, the different $`x_i`$ along the trajectories are independent variables. Note that the annealed trajectories ‘interpolate’ between the Poissonian case (all $`x_i`$ are equal) and the correlated trajectories obtained through the SCN model: (i) The comparison between the annealed trajectories and the simple Poissonian case allows us to identify the effect of the broad distribution of the fraction of neutral neighbors, whereas (ii) the comparison between the actual and the annealed trajectories allows us to identify the effect of correlations.

In the annealed case, the time spans $`\tau `$ between subsequent substitutions are independent variables distributed with the density $`\overline{D(\tau )}=_0^1P(x)(\mu x)^1\mathrm{exp}(\mu x\tau )dx`$, whose average value is $`\overline{\tau }=_0^1P(x)(\mu x)^1dx`$. Thus, the average substitution rate $`\overline{<S_t>}/t`$ is not constant in time as in the Poissonian case. Initially, $`S_t`$ is a Poissonian variable with average rate $`\mu \overline{x}`$. At large time, however, the rate converges to the smaller value $`\overline{<S_t>}/t1/\overline{\tau }`$, since the process spends more and more time in sequences with small $`x`$. Hence, the standard deviation of $`\tau `$ is slightly larger than its average value, and the normalized mutation variance $`R_\mu `$ is larger than the Poissonian value $`R_\mu =1`$, although the actual difference is small (see Fig. 2(b)). The normalized trajectory variance $`R_x`$ is very small, of the order of the ratio between variance and average value of the $`x_i`$.

Using the actual evolutionary trajectories (see Fig. 2(a)), we note that the presence of correlations has only a weak effect both on the average number of substitutions $`\overline{<S_t>}`$ and on the normalized mutation variance $`R_\mu `$. However, the normalized trajectory variance increases considerably in response to the correlations, as $`R_x1`$ for $`\mu t=200`$. It even grows with time, although more and more sequences are used to compute the mutational averages and one could naively expect that such averages approach typical values. Hence, the large fluctuations between different trajectories caused by the strong auto-correlations result in the breakdown of self-averaging in the substitution process, in the sense that even averaging over an arbitrary long trajectory does not give values representative of typical trajectories. (We note, however, that the variable $`<S_t>/t`$ is still self-averaging as its variance vanishes in the long $`t`$ limit, but the normalized variable $`<S_t>/\sqrt{t}`$ has a variance which increases with time.) The normalized total variance $`R`$ becomes larger than $`2`$ already for $`\mu t100`$, in agreement with empirical estimates, varying in the range $`1<R5`$ for most of the proteins studied Ohta/Kimura:1971+Langley/Fitch:1973 ; Gillespie:1989 . Hence, these large dispersion indices may be to a large extent due to the correlations present in the evolutionary process.

Conclusions. – We have shown that the evolutionary trajectories in sequence space generated by the SCN model are characterized both by a broad distribution of neutral connectivities and by strong correlations, which deeply influence the statistics of neutral evolution and the substitution process. For example, fluctuations of the evolutionary rate from one branch of the evolutionary tree to the other can obscure lineage effects, i.e. variations of the substitution rate among different taxonomic groups. One such effect is the generation time effect: Since the natural time unit for measuring substitution events is the generation time (at which reproduction takes place), the longer the latter the slower the substitution rate is expected to be. This has been verified by comparing for instance substitution rates in rodents and primates Britten:1986+Li/etal:1987+Grauer/Li:2000 . However, the effect is significantly larger for synonymous substitutions (DNA changes which still code for the same amino acid) than for non-synonymous ones. Non-synonymous substitutions are superimposed with the large and correlated fluctuations in the substitution rate that we just described, while synonymous ones are not. Thus the statistics of neutral substitutions could explain this quantitative difference. Additionally, a better understanding of the mechanisms of neutral evolution will help to single out the more interesting cases of positive selection as, for instance, changes in the protein function and responses to variations of the environment. The best current bioinformatics method to identify such cases of positive selection McDonald/Kreitman:1991 , recently used to study the evolution of Drosophila genes Smith/Eyre-Walker:2002+Fay/etal:2002 , assumes a neutral substitution process with constant fraction of neutral neighbors $`x`$. The broad distribution and the correlations that we observe can mimick the presence of positive selection, and they should be taken into account to improve the performance of such methods.


text:
# Inelastic relaxation and noise temperature in S/N/S junctions

11institutetext: Département de Recherche Fondamentale sur la Matière Condensée / SPSMS / LCP, CEA-Grenoble,
17 rue des Martyrs, 38054 Grenoble cedex 09, FRANCE

## Abstract

We studied electronic relaxation in long diffusive superconductor / normal metal / superconductor (S/N/S) junctions by means of current noise and transport measurements down to very low temperature ($`100mK`$). Samples with normal metal lengths of $`4,10`$ and $`60\mu m`$ have been investigated. In all samples the shot noise increases very rapidly with the voltage. This is interpreted in terms of enhanced heating of the electron gas confined between the two S/N interfaces. Experimental results are analyzed quantitatively taking into account electron-phonon interaction and heat transfer through the S/N interfaces. Transport measurements reveal that in all samples the two S/N interfaces are connected incoherently, as shown by the reentrance of the resistance at low temperature. The complementarity of noise and transport measurements allows us to show that the energy dependence of the reentrance at low voltage is essentially due to the increasing effective temperature of the quasiparticles in the normal metal.

The profound comprehension of the current transport in metals is a topic of permanent interest LesHouches ; Quirion . With the progress in thin film technology a great number of studies deal with coherence phenomena at low temperature in metallic samples of length $`L`$ shorter than the phase coherence length $`L_\varphi `$ (mesoscopic regime). In this context, a lot of works focus on inelastic processes and for instance, on the apparent saturation of the phase breaking length Mohanty ; Pierre . Current noise measurement is particularly well suited to such investigations because it is sensitive to energy relaxation processes and gives access to the involved inelastic scattering length $`L_{in}`$ Naveh ; Blanter .
Current noise in diffusive mesoscopic normal metals connected to two normal reservoirs (N/N/N-case) has been studied by various groups Steinbach ; Schoelkopf ; Henny . Using a Boltzmann-Langevin approach the current noise is given by Nagaev :

$$S_I=\frac{4}{RL}_{L/2}^{L/2}dxdϵf(ϵ,x)[1f(ϵ,x)],$$
(1)
where $`f`$ is the distribution function of the electrons and $`R`$ the resistance of the sample. In the regime $`L<<L_{in}`$, the noise is reduced by a factor 3 compared to the Schottky value $`2eI`$. If the length of the sample exceeds $`L_{in}`$, the electron gas can be described by a Fermi distribution with an effective temperature $`T_e`$ and Eq. (1) simplifies to:

$$S_I=\frac{4k_B}{RL}_{L/2}^{L/2}dxT_e(x)=\frac{4k_B\overline{T_e}}{R}.$$
(2)
Electron-electron (e-e) and electron-phonon (e-ph) scattering affect the noise differently. On one hand, the e-e interaction redistributes the energy of the electron system and $`T_e`$ increases due to an increasing number of electronic states that contribute to the noise. On the other hand, the power injected by the bias current in the sample can be dissipated to the lattice through e-ph interaction and $`T_e`$ decreases.
In S/N/S junctions, coherent electron-hole pairs penetrate from the superconductor into the normal metal over a distance $`L_c=min(L_\varphi ,\xi _ϵ=\sqrt{\mathrm{}D/ϵ})`$ with $`ϵ=max(k_BT,eV)`$ and $`D`$ the diffusion constant of the normal metal. If the sample length is smaller than $`L_c`$ (coherent case), the phase coherence covers the entire normal region and the Josephson effects determine conductance and noise behavior at low voltage. In this case, the conductance exhibits clear subgap structures (SGS) and the noise is strongly enhanced compared to the normal case due to the coherent transfer of large charge quanta SNScoherent ; Hoss .
In this paper we consider the incoherent case where $`L>>L_c`$. This regime has been recently studied theoretically by Bezuglyi et al. Bezuglyi  and Nagaev Nagaev2001 . They show that the noise is enhanced compared to normal junctions (N/N/N) because of the confinement of the subgap electrons in the normal part between the superconducting electrodes. If no inelastic processes take place (“collisionless regime”) and in the zero temperature limit the noise increases linearly with the bias voltage:

$$S_I(V)=\frac{2}{3R}(eV+2\mathrm{\Delta }),$$
(3)
where $`\mathrm{\Delta }`$ is the gap of the superconductor.
A simplified model to illustrate this behavior is the following: an electron with an energy $`eV<<\mathrm{\Delta }`$ can not escape in the superconducting reservoirs due to the absence of electronic states in the gap. Instead, it enters the superconductor together with a second electron as a Cooper pair and a hole is retroreflected in the normal metal (Andreev reflection). The reflected hole travels the normal region a second time and is retroreflected as an electron at the other S/N interface and so on. In the incoherent case the phase information between two subsequent Andreev reflections is lost. Therefore the quasiparticles experience Incoherent Multiple Andreev Reflections (IMAR). During each travel across the junction, the gain in energy is equal to $`eV`$, where $`V`$ is the applied voltage. Therefore the quasiparticles travel the normal part of the junction $`N`$ times with $`N=int[2\mathrm{\Delta }/(eV)+1]`$ before acquiring enough energy to escape in the superconducting electrode. Within this description, the total noise is the shot noise in a diffusive normal metal at zero temperature $`\frac{1}{3}2eI`$ times $`N`$.

At low voltage and finite temperature, the effective length of the junction for the multiply retroreflected particles $`L_{eff}=NLL\mathrm{\Delta }/V`$ exceeds the inelastic length $`L_{in}`$. In this “interacting regime” e-e-collisions interrupt the IMAR cycle before the quasiparticles reach the gap. In the case of strong interaction a Fermi distribution with an effective temperature $`T_e`$ is restored. $`T_e`$ decreases with decreasing voltage and reaches equilibrium (lattice temperature) at zero bias. Simultaneously the noise drops from the strongly enhanced level described by Eq. (3) to the Johnson-Nyquist level. Contrary to the N/N/N-case, e-e-interactions reduce the energy window of the involved electronic states.
Note that the voltage dependence of the effective length $`L_{eff}`$ in S/N/S junctions provides a unique way to study inelastic interactions in a normal metal because the same sample can be tuned from the strongly thermalized regime to the collisionless regime simply by changing the applied voltage.
On the experimental side only a few results are reported on current noise in incoherent diffusive S/N/S junctions at present. Besides measurements on short (coherent) junctions, Hoss et al. Hoss  also addressed the incoherent case. They studied $`2\mu m`$ long Nb/Au/Nb junctions in the interacting regime over the entire voltage range because of the large superconducting gap of Nb. In a paper by Roche et al. Roche  the length of the normal part is also about $`2\mu m`$, and despite the use of a high mobility 2DEG as the normal part, the e-ph scattering rate is high enough to drive the junction in the interacting regime. Finally Jehl et al. studied Nb/Al/Nb junctions at relatively high temperatures which behave as two S/N contacts in series JehlPRL . Using a different technique Pierre et al. measured directly the distribution function in long S/N/S junction but only at energies above the gap PierrePRL .
Whereas all these works investigated contacts with a length $`L5\mu m`$, we deliberately choose to study longer junctions ($`4,10`$ and $`60\mu m`$) to determine the respective role of e-e collisions, e-ph collisions and heat transfer through the S/N interfaces on the shot noise in S/N/S junctions. Moreover, to see the crossover from the interacting regime to the collisionless regime, we used aluminium because of its small superconducting gap. Because of the very small phase coherence length in the normal metal we used (see below), we have an almost perfect realization of the regime of IMAR.
To measure the current fluctuations we used a SQUID-based experimental setup JehlRSI . The intrinsic noise is about $`10\mu \mathrm{\Phi }_0/\sqrt{Hz}`$ which is equivalent to $`2pA/\sqrt{Hz}`$ in the input coil of the SQUID. The same experimental setup has been used to perform transport measurements.
Samples were fabricated by DC magnetron sputtering and optical lithography. First, a bilayer of Al/Cu ($`130nm/30nm`$) is deposited in situ to ensure a good contact between the two metals. This bilayer is etched to define the electrodes. Then the copper bridge (thickness $`130nm`$, purity of the Cu target: $`99.9999\%`$) is deposited by lift-off, preceded by a short backsputtering to clean the copper surface. Finally the whole sample is etched to remove the copper film from the aluminium electrodes. The resulting thickness of the copper bridge is $`90nm`$. A typical sample and a schematic cross section are shown in Fig. 1.
We studied junctions with 3 different lengths $`L`$: $`4`$, $`10`$ and $`60\mu m`$, but with the same width/length ratio: $`w/L3/2`$. All samples originate from the same wafer. In spite of very different overlap surfaces between aluminium and copper (because all dimensions scale with the length $`L`$), all samples have roughly the same resistance of about $`0.65\mathrm{\Omega }`$, which indicates a good interface (small barrier resistance). Below the transition of the Al electrodes at $`1.5K`$ and the transition of the Al/Cu bilayer (overlap region) at about $`1K`$ (see right inset of Fig. 2), the measured resistance is therefore essentially that of the normal part of the junctions. We then deduce the diffusion constant in copper: $`D=30cm^2s^1`$. We confirmed this value by measuring the resistance of a meander line consisting of $`700`$ squares in series, which was cosputtered on the same wafer.

In the temperature range from $`1K`$ to $`0.3K`$ the resistance decreases as expected in the classical proximity-effect. This fact and the reduced transition temperature of the bilayer are other indications of clean S/N interfaces. The resistance does not go down to zero since the sample lengths are much greater than the thermal length $`L_T=\sqrt{\frac{\mathrm{}D}{k_BT}}0.9\mu m`$ at $`T=30mK`$. Therefore no supercurrent is observed. On the contrary, at low temperatures ($`T<0.3K`$), the resistance increases again. This behavior is very similar to the reentrance in S/N junctions first observed by Charlat et al. Charlat . It means the two S/N interfaces of the junctions are connected incoherently. The temperature $`T_r`$ at which the resistance behavior changes is related to the phase breaking length $`L_\varphi `$ by $`k_BT_r\mathrm{}D/L_\varphi ^2`$ Courtois  (leading to $`L_\varphi 0.3\mu m`$) and is independent of $`L`$ (see left inset of Fig. 2). Such a large phase breaking rate is usual for Cu layers probably because of paramagnetic centers in the Cu oxide at the surface Vranken ; Gougam . The relative amplitude of the reentrance increases as the length $`L`$ decreases since the relative volume of the sample which is affected by coherent pairs increases.
The voltage dependence of the differential resistance $`dV/dI`$ measured at $`T=100mK`$ (Fig. 2) confirms the reentrance behavior. As for the temperature dependence, the amplitude and width of the reentrance peak depend on the width of the sample. However, we clearly see that the voltage needed to destroy the effect ($`2\mu V`$ and $`4\mu V`$ for the $`4\mu m`$ and $`10\mu m`$ samples) is much smaller than $`\frac{k_B}{e}T_r`$. As suggested in other experiments Courtois ; Black , this apparent discrepancy is due to a heating effect. Because we performed conductance and noise measurements simultaneously, we know the electronic temperature at low voltage. Anticipating the detailed description of noise measurements below, we found that the electron gas reaches a temperature of $`0.23K`$ at $`2\mu V`$ for the $`4\mu m`$ sample and $`0.2K`$ at $`4\mu V`$ for the $`10\mu m`$ sample. Consequently the rapid destruction of the reentrance peak is essentially due to the dramatic increase of $`T_e`$ at low bias voltage.
At higher voltage ($`V70\mu V`$) the bilayer (or at least part of it) is driven in the normal state, and a large peak occurs in the differential resistance (not shown). $`dV/dI`$ changes by only approximately $`10\%`$ up to $`V=70\mu V`$ and over this voltage range, we do not observe subharmonic gap structures (SGS) as expected in our experimental situation.
SGS appear at voltages $`V=2\mathrm{\Delta }/ne`$ due to the singularity in the DOS of the superconductor at the gap edges and, especially for junctions with high transparent interfaces, to the strongly enhanced probability of Andreev reflection at low energy. The latter is caused by the proximity effect in the normal metal Bezuglyi2000 . In our samples the DOS singularities are smeared by the bilayer structure of the electrodes (Al/Cu with good interface) and the proximity corrections are very small because of the very short correlation length $`L_c`$. Moreover the SGS of high order ($`n4`$ because $`70\mu eV\mathrm{\Delta }/2`$) are usually very weak and additionally smeared out by inelastic scattering at low voltage.
The results of the noise measurements at $`T=100mK`$ are shown in Fig. 3 where we have plotted the current noise density $`S_I`$ times the resistance $`R=V/I`$ versus DC voltage. The noise increases more rapidly with the bias voltage than what is expected for two independent S/N junctions in series with a reservoir in between (see dot-dashed line in Fig. 3), and we do not see the thermal crossover towards the Johnson-Nyquist noise level at $`eVk_BT9\mu V`$. The noise enhancement in these S/N/S junctions is due to the confinement of the quasiparticles between the two superconducting electrodes. With increasing sample length the confinement is relaxed by e-ph interaction and the noise slope at low voltage becomes less important.
At $`V50\mu V`$, the noise of the $`4\mu m`$ sample approaches a straight line compatible with the prediction in the collisionless regime (Eq. (3)) with $`\mathrm{\Delta }=135\mu eV`$. The gap value is reduced compared to the pure Al film ($`\mathrm{\Delta }200\mu eV`$) because of the bilayer structure of the superconducting electrodes. The reduction factor is the same as for the transition temperature.

At higher voltages ($`V70\mu V`$) the noise shows an irregular behavior which reflects the transition of the bilayer. In this voltage range, all IMAR cycles end by injecting a quasiparticle above the gap in the superconductor. Consequently, many quasiparticles arrive in the superconducting electrodes and weaken the superconductivity. Therefore, the voltage driven transition of the bilayer could be related to the collisionless regime itself.
In the following we focus on the interacting regime at low voltage where the electron gas can be described by a Fermi distribution function with an effective temperature $`T_e`$. Unlike normal junctions where the injected power can be evacuated in the reservoirs, in S/N/S junctions only quasiparticles with energy above the gap can dissipate power. When the sample length is long enough, some heat can be transferred to the phonons through e-ph-scattering. These two mechanisms are covered with our experiments since we studied long to very long junctions.
For the longest junction ($`60\mu m`$) we suppose cooling by phonons to be dominant. The electron temperature $`T_e`$ is then nearly constant over the whole sample length and can be calculated by a heat-diffusion equation Henny ; Wellstood . It yields:

$$T_e=(\frac{P}{\mathrm{\Sigma }\mathrm{\Omega }}+T_{ph}^5)^{1/5},$$
(4)
where $`P`$ is the power injected in the sample, $`\mathrm{\Omega }`$ its volume and $`T_{ph}`$ the temperature of the phonon bath (equal to the temperature of the mixing chamber). The parameter $`\mathrm{\Sigma }`$ is the e-ph coupling constant Wellstood . The noise is then given by $`4k_BT_e/R`$ and we obtain the best fit (see Fig. 3 and 4) for $`\mathrm{\Sigma }=\mathrm{2.4\hspace{0.17em}10}^9Wm^3K^5`$, a value which is of the same order of magnitude as reported for Au, Ag, Cu and AuCu Henny ; PierrePRL ; Roukes ; Wellstood .

For the shortest junctions ($`4\mu m`$) e-ph-scattering is not very efficient because the volume $`\mathrm{\Omega }`$ is much smaller and the main cooling mechanism is the heat transfer by quasiparticles through the S/N interfaces outside the gap region. Following Bezuglyi et al. Bezuglyi  the noise at low voltage is given by the Nyquist formula with a cut off at $`T_e`$ of the order of $`\mathrm{\Delta }`$:

$$S_I=\frac{4k_BT_e}{R}(12exp[\frac{\mathrm{\Delta }}{k_BT_e}]).$$
(5)
Here the reservoirs are supposed at $`T=0K`$. In fact, we can neglect the finite temperature of the electrodes because the quasiparticle temperature is rapidly much larger than the base temperature. From the measured noise $`S_I(V)`$ we obtain the effective electron temperature in the normal metal as a function of voltage $`T_e(V)`$ solving Eq. (5) numerically with $`\mathrm{\Delta }=135\mu eV`$.
Fig. 4 shows the total power $`V^2/R`$ as a function of the effective temperature of the quasiparticles. For the $`4\mu m`$ sample, we then compare the experimental results to the theoretical prediction obtained by Bezuglyi et al. Bezuglyi  who derived the power dissipated through the S/N interfaces as a function of the effective electron temperature at small voltage $`eV<<\mathrm{\Delta }`$ in good agreement with numerical simulations by Nagaev Nagaev2001 . The fit in Fig. 4 (thick solid line) is given by:

$$P_{NS}\frac{V^2}{R}=\frac{k_BT_e\mathrm{\Delta }}{e^2W_ϵR}\left(1+\frac{k_BT_e}{\mathrm{\Delta }}\right)exp[\frac{\mathrm{\Delta }}{k_BT_e}].$$
(6)
with $`W_ϵ=\tau _{ee}(\mathrm{\Delta })/\tau _D`$ where $`\tau _{ee}(\mathrm{\Delta })`$ is the e-e-scattering time at the gap energy and $`\tau _D`$ the diffusion time. We obtain good agreement between the experimental data and the prediction with $`W_ϵ=0.7\pm 0.2`$. This value is of the same order as the theoretical estimation of $`W_ϵ`$ ($`W_ϵ=2.3`$) using the standard theory of e-e-interaction in a 2D geometry Altshuler . A small but finite interface resistance would renormalize this estimation Bezuglyi2000  and could give a better agreement. The fact that $`W_ϵ1`$ indicates that the e-e-scattering time at the gap energy is of the order of the diffusion time $`\tau _D5ns`$ and is two orders of magnitude larger than the phase coherence time found above ($`L_\varphi 0.3\mu m`$ corresponds to $`\tau _\varphi 0.03ns`$).
Note that the theoretical model applied here is especially dedicated to incoherent S/N/S junctions and therefore more appropriate than the usually used Blonder-Tinkham-Klapwijk (BTK) model BTK ; Hoss ; QuirionJLT  and its derivatives Bardas .
For comparison we also plotted in Fig. 4 the power $`P_{ph}`$ dissipated by the phonons in the $`4\mu m`$ sample using the value of $`\mathrm{\Sigma }`$ obtained for the $`60\mu m`$ sample (thin solid line). In the range $`200`$ to $`600mK`$, $`P_{ph}`$ is about five times smaller than $`P_{NS}`$ and the error that we make neglecting this contribution is covered by the uncertainty on $`W_ϵ`$ given above.
Concerning the intermediate sample of length $`10\mu m`$ the power dissipated by phonons and through the interfaces is of the same order of magnitude. Extracting $`T_e`$ from Eq. (5), we fit the total power by adding the contributions of the two cooling mechanisms, treated separately according to Eq. (4) and (6). We obtain good agreement (see dashed line in Fig. 4) with the following parameters: $`\mathrm{\Delta }=135\mu eV`$, the same value of $`\mathrm{\Sigma }`$ as for the $`60\mu m`$ sample and $`W_ϵ=0.3\pm 0.1`$ (theoretical estimation $`W_ϵ0.4`$).
In conclusion, we investigated IMAR enhanced current noise in long S/N/S junctions of very different lengths ($`4,10`$ and $`60\mu m`$). We found that the noise temperature increases very rapidly at low voltage. We deduce the energy dependence of the thermal conductivity of the S/N interfaces which is in good agreement with recent semiclassical theory Bezuglyi ; Nagaev2001 . The noise behavior of the longest sample can be well fitted taking into account only phonon cooling. The inelastic scattering times we deduced are in agreement with standard description of e-e and e-ph interaction. With the same experimental setup we performed transport measurements. They reveal that in all the samples the two S/N interfaces are connected incoherently, indicating $`L_\varphi <<L_{in}`$. The complementarity of transport and noise measurements provided a direct analysis of the voltage dependence of the reentrance in terms of an effective electron temperature.

###### Acknowledgements.

We would like to acknowledge E.V. Bezuglyi and V.S. Shumeiko for valuable discussions and J.L. Thomassin for technical support.