text:
# PLANET observations of microlensing event OGLE-1999-BUL-23: limb darkening measurement of the source star

## 1 Introduction

In point-source-point-lens (PSPL) microlensing events, the light curve yields only one physically interesting parameter, the characteristic time scale of the event, $`t_\mathrm{E}`$, which is a combination of the mass of the lens and the source-lens relative parallax and proper motion. However, more varieties than PSPL events have been observed in reality, and using deviations from the standard light curve, one can deduce more information about the lens and the source. The Probing Lensing Anomalies NETwork (PLANET) is an international collaboration that monitors events in search of such anomalous light curves using a network of telescopes in the southern hemisphere (Albrow et al., 1998).

One example of information that can be extracted from anomalous events is the surface brightness profile of the source star (Witt, 1995). In a binary or multiple lens system, the caustic is an extended structure. If the source passes near or across the caustic, drastic changes in magnification near the caustics can reveal the finite size of the source (Gould, 1994; Nemiroff & Wickramasinghe, 1994; Witt & Mao, 1994; Alcock et al., 1997), and one can even extract its surface-brightness profile (Bogdanov & Cherepashchuk, 1996; Gould & Welch, 1996; Sasselov, 1997; Valls-Gabaud, 1998).

The fall-off of the surface brightness near the edge of the stellar disk with respect to its center, known as limb darkening, has been extensively observed in the Sun. Theories of stellar atmospheres predict limb darkening as a general phenomenon and give models for different types of stars. Therefore, measurement of limb darkening in distant stars other than the Sun would provide important observational constraints on the study of stellar atmospheres. However, such measurements are very challenging with traditional techniques and have usually been restricted to relatively nearby stars or extremely large supergiants. As a result, only a few attempts have been made to measure limb darkening to date. The classical method of tracing the stellar surface brightness profile is the analysis of the light curves of eclipsing binaries (Wilson & Devinney, 1971; Twigg & Rafert, 1980). However, the current practice in eclipsing-binary studies usually takes the opposite approach to limb darkening (Claret, 1998a) – constructing models of light curves using theoretical predictions of limb darkening. This came to dominate after Popper (1984) demonstrated that the uncertainty of limb darkening measurements from eclipsing binaries is substantially larger than the theoretical uncertainty. Since the limb-darkening parameter is highly correlated with other parameters of the eclipsing binary, fitting for limb darkening could seriously degrade the measurement of these other parameters. Multi-aperture interferometry and lunar occultation, which began as measurements of the angular sizes of stars, have also been used to resolve the surface structures of stars (Hofmann & Scholz, 1998). In particular, a large wavelength dependence of the interferometric size of a stellar disk has been attributed to limb darkening, and higher order corrections to account for limb darkening have been widely adopted in the interferometric angular size measurement of stars. Several recent investigations using optical interferometry extending beyond the first null of the visibility function have indeed confirmed that the observed patterns of the visibility function contradict a uniform stellar disk model and favor a limb-darkened disk (Quirrenbach et al., 1996; Hajian et al., 1998) although these investigations have used a model prediction of limb darkening inferred from the surface temperature rather than deriving the limb darkening from the observations. However, at least in one case, Burns et al. (1997) used interferometric imaging to measure the stellar surface brightness profile with coefficients beyond the simple linear model. In addition, developments of high resolution direct imaging in the last decade using space telescopes (Gilliland & Dupree, 1996) or speckle imaging (Kluckers et al., 1997) have provided a more straightforward way of detecting stellar surface irregularities. However, most studies of this kind are still limited to a few extremely large supergiants, such as $`\alpha `$ Ori. Furthermore, they seem to be more sensitive to asymmetric surface structures such as spotting than to limb darkening.

By contrast, microlensing can produce limb-darkening measurements for distant stars with reasonable accuracy. To date, limb darkening (more precisely, a set of coefficients of a parametrized limb-darkened profile) has been measured for source stars in three events, two K giants in the Galactic bulge and an A dwarf in the Small Magellanic Cloud (SMC). MACHO 97-BLG-28 was a cusp-crossing event of a K giant source with extremely good data, permitting Albrow et al. (1999a) to make a two-coefficient (linear and square-root) measurement of limb darkening. Afonso et al. (2000) used data from five microlensing collaborations to measure linear limb darkening coefficients in five filter bandpasses for MACHO 98-SMC-1, a metal-poor A star in the SMC. Although the data for this event were also excellent, the measurement did not yield a two-parameter determination because the caustic crossing was a fold-caustic rather than a cusp, and these are less sensitive to the form of the stellar surface brightness profile. Albrow et al. (2000a) measured a linear limb-darkening coefficient for MACHO 97-BLG-41, a complex rotating-binary event with both a cusp crossing and a fold-caustic crossing. In principle, such an event could give very detailed information about the surface brightness profile. However, neither the cusp nor the fold-caustic crossing was densely sampled, so only a linear parameter could be extracted.

In this paper, we report a new limb-darkening measurement of a star in the Galactic bulge by a fold-caustic crossing event, OGLE-1999-BUL-23, based on the photometric monitoring of PLANET.

## 2 OGLE-1999-BUL-23

OGLE-1999-BUL-23 was originally discovered towards the Galactic bulge by the Optical Gravitational Lensing Experiment (OGLE) 111The OGLE alert for this event is posted at http://www.astrouw.edu.pl/~ftp/ogle/ogle2/ews/bul-23.html (Udalski et al., 1992; Udalski, Kubiak, & Szymański, 1997). The PLANET collaboration observed the event as a part of our routine monitoring program after the initial alert, and detected a sudden increase in brightness on 12 June 1999.222 The PLANET anomaly and caustic alerts are found at http://www.astro.rug.nl/~planet/OB99023cc.html Following this anomalous behavior, we began dense (typically one observation per hour) photometric sampling of the event. Since the source lies close to the (northern) winter solstice ($`\alpha =18^h07^m45\stackrel{\mathrm{s}}{\mathrm{.}}14`$, $`\delta =27\mathrm{°}33\mathrm{}15\stackrel{}{\mathrm{.}}4`$), while the caustic crossing occurred nearly at the summer solstice (19 June 1999), and since good weather at all four of our southern sites prevailed throughout, we were able to obtain nearly continuous coverage of the second caustic crossing without any significant gaps. Visual inspection and initial analysis of the light curve revealed that the second crossing was due to a simple fold caustic crossing (see §2.2).

### 2.1 Data

We observed OGLE-1999-BUL-23 with I and V band filters at four participant telescopes: the Elizabeth 1 m at South African Astronomical Observatory (SAAO), Sutherland, South Africa; the Perth/Lowell 0.6 m telescope at Perth, Western Australia; the Canopus 1 m near Hobart, Tasmania, Australia; and the Yale/AURA/Lisbon/OSU 1 m at Cerro Tololo Inter-American Observatory (CTIO), La Serena, Chile. From June to August 1999 ($`1338<\mathrm{HJD}^{}<1405`$), PLANET obtained almost 600 images of the field of OGLE-1999-BUL-23. In addition, baseline points were taken at SAAO ($`\mathrm{HJD}^{}1440`$) and Perth ($`\mathrm{HJD}^{}1450`$; $`\mathrm{HJD}^{}1470`$). Here $`\mathrm{HJD}^{}\mathrm{HJD}2450000`$, where HJD is Heliocentric Julian Date at center of exposure. The data reduction and photometric measurements of the event were performed relative to non-variable stars in the same field using DoPHOT. After several re-reductions, we recovered the photometric measurements from a total of 476 frames.

We assumed independent photometric systems for different observatories and thus explicitly included the determination of independent (unlensed) source and background fluxes for each different telescope and filter band in the analysis. This provides both determinations of the photometric offsets between different systems and independent estimates of the blending factors. The final results demonstrate satisfactory alignment among the data sets (see §2.3), and we therefore believe that we have reasonable relative calibrations. Our previous studies have shown that the background flux (or blending factors) may correlate with the size of seeing disks in some cases (Albrow et al., 2000a, b). To check this, we introduced linear seeing corrections in addition to constant backgrounds.

From previous experience, it is expected that the formal errors reported by DoPHOT underestimate the actual errors (Albrow et al., 1998), and consequently that $`\chi ^2`$ is overestimated. Hence, we renormalize photometric errors to force the final reduced $`\chi ^2/dof=1`$ for our best fit model. Here, $`dof`$ is the number of degrees of freedom (the number of data points less the number of parameters). We determine independent rescaling factors for the photometric uncertainties from the different observatories and filters. The process involves two steps: the elimination of bad data points and the determination of error normalization factors. In this as in all previous events that we have analyzed, there are outliers discrepant by many $`\sigma `$ that cannot be attributed to any specific cause even after we eliminate some points whose source of discrepancy is identifiable. Although, in principle, whether particular data points are faulty or not should be determined without any reference to models, we find that the light curves of various models that yield reasonably good fits to the data are very similar to one another, and furthermore, there is no indication of temporal clumping of highly discrepant points. We, therefore, identify outlier points with respect to our best model and exclude them from the final analysis.

For the determination of outliers, we follow an iterative approach using both steps of error normalization. First, we calculate the individual $`\chi ^2`$’s of data sets from different observatories and filter bands with reference to our best model without any rejection or error scaling. Then, the initial normalization factors are determined independently for each data set using those individual $`\chi ^2`$’s and the number of data points in each set. If the deviation of the most discrepant outlier is larger than what is predicted based on the number of points and the assumption of a normal distribution, we classify the point as bad and calculate the new $`\chi ^2`$’s and the normalization factors again. We repeat this procedure until the largest outlier is comparable with the prediction of a normal distribution. Although the procedure appears somewhat arbitrary, the actual result indicates that there exist rather large decreases of $`\sigma `$ between the last rejected and included data points. After rejection of bad points, 429 points remain (see Table 1 and Fig. 1).

### 2.2 Analysis: searching for $`\chi ^2`$ minima

We use the method of Albrow et al. (1999b, hereafter Paper I), which was devised to fit the light curve of fold-caustic crossing binary-lens events, to analyze the light curve of this event and find an appropriate binary-lens solution. This method consists of three steps: (1) fitting of caustic-crossing data using an analytic approximation of the magnification, (2) searching for $`\chi ^2`$ minima over the whole parameter space using the point-source approximation and restricted to the non-caustic crossing data, and (3) $`\chi ^2`$ minimization using all data and the full binary-lens equation in the neighborhood of the minima found in the second step.

For the first step, we fit the I-band caustic crossing data ($`1348.5\mathrm{HJD}^{}1350`$) to the six-parameter analytic curve shown in equation (1) that characterizes the shape of the second caustic crossing (Paper I; Afonso et al. 2000),

$$F(t)=\left(\frac{Q}{\mathrm{\Delta }t}\right)^{1/2}\left[G_0\left(\frac{tt_{\mathrm{cc}}}{\mathrm{\Delta }t}\right)+\mathrm{\Gamma }H_{1/2}\left(\frac{tt_{\mathrm{cc}}}{\mathrm{\Delta }t}\right)\right]+F_{\mathrm{cc}}+(tt_{\mathrm{cc}})\stackrel{~}{\omega },$$
(1a)
$$G_n(\eta )\pi ^{1/2}\frac{(n+1)!}{(n+1/2)!}_{\mathrm{max}(\eta ,1)}^1𝑑x\frac{(1x^2)^{n+1/2}}{(x\eta )^{1/2}}\mathrm{\Theta }(1\eta ),$$
(1b)
$$H_{1/2}(\eta )G_{1/2}(\eta )G_0(\eta ).$$
(1c)
Figure 2 shows the best-fit curve and the data points used for the fit. This caustic-crossing fit essentially constrains the search for a full solution to a four-dimensional hypersurface instead of the whole nine-dimensional parameter space (Paper I).

We construct a grid of point-source light curves with model parameters spanning a large subset of the hypersurface and calculate $`\chi ^2`$ for each model using the I-band non-caustic crossing data. After an extensive search for $`\chi ^2`$-minima over the four-dimensional hypersurface, we find positions of two apparent local minima, each in a local valley of the $`\chi ^2`$-surface. The smaller $`\chi ^2`$ of the two is found at $`(d,q,\alpha )(2.4,\mathrm{\hspace{0.17em}0.4},\mathrm{\hspace{0.17em}75}\mathrm{°})`$, where $`d`$ is the projected binary separation in units of the Einstein ring radius, $`q`$ is the mass ratio of the binary system, and $`\alpha `$ is the angle between the binary axis and the path of the source, defined so that the geometric center of the lens system lies on the right hand side of the moving source. The other local minimum is $`(d,q,\alpha )(0.55,\mathrm{\hspace{0.17em}0.55},\mathrm{\hspace{0.17em}260}\mathrm{°})`$. The results appear to suggest a rough symmetry of $`dd^1`$ and $`(\alpha <\pi )(\alpha >\pi )`$, as was found for MACHO 98-SMC-1 (Paper I; Afonso et al. 2000). Besides these two local minima, there are several isolated $`(d,q)`$-grid points at which $`\chi ^2`$ is smaller than at neighboring grid points. However, on a finer grid they appear to be connected with one of the two local minima specified above. We include the two local minima and some of the apparently isolated minimum points as well as points in the local valley around the minima as starting points for the refined search of $`\chi ^2`$-minimization in the next step.

### 2.3 Solutions: $`\chi ^2`$ minimization

Starting from the local minima found in §2.2 and the points in the local valleys around them, we perform a refined search for the $`\chi ^2`$ minimum. The $`\chi ^2`$ minimization includes all the I and V data points for successive fitting to the full expression for magnification, accounting for effects of a finite source size and limb darkening.

As described in Paper I, the third step makes use of a variant of equation (1) to evaluate the magnified flux in the neighborhood of the caustic crossing. Paper I found that, for MACHO 98-SMC-1, this analytic expression was an extremely good approximation to the results of numerical integration and assumed that the same would be the case for any fold crossing. Unfortunately, we find that, for OGLE-1999-BUL-23, this approximation deviates from the true magnification as determined using the method of Gould & Gaucherel (1997) as much as 4%, which is larger than our typical photometric uncertainty in the region of caustic crossing. To maintain the computational efficiency of Paper I, we continue to use the analytic formula (1), but correct it by pre-tabulated amounts given by the fractional difference (evaluated close to the best solution) between this approximation and the values found by numerical integration. We find that this correction works quite well even at the local minimum for the other (close-binary) solution – the error is smaller than 1%, and in particular, the calculations agree within 0.2% for the region of primary interest. The typical (median) photometric uncertainties for the same region are 0.015 mag (Canopus after the error normalization) and 0.020 mag (Perth). In addition, we test the correction by running the fitting program with the exact calculation at the minimum found using the corrected approximation, and find that the measured parameters change less than the precision of the measurement. In particular, the limb-darkening coefficients change by an order of magnitude less than the measurement uncertainty due to the photometric errors.

The results of the refined $`\chi ^2`$ minimization are listed in Table 2 for three discrete “solutions” and in Table 3 for grid points neighboring the best-fit solution whose $`\mathrm{\Delta }\chi ^2`$ is less than one. The first seven columns describe seven of the standard parameters of the binary-lens model (the remaining two parameters are the source and background flux). The eighth column is the time of the second caustic crossing ($`t_{\mathrm{cc}}`$) – the time when the center of the source crossed the caustic. The limb darkening coefficients for I and V bands are shown in the next two columns. The final column is $`\mathrm{\Delta }\chi ^2`$,

$$\mathrm{\Delta }\chi ^2\frac{\chi ^2\chi _{\mathrm{best}}^2}{\chi _{\mathrm{best}}^2/dof},$$
(2)
as in Paper I. The light curve (in magnification) of the best-fit model is shown in Figure 3 together with all the data points used in the analysis.

#### 2.3.1 “Degeneracy”

For typical binary-lens microlensing events, more than one solution often fits the observations reasonably well. In particular, Dominik (1999) predicted a degeneracy between close and wide binary lenses resulting from a symmetry in the lens equation itself, and such a degeneracy was found empirically for MACHO 98-SMC-1 (Paper I; Afonso et al. 2000).

We also find two distinct local $`\chi ^2`$ minima (§2.2) that appear to be closely related to such degeneracies. However, in contrast to the case of MACHO 98-SMC-1, our close-binary model for OGLE-1999-BUL-23 has substantially higher $`\chi ^2`$ than the wide-binary model ($`\mathrm{\Delta }\chi ^2=\mathrm{\hspace{0.17em}127.86}`$). Figure 4 shows the predicted light curves in SAAO instrumental I band. The overall geometries of these two models are shown in Figures 5 and 6. The similar morphologies of the caustics with respect to the path of the source is responsible for the degenerate light curves near the caustic crossing (Fig. 6). However, the close-binary model requires a higher blending fraction and lower baseline flux than the wide-binary solution because the former displays a higher peak magnification ($`A_{\mathrm{max}}50`$ vs. $`A_{\mathrm{max}}30`$). Consequently, a precise determination of the baseline can significantly contribute to discrimination between the two models, and in fact, the actual data did constrain the baseline well enough to produce a large difference in $`\chi ^2`$.

A fair number of pre-event baseline measurements are available via OGLE, and those data can further help discriminate between these two “degenerate” models. We fit OGLE measurements to the two models with all the model parameters being fixed and allowing only the baseline and the blending fraction as free parameters. We find that the PLANET wide-binary model produces $`\chi ^2=306.83`$ for 169 OGLE points ($`\chi ^2/dof=1.83`$, compare Table 1) while $`\chi ^2=608.22`$ for the close-binary model for the same 169 points (Fig. 7). That is, $`\mathrm{\Delta }\chi ^2=164.04`$, so that the addition of OGLE data by itself discriminates between the two models approximately as well as all the PLANET data combined. The largest contribution to this large $`\mathrm{\Delta }\chi ^2`$ appears to come from the period about a month before the first caustic crossing which is well covered by the OGLE data but not by the PLANET data. In particular, the close-binary model predicts a bump in the light curve around $`\mathrm{HJD}^{}1290`$ due to a triangular caustic (see Fig. 5), but the data do not show any abnormal feature in the same region, although it is possible that rotation of the binary moved the caustic far from the source trajectory (eg. Afonso et al. 2000). In brief, the OGLE data strongly favor the wide-binary model.

#### 2.3.2 Limb-Darkening Coefficients

The limb darkening of the source is parametrized using a normalized linear model of the source surface brightness profile, which was introduced in Appendix B of Paper I,

$`S_\lambda (\vartheta )=\overline{S}_\lambda \left[1\mathrm{\Gamma }_\lambda \left(1{\displaystyle \frac{3}{2}}\mathrm{cos}\vartheta \right)\right]=\overline{S}_\lambda \left[(1\mathrm{\Gamma }_\lambda )+{\displaystyle \frac{3}{2}}\mathrm{\Gamma }_\lambda \mathrm{cos}\vartheta \right],`$
$`\mathrm{where}\mathrm{sin}\vartheta {\displaystyle \frac{\theta }{\theta _{}}}\mathrm{and}\overline{S}_\lambda {\displaystyle \frac{F_{s,\lambda }}{\pi \theta _{}^2}},`$ (3)
while linear limb darkening is usually parametrized by,

$$S_\lambda (\vartheta )=S_\lambda (0)\left[1c_\lambda (1\mathrm{cos}\vartheta )\right]=S_\lambda (0)\left[(1c_\lambda )+c_\lambda \mathrm{cos}\vartheta \right].$$
(4)
The relationship between the two expressions of linear limb-darkening coefficients is

$$c_\lambda =\frac{3\mathrm{\Gamma }_\lambda }{2+\mathrm{\Gamma }_\lambda }.$$
(5)
Amongst our six data sets, data from SAAO did not contain points that were affected by limb darkening, i.e. caustic crossing points. Since the filters used at different PLANET observatories do not differ significantly from one another, we use the same limb-darkening coefficient for the three remaining I-band data sets. The V-band coefficient is determined only from Canopus data, so that a single coefficient is used automatically.

For the best-fit lens geometry, the measured values of linear limb-darkening coefficients are $`\mathrm{\Gamma }_I=\mathrm{\hspace{0.17em}0.534}\pm 0.020`$ and $`\mathrm{\Gamma }_V=\mathrm{\hspace{0.17em}0.711}\pm 0.089`$, where the errors include only uncertainties in the linear fit due to the photometric uncertainties at *fixed* binary-lens model parameters. However, these errors underestimate the actual uncertainties of the measurements because the measurements are correlated with the determination of the seven lens parameters shown in Tables 2 and 3. Incorporating these additional uncertainties in the measurement (see the next section for a detailed discussion of the error determination), our final estimates are

$$\mathrm{\Gamma }_I=\mathrm{\hspace{0.17em}0.534}_{0.040}^{+0.050}\left(c_I=\mathrm{\hspace{0.17em}0.632}_{0.037}^{+0.047}\right),$$
(6a)
$$\mathrm{\Gamma }_V=\mathrm{\hspace{0.17em}0.711}_{0.095}^{+0.098}\left(c_V=\mathrm{\hspace{0.17em}0.786}_{0.078}^{+0.080}\right).$$
(6b)
This is consistent with the result of the caustic-crossing fit of §2.2 ($`\mathrm{\Gamma }_I=\mathrm{\hspace{0.17em}0.519}\pm 0.043`$). Our result suggests that the source is more limb-darkened in $`V`$ than in $`I`$, which is generally expected by theories. Figure 8 shows the I-band residuals (in magnitudes) at the second caustic crossing from our best-fit models for a linearly limb-darkened and a uniform disk model. It is clear that the uniform disk model exhibits larger systematic residuals near the peak than the linearly limb-darkened disk. From the residual patterns – the uniform disk model produces a shallower slope for the most of the falling side of the second caustic crossing than the data require, one can infer that the source should be more centrally concentrated than the model predicts, and consequently the presence of limb darkening. The linearly limb-darkened disk reduces the systematic residuals by a factor of $``$ 5. Formally, the difference of $`\chi ^2`$ between the two models is 172.8 with two additional parameters for the limb-darkened disk model, i.e. the data favor a limb-darkened disk over a uniform disk at very high confidence.

## 3 Error Estimation for Limb Darkening Coefficients

Due to the multi-parameter character of the fit, a measurement of any parameter is correlated with other parameters of the model. The limb-darkening coefficients obtained with the different model parameters shown in Table 3 exhibit a considerable scatter, and in particular, for the I-band measurement, the scatter is larger than the uncertainties due to the photometric errors. This indicates that, in the measurement of the limb-darkening coefficients, we need to examine errors that correlate with the lens model parameters in addition to the uncertainties resulting from the photometric uncertainties at fixed lens parameters. This conclusion is reinforced by the fact that the error in the estimate of $`\mathrm{\Gamma }`$ from the caustic-crossing fit (see Fig. 2), which includes the correlation with the parameters of the caustic-crossing, is substantially larger than the error in the linear fit, which does not.

Since limb darkening manifests itself mainly around the caustic crossing, its measurement is most strongly correlated with $`\mathrm{\Delta }t`$ and $`t_{\mathrm{cc}}`$. To estimate the effects of these correlations, we fit the data to models with $`\mathrm{\Delta }t`$ or $`t_{\mathrm{cc}}`$ fixed at several values near the best fit – the global geometry of the best fit, i.e. $`d`$ and $`q`$ being held fixed as well. The resulting distributions of $`\mathrm{\Delta }\chi ^2`$ have parabolic shapes as a function of the fit values of the limb-darkening coefficient and are centered at the measurement of the best fit. (Both, $`\mathrm{\Delta }t`$ fixed and $`t_{\mathrm{cc}}`$ fixed, produce essentially the same parabola, and therefore, we believe that the uncertainty related to each correlation with either $`\mathrm{\Delta }t`$ or $`t_{\mathrm{cc}}`$ is, in fact, same in its nature.) We interpret the half width of the parabola at $`\mathrm{\Delta }\chi ^2=1`$ ($`\delta \mathrm{\Gamma }_I=0.031,\delta \mathrm{\Gamma }_V=0.032`$) as the uncertainty due to the correlation with the caustic-crossing parameters at a given global lens geometry of a fixed $`d`$ and $`q`$.

Although the global lens geometry should not *directly* affect the limb darkening measurement, the overall correlation between local and global parameters can contribute an additional uncertainty to the measurement. This turns out to be the dominant source of the scatter found in Table 3. To incorporate this into our final determination of errors, we examine the varying range of the measured coefficients over $`\mathrm{\Delta }\chi ^21`$. The result is apparently asymmetric between the direction of increasing and decreasing the amounts of limb darkening. We believe that this is real, and thus we report asymmetric error bars for the limb-darkening measurements.

The final errors of the measurements reported in §2.3.2 are determined by adding these two sources of error to the photometric uncertainty in quadrature. The dominant source of errors in the I-band coefficient measurement is the correlation between the global geometry and the local parameters whereas the photometric uncertainty is the largest contribution to the uncertainties in the V-band coefficient measurement.

Although the measurements of V and I band limb darkening at fixed model parameters are independent, the final estimates of two coefficients are not actually independent for the same reason discussed above. (The correlation between V and I limb-darkening coefficients is clearly demonstrated in Table 3.) Hence, the complete description of the uncertainty requires a covariance matrix.

$$C=C_{\mathrm{phot}}+\stackrel{~}{C}_{\mathrm{cc}}^{1/2}\left(\begin{array}{cc}1& \xi \\ \xi & 1\end{array}\right)\stackrel{~}{C}_{\mathrm{cc}}^{1/2}+\stackrel{~}{C}_{\mathrm{geom}}^{1/2}\left(\begin{array}{cc}1& \xi \\ \xi & 1\end{array}\right)\stackrel{~}{C}_{\mathrm{geom}}^{1/2},$$
(7a)
$$C_{\mathrm{phot}}\left(\begin{array}{cc}\sigma _{V,\mathrm{phot}}^2& 0\\ 0& \sigma _{I,\mathrm{phot}}^2\end{array}\right),$$
(7b)
$$\stackrel{~}{C}_{\mathrm{cc}}^{1/2}\left(\begin{array}{cc}\sigma _{V,\mathrm{cc}}& 0\\ 0& \sigma _{I,\mathrm{cc}}\end{array}\right),$$
(7c)
$$\stackrel{~}{C}_{\mathrm{geom}}^{1/2}\left(\begin{array}{cc}\overline{\sigma }_{V,\mathrm{geom}}& 0\\ 0& \overline{\sigma }_{I,\mathrm{geom}}\end{array}\right),$$
(7d)
where the subscript (phot) denotes the uncertainties due to the photometric errors; (cc), the correlation with $`\mathrm{\Delta }t`$ and $`t_{\mathrm{cc}}`$ at a fixed $`d`$ and $`q`$; (geom), the correlation with the global geometry, and $`\xi `$ is the correlation coefficient between $`\mathrm{\Gamma }_V`$ and $`\mathrm{\Gamma }_I`$ measurement. We derive the correlation coefficient using each measurement of $`\mathrm{\Gamma }_V`$ and $`\mathrm{\Gamma }_I`$, and the result indicates that two measurements are almost perfectly correlated $`(\xi =0.995)`$. We accommodate asymmetry of the errors by making the error ellipse off-centered with respect to the best estimate. (See §5 for more discussion on the error ellipses.)

## 4 Physical Properties of the Source Star

Figure 9 shows color-magnitude diagrams (CMDs) derived from a $`2^{}\times 2^{}`$ SAAO field and a $`4^{}\times 4^{}`$ Canopus field centered on OGLE-1999-BUL-23 with positions marked for the unmagnified source (S), the baseline (B), blended light (BL) at median seeing, and the center of red clump giants (RC). The source position in these CMDs is consistent with a late G or early K subgiant in the Galactic bulge (see below). Using the color and magnitude of red clump giants in the Galactic bulge reported by Pacyński et al. (1999) ($`I_{\mathrm{RC}}=\mathrm{\hspace{0.17em}14.37}\pm 0.02,[VI]_{\mathrm{RC}}=\mathrm{\hspace{0.17em}1.114}\pm 0.003`$), we measure the reddening-corrected color and magnitude of the source in the Johnson-Cousins system from the relative position of the source with respect to the center of red clump in our CMDs, and obtain:

$$(VI)_{\mathrm{S},0}=\mathrm{\hspace{0.17em}1.021}\pm 0.044,$$
(8a)
$$V_{\mathrm{S},0}=\mathrm{\hspace{0.17em}18.00}\pm 0.06,$$
(8b)
where the errors include the difference of the source positions in the two CMDs, but may still be somewhat underestimated because the uncertainty in the selection of red clump giants in our CMDs has not been quantified exactly.

From this information, we derive the surface temperature of the source; $`T_{\mathrm{eff}}=(4830\pm 100)`$ K, using the color calibration in Bessell, Castelli, & Plez (1998) and assuming $`\mathrm{log}g=\mathrm{\hspace{0.17em}3.5}`$ and the solar abundance. This estimate of temperature is only weakly dependent on the assumed surface gravity and different stellar atmospheric models. To determine the angular size of the source, we use equation (4) of Albrow et al. (2000a), which is derived from the surface brightness-color relation of van Belle (1999). We first convert $`(VI)_{\mathrm{S},0}`$ into $`(VK)_{\mathrm{S},0}=\mathrm{\hspace{0.17em}2.298}\pm 0.113`$ using the same color calibration of Bessell et al. (1998) and then obtain the angular radius of the source of

$`\theta _{}`$ $`=`$ $`(1.86\pm 0.13)\mu \mathrm{as}`$ (9)
$`=`$ $`(0.401\pm 0.027)R_{\mathrm{}}\mathrm{kpc}^1.`$
If the source is at the Galactocentric distance (8 kpc), this implies that the radius of the source is roughly 3.2$`R_{\mathrm{}}`$, which is consistent with the size of a $`1M_{\mathrm{}}`$ subgiant ($`\mathrm{log}g=\mathrm{\hspace{0.17em}3.4}`$).

Combining this result with the parameters of the best-fit model yields

$`\mu =\theta _{}/(\mathrm{\Delta }t\mathrm{sin}\varphi )`$ $`=`$ $`(13.2\pm 0.9)\mu \mathrm{as}\mathrm{day}^1`$ (10)
$`=`$ $`(22.8\pm 1.5)\mathrm{km}\mathrm{s}^1\mathrm{kpc}^1,`$
$`\theta _\mathrm{E}=\mu t_\mathrm{E}`$ $`=`$ $`(0.634\pm 0.043)\mathrm{mas},`$ (11)
where $`\varphi =\mathrm{\hspace{0.17em}123}\stackrel{}{\mathrm{.}}9`$ is the angle that the source crosses the caustic (see Fig. 6). This corresponds to a projected relative velocity of $`(182\pm 12)\mathrm{km}\mathrm{s}^1`$ at the Galactocentric distance, which is generally consistent with what is expected in typical bulge/bulge or bulge/disk (source/lens) events, but inconsistent with disk/disk lensing. Hence we conclude that the source is in the bulge. As for properties of the lens, the projected separation of the binary lens is $`(1.53\pm 0.10)`$ AU $`\mathrm{kpc}^1`$, and the combined mass of the lens is given by

$$M_\mathrm{L}=\frac{c^2D_\mathrm{S}D_\mathrm{L}}{4G(D_\mathrm{S}D_\mathrm{L})}\theta _\mathrm{E}^2=(0.395\pm 0.053)\left(\frac{x}{1x}\right)\left(\frac{D_\mathrm{S}}{8\mathrm{kpc}}\right)M_{\mathrm{}},$$
(12)
where $`xD_\mathrm{L}/D_\mathrm{S}`$, $`D_\mathrm{L}`$ is the distance to the lens, and $`D_\mathrm{S}`$ is the distance to the source.

## 5 Limb Darkening of the Source

We compare our determination of the linear limb-darkening coefficients to model calculations by Claret, Díaz-Cordovés, & Giménez (1995) and Díaz-Cordovés, Claret, & Giménez (1995). For an effective temperature of $`T_{\mathrm{eff}}=(4830\pm 100)`$ K and a surface gravity of $`\mathrm{log}g=\mathrm{\hspace{0.17em}3.5}`$, the interpolation of the V-band linear limb-darkening coefficients, $`c_V`$, of Díaz-Cordovés et al. (1995) predicts a value $`c_V=\mathrm{\hspace{0.17em}0.790}\pm 0.012`$, very consistent with our measurement. However, for the I-band coefficient, the prediction of Claret et al. (1995), $`c_I=\mathrm{\hspace{0.17em}0.578}\pm 0.008`$, is only marginally consistent with our measurement, at the $`1.46\sigma `$ level. Adopting a slightly different gravity does not qualitatively change this general result. Since we believe that the uncertainty in the color of the source is larger than in the limb-darkening coefficients, we also examine the opposite approach to the theoretical calculations $``$ using the measured values of limb-darkening coefficients to derive the effective temperature of the source. If the source is a subgiant ($`\mathrm{log}g\mathrm{\hspace{0.17em}3.5}`$) as our CMDs suggest, the measured values of the limb-darkening coefficients are expected to be observed in stars of the effective temperature, $`T_{\mathrm{eff}}=(4850_{670}^{+650})`$ K for $`c_V`$ or $`T_{\mathrm{eff}}=(4200_{490}^{+390})`$ K for $`c_I`$. As before, the estimate from the V-band measurement shows a better agreement with the measured color than the estimate from the I-band. Considering that the data quality of I band is better than V band (the estimated uncertainty is smaller in $`I`$ than in $`V`$), this result needs to be explained.

In Figure 10, we plot theoretical calculations of $`(c_I,c_V)`$ together with our measured values. In addition to Díaz-Cordovés et al. (1995) and Claret et al. (1995) (A), we also include the calculations of linear limb-darkening coefficients by van Hamme (1993) (B) and Claret (1998b) (C). For all three calculations, the V-band linear coefficients are generally consistent with the measured coefficients and the color, although van Hamme (1993) predicts a slightly smaller amount of limb darkening than the others. On the other hand, the calculations of the I-band linear coefficients are somewhat smaller than the measurement except for Claret (1998b) with $`\mathrm{log}g=\mathrm{\hspace{0.17em}4.0}`$. (However, to be consistent with a higher surface gravity while maintaining its color, the source star should be in the disk, which is inconsistent with our inferred proper motion.) Since $`c_V`$ and $`c_I`$ are not independent (in both the theories and in our measurement), it is more reasonable to compare the I and V band measurements to the theories simultaneously. Using the covariance matrix of the measurement of $`\mathrm{\Gamma }_I`$ and $`\mathrm{\Gamma }_V`$ (see §3), we derive error ellipses for our measurements in the $`(c_I,c_V)`$ plane and plot them in Figure 10. Formally, at the $`1\sigma `$ level, the calculations of the linear limb-darkening coefficients in any of these models are not consistent with our measurements. In principle, one could also constrain the most likely stellar types that are consistent with the measured coefficients, independent of *a priori* information on the temperature and the gravity, with a reference to a model. If we do this, the result suggests either that the surface temperature is cooler than our previous estimate from the color or that the source is a low-mass main-sequence ($`\mathrm{log}g\mathrm{\hspace{0.17em}4.0}`$) star. However, the resulting constraints are not strong enough to place firm limits on the stellar type even if we assume any of these models to be “correct”.

One possible explanation of our general result – the measured V-band coefficients are nearly in perfect agreement with theories while the I-band coefficients are only marginally consistent – is non-linearity of stellar limb darkening. Many authors have pointed out the inadequacy of the linear limb darkening in producing a reasonably high-accuracy approximation of the real stellar surface brightness profile (Wade & Rucinski, 1985; Díaz-Cordovés & Giménez, 1992; van Hamme, 1993; Claret, 1998b). Indeed, Albrow et al. (1999a) measured the two-coefficient square-root limb darkening for a cusp-crossing microlensing event and found that the single-coefficient model gives a marginally poorer fit to the data. The quality of the linear parameterization has been investigated for most theoretical limb-darkening calculations, and the results seem to support this explanation. van Hamme (1993) defined the quality factors (Q in his paper) for his calculations of limb-darkening coefficients, and for 4000 K $`T_{\mathrm{eff}}`$ 5000 K and 3.0 $`\mathrm{log}g`$ 4.0, his results indicate that the linear parameterization is a better approximation for V band than for I band. Similarly, Claret (1998b) provided plots of summed residuals ($`\sigma `$ in his paper) for his fits used to derive limb-darkening coefficients showing that the V-band linear limb-darkening has lower $`\sigma `$ than I-band and is as good as the V-band square-root limb-darkening near the temperature range of our estimate for the source of OGLE-1999-BUL-23. In fact, Díaz-Cordovés et al. (1995) reported that the V-band limb darkening is closest to the linear law in the temperature range $`T_{\mathrm{eff}}=\mathrm{\hspace{0.17em}4500}4750`$ K. In summary, the source happens to be very close to the temperature at which the linear limb darkening is a very good approximation in $`V`$, but is less good in $`I`$.

The actual value of the coefficient in the linear parameterization of a non-linear profile may vary depending on the method of calculation and sampling. In order to determine the linear coefficients, models (A) and (C) used a least square fit to the theoretical (non-parametric) profile by sampling uniformly over $`\mathrm{cos}\vartheta `$ (see eq. ), while model (B) utilized the principle of total flux conservation between parametric and non-parametric profiles. On the other hand, a fold-caustic crossing event samples the stellar surface brightness by convolving it with a rather complicated magnification pattern (Gaudi & Gould, 1999). Therefore, it is very likely that neither of the above samplings and calculations is entirely suitable for the representation of the limb-darkening measurement by microlensing unless the real intensity profile of the star is actually same as the assumed parametric form (the linear parameterization, in this case). In fact, the most apropriate way to compare the measurement to the stellar atmospheric models would be a direct fit to the (non-parametric) theoretical profile after convolution with the magnification patterns near the caustics. In the present paper, this has not been done, but we hope to make such a direct comparison in the future.

We thank A. Udalski for re-reducing the OGLE data on very short notice after we noticed an apparent discrepancy between the PLANET data and the original OGLE reductions. This work was supported by grants AST 97-27520 and AST 95-30619 from the NSF, by grant NAG5-7589 from NASA, by a grant from the Dutch ASTRON foundation through ASTRON 781.76.018, by a Marie Curie Fellowship from the European Union, and by “coup de pouce 1999” award from the Ministère de l’Éducation nationale, de la Recherche et de la Technologie.


format: text

text:
# No CPT Violation from Tilted Brane in Neutral Meson-Antimeson Systems

## Acknowledgement

The work of G.C. was supported by the KOSEF Brainpool Program. The work of C.S.K. was supported in part by grant No. 1999-2-111-002-5 from the Interdisciplinary Research Program of the KOSEF, in part by the BSRI Program, Ministry of Education, Project No. 99-015-DI0032, and in part by Sughak program of Korean Reasearch Foundation, Project No. 1997-011-D00015. The work of K.L. was supported by the KOSEF 1998 Interdisciplinary Research Program.


text:
# On the late spectral types of cataclysmic variable secondaries

## 1 Introduction

Over last two decades, numerous studies have considered the evolution of the low–mass secondaries in cataclysmic variable (CV) systems. The calculations aimed mainly at reproducing the observed distribution of orbital periods $`P`$ of CVs, in particular the minimum period at 80 min and the dearth of systems in the 2-3 h period range (see e.g. King 1988 for a review). Based on the most popular explanation for this period gap, the disrupted magnetic braking model (cf. Rappaport, Verbunt and Joss 1983; Spruit and Ritter 1983), evolutionary sequences constructed with full stellar models or simplified bipolytrope models (cf. Hameury et al. 1988; Hameury 1991; Kolb and Ritter 1992; and references therein) reproduce the broad features of the observed $`P`$ distribution. Despite the significant uncertainties in the underlying stellar input physics (e.g. opacities, equation of state, etc…), the lack of a reliable description of angular momentum losses $`\dot{J}`$ that drive the mass transfer provided sufficient freedom that it was always possible to find a set of secondary star models which reasonably fit the period gap.

Input physics like atmospheric opacities and the outer boundary condition are essential for a comparison with observable quantities such as colours or spectral types.

The spectral type and the orbital period ($`P(R^3/M)^{1/2}`$; $`R`$ and $`M`$ are the secondary radius and mass) provide a unique set of constraints that test the structure both of the secondary’s interior and its outermost layers. As we shall show below, these can be used to find constraints on the actual orbital braking strength $`\dot{J}`$ as well.

Previously, an analysis of this kind was hampered by the lack of knowledge, both theoretical and observational, of stars in the low–mass range populated by CV secondaries. Here we make use of the marked and independent progress in the study of low–mass stars over the past decade. A wealth of ground-based and spaced-based observations has significantly improved our knowledge of the photometric and spectroscopic properties of stars at the bottom of the main sequence. In parallel, recent theoretical work has demonstrated the necessity to use accurate internal physics and outer boundary conditions based on non–grey atmosphere models to correctly describe the structure and evolution of low–mass objects (Burrows et al. 1993; Baraffe et al. 1995, 1997, 1998; Chabrier and Baraffe 1997, 2000, and references therein). These interior models, combined with recent, much improved models of cool atmospheres and synthetic spectra (see the review of Allard et al. 1997, and references therein) allow one to calculate self–consistent magnitudes and synthetic colours of low–mass stars. These can be compared directly to observed quantities, avoiding the use of uncertain empirical effective temperature $`T_{\mathrm{eff}}`$ and bolometric correction scales. The agreement between models and observations for (i) eclipsing binary systems (Chabrier and Baraffe 1995), (ii) mass–magnitude and mass–spectral type relations (Baraffe and Chabrier, 1996; Chabrier et al. 1996), (iii) colour–magnitude diagrams (Baraffe et al. 1997, 1998) and (iv) synthetic spectra (Leggett et al. 1996; Allard et al. 1997), demonstrate how reliable the theory of low–mass stars already is.

Recently, Beuermann et al. (1998) used our preliminary calculations based on these next–generation low–mass star models and compared calculated and observed spectral types (SpT) of CV secondaries. Beuermann et al. (1998) found that CVs with donors that are nuclearly evolved, with a substantial amount of central hydrogen depletion, can account for the late spectral types in long–period systems ($`P/\mathrm{h}6`$), but not in those close to the upper edge of the period gap ($`3P/\mathrm{h}6`$). Higher mass transfer rates as usually assumed can account for the latter.

The purpose of the present paper is to give a detailed account of the full set of our CV evolutionary calculations, a preliminary subset of which has been used by Beuermann et al. 1998. We discuss implications of the need for both high mass transfer rates and a significant fraction of systems with nuclearly evolved donors. We focus on CVs where the spectral type can be determined fairly accurately, i.e. on systems with orbital period above 3 h. In a previous paper (Kolb and Baraffe 1999), we considered aspects specific to the evolution of systems below the period gap ($`P2`$ h).

In §2 we summarize the input physics and present evolutionary sequences with initial ZAMS donors. Nuclearly evolved sequences are presented in §3. A discussion and conclusions follow in §4.

## 2 Evolutionary models with initial ZAMS donor

A brief description of the input physics for our stellar code can be found in Kolb and Baraffe (1999) and Baraffe et al. (1998). More details are given by Chabrier and Baraffe (1997) for the interior physics and by Hauschildt et al. (1999) for the atmosphere models. Here we just recall the main improvements compared to earlier models applied to CV donors: (i) the new models are based on the equation of state (EOS) of Saumon et al. (1995) which is specially devoted to low–mass stars and brown dwarfs, and to the description of strong non–ideal effects characterizing the interior stellar plasma; (ii) the outer boundary condition is based on non-grey atmosphere models (the use of a grey boundary condition overestimates the effective temperature $`T_{\mathrm{eff}}`$ and luminosity $`L`$ for a given mass $`M`$); (iii) mass–colour and mass–magnitude relations are derived self–consistently from the synthetic spectra of the same atmosphere models which provide the boundary condition.

We use the empirical SpT – $`(IK)`$ relation established by Beuermann et al. (1998) to determine the spectral type SpT for a given model ($`M`$, $`T_{\mathrm{eff}}`$, $`L`$) from the calculated colour $`IK`$. The accuracy and relevance of this conversion is discussed in Beuermann et al. (1998).

Because of the uncertainty of current descriptions for orbital angular momentum losses by magnetic braking and the substantial computer time (cf. Kolb and Ritter, 1990; Hameury 1991) required for evolutionary sequences if mass transfer is treated explicitly and self–consistently, we restrict our analysis to sequences calculated with constant mass transfer rate. We have checked that sequences with mass transfer driven by angular momentum losses $`\dot{J}`$ with either constant $`\dot{J}`$, constant $`\dot{J}/J`$, or $`\dot{J}`$ according to Verbunt & Zwaan (1981), have the same main properties as the sequences with constant $`\dot{M}`$ on which we base our conclusions.

In this section we consider sequences where mass transfer starts with donors of mass $`M_2=1\mathrm{M}_{}`$ that are initially on the ZAMS, with solar composition (X = 0.70, Z=0.02). We define a “standard sequence” which reproduces the width and location of the period gap (2.1 - 3.2 h) in the framework of the disrupted magnetic braking model. This sequence is calculated with $`\dot{M}=1.5\times 10^9M_{}\mathrm{yr}^1`$ until the donor becomes fully convective at $`M_2`$ = 0.21 $`\mathrm{M}_{}`$ and $`P`$ = 3.2 h. At this point mass transfer ceases and the secondary shrinks back to its thermal equilibrium configuration in $`23\times 10^8`$ years. Mass transfer resumes when the Roche lobe radius catches up with the donor radius at a period $`P`$ = 2.1 h. The transfer rate in this second phase is $`5\times 10^{11}M_{}\mathrm{yr}^1`$, a value typical for systems driven by gravitational wave emission (cf. Kolb and Baraffe 1999). This standard sequence is shown in Fig. 1 (thick solid line), together with sequences calculated with higher transfer rates ($`3\times 10^9M_{}\mathrm{yr}^1`$, $`10^8M_{}\mathrm{yr}^1`$ and $`10^7M_{}\mathrm{yr}^1`$).

Mass–losing stars deviate from thermal equilibrium. The deviation is large if the secondary’s thermal timescale $`t_{\mathrm{KH}}GM^2/RL`$ is long compared to the mass transfer timescale $`t_M=M/(\dot{M})`$. As shown in Fig. 2a, mass transfer causes a contraction of the donor radius (e.g. Whyte and Eggleton 1980; Stehle et al. 1996) relative to the corresponding ZAMS radius for predominantly radiative objects $`M0.6\mathrm{M}_{}`$, where the radiative core exceeds 80% of the total mass (cf. Chabrier and Baraffe 1997, their Fig. 9). Predominantly convective stars ($`M0.6\mathrm{M}_{}`$) expand relative to the ZAMS. Although a large mass transfer rate causes a significant departure of R from its equilibrium value (see Fig. 2a), Fig. 2b shows that $`T_{\mathrm{eff}}`$ is rather insensitive to $`\dot{M}`$ (see also Singer et al. 1993; Kolb et al. 2000).

The outer boundary condition of the stellar model is more important than $`\dot{M}`$ for the determination of $`T_{\mathrm{eff}}`$, as indicated by the dotted curve in Fig. 1 and Fig. 2b. This curve corresponds to a sequence with $`\dot{M}=1.5\times 10^9M_{}\mathrm{yr}^1`$, but calculated with the Eddington approximation, which implicitly assumes greyness and the diffusion approximation for radiative transfer in the atmosphere. Chabrier and Baraffe (1997) have shown that such an approximation overestimates $`T_{\mathrm{eff}}`$ for a given mass when molecules form in the atmosphere at $`T_{\mathrm{eff}}40004500`$ K. Note that the Eddington approximation sequence reproduces the period gap as well as the corresponding non-grey sequence, with a width and location of 2.2 - 3 h. However, the Eddington sequence consistently gives earlier spectral types than observed, and dramatically fails to reproduce the location of objects with $`P5`$ h in Fig. 1. This example highlights the significant improvements of our new generation of low–mass star models.

Note that the period $`P_{\mathrm{turn}}`$ where the period derivative changes from negative to positive due to departure from thermal equilibrium increases with mass transfer rate. Remarkably, the standard sequence reaches period bounce at $`P`$ = 3.2 h, the same period where the secondary star becomes fully convective. Sequences with higher mass transfer rate bounce before the donor becomes fully convective. In other words, for mass transfer rates $`\dot{M}1.5\times 10^9M_{}\mathrm{yr}^1`$, the minimum period reached by CVs just corresponds to or exceeds the observed upper edge of the period gap at $`3`$ h.

As mentioned in Beuermann et al. (1998), a large spread of the secular mean mass transfer rates from $`1.5\times 10^9M_{}\mathrm{yr}^1`$ to $`10^8M_{}\mathrm{yr}^1`$ could account for the observed range of late spectral type objects with $`P5`$ h. It is well known that such a spread of $`\dot{M}`$ in systems above the period gap cannot be excluded by observations (cf. e.g. Warner 1995, his Fig. 9.8; Sproats et al. 1996).

Although the sharp observed gap boundaries seem to imply uniform values of the secular mean $`\dot{M}`$ near the upper edge of the period gap (e.g. Ritter 1996), a spread of $`\dot{M}`$ cannot be dismissed unambiguously. The sharp upper edge is preserved if disrupted magnetic braking holds and $`\dot{M}`$ adopts only values $`1.5\times 10^9M_{}\mathrm{yr}^1`$. With increasing $`\dot{M}`$ the critical mass where the secondary becomes fully convective becomes smaller (cf. Fig. 2a). Thus, if the system detaches at this point, by analogy with the disrupted magnetic braking model, mass transfer resumes after the secondary has reached thermal equilibrium at a period shorter than 2 h. As an example, for the sequence with $`\dot{M}=10^8M_{}\mathrm{yr}^1`$ the donor becomes fully convective at $`M_20.15\mathrm{M}_{}`$, and the system reappears as CV at $`P1.7`$ h. Test calculations with a population synthesis code (e.g. Kolb 1993) show that such a spread of systems reappearing below the gap still gives a reasonably sharp lower edge of the period gap, with only a mild over–accumulation of systems in a period interval $`1.72.1`$ h. This is perfectly consistent with the observed distribution which does show a slight, although statistically probably not significant, accumulation of systems there.

Despite this consistency, accepting the $`\dot{M}`$ spread is not an attractive solution to the late SpT problem. It is hard to see why a presumably global angular momentum loss mechanism should drive such different transfer rates in otherwise similar systems.

Finally we note that there are no CV secondaries above the ZAMS line in the $`P`$SpT diagram for periods longer than $`6`$ h. This suggests that in these systems the secular mean $`\dot{M}`$ must be smaller than $`10^8M_{}\mathrm{yr}^1`$. But whatever the value of $`\dot{M}`$, sequences with unevolved donors never drop below the ZAMS line for $`P6`$ h.

## 3 Evolved sequences

We now consider sequences with initial donors which have evolved off the ZAMS prior to mass transfer. When mass transfer turns on, the star is old enough (cf. Table 1) to have burned a substantial amount of hydrogen in the centre. Only stars with mass $`M>0.35\mathrm{M}_{}`$, which develop a radiative core (cf. Chabrier and Baraffe 1997), start to deplete central H within a Hubble time. But only if $`M>0.8\mathrm{M}_{}`$ is more than $`50\%`$ of the initial H depleted in less than 10 Gyr. Therefore, this donor type is restricted to systems that form with a fairly massive secondary $`1\mathrm{M}_{}`$.

Figure 3 depicts sequences with such nuclearly evolved donors — hereafter “evolved sequences” — in the $`P`$–SpT diagram, for different central H mass fractions $`X_\mathrm{c}`$ at the time mass transfer starts. The initial secondary mass is $`M_2=1\mathrm{M}_{}`$ or $`M_2=1.2\mathrm{M}_{}`$ and mass transfer rate is fixed to $`1.5\times 10^9M_{}\mathrm{yr}^1`$. Table 2 gives for each of these sequences the secondary mass $`M`$, spectral type, $`T_{\mathrm{eff}}`$ and radius $`R`$ as a function of $`P`$, for $`3P/h10`$. The evolutionary track of ZAMS donor sequences in the $`P`$–SpT diagram is insensitive to the initial donor mass, while for evolved sequences with the same initial $`X_\mathrm{c}`$ $`T_{\mathrm{eff}}`$ decreases (the spectral type becomes later) at a given $`P`$ with increasing initial donor mass. This is illustrated in Fig. 3 where we plot sequences starting with $`M_2=1\mathrm{M}_{}`$ (long–dashed curve) and $`M_2=1.2\mathrm{M}_{}`$ (dash–dotted curve), both for $`X_\mathrm{c}=0.05`$.

As already emphasized by Beuermann et al. (1998), evolved sequences result in later spectral types for a given $`P5`$ h than the standard sequence and can explain amazingly well the spread of data observed above 6 h. The reason for this is illustrated in Figures 4 and 5, which display several diagnostic quantities along the sequences as a function of secondary mass (Fig. 4) and orbital period (Fig. 5).

When mass transfer turns on, more evolved donors with smaller $`X_\mathrm{c}`$ have larger $`R`$, $`T_{\mathrm{eff}}`$ and $`T_\mathrm{c}`$ than less evolved donors. This follows from the well known properties of the H burning phase of solar–type stars with a central radiative core: as hydrogen is depleted in the centre, the increase of the central molecular weight $`\mu `$ yields an increase of the central temperature and thus $`L`$, causing an expansion of the star. Hence the main–sequence evolution of a star with constant mass proceeds towards larger $`L`$ and $`T_{\mathrm{eff}}`$. The same effect exists for mass–losing main–sequence stars. Consequently, for a given mass, $`R`$ is larger for more evolved sequences, and $`P`$ is longer, while for a given $`P`$, the secondary mass and $`T_{\mathrm{eff}}`$ is smaller, hence the spectral type is later (cf. Fig. 3). As a consequence, the donor mass at any given $`P`$ is smallest for the most evolved sequence, cf. Table 2.

Although the evolved sequences shown in Fig 3 provide an attractive explanation for the observed data above 5-6 h, they are obviously in conflict with the standard explanation of the period gap. Indeed, with the adopted value $`\dot{M}=1.5\times 10^9M_{}\mathrm{yr}^1`$ we find that these sequences enter the period gap before the donor becomes fully convective or before the system bounces due to departure from thermal equilibrium. Table 1 gives the mass of the secondary and the mass $`M_{\mathrm{core}}`$ of its radiative core at $`P`$ = 3 h for the different sequences. Period bounce occurs for the ($`M_2=1\mathrm{M}_{}`$, $`X_\mathrm{c}=0.16`$) sequence at $`P`$ = 2.60 h, corresponding to $`M=0.1\mathrm{M}_{}`$ and $`M_{\mathrm{core}}=0.04\mathrm{M}_{}`$. The donor becomes fully convective only at $`M0.05\mathrm{M}_{}`$ ($`P`$ = 2.9 h). For more evolved sequences, period bounce would occur at even smaller periods.

This difference in the size of the radiative core is caused by the small central H abundance. The smaller $`X_\mathrm{c}`$ and the higher central temperature of evolved donors compared to unevolved donors (cf. Figs. 4 and 5) imply lower radiative opacities in the central region, and consequently favour radiative transport. As a result, the inward progression of the bottom of the donor’s convective envelope as mass decreases proceeds at a slower rate. Hence in evolved sequences the donor has a larger radiative core $`M_{\mathrm{core}}`$ than a donor with the same mass in an unevolved sequence (see Fig. 4). In terms of $`P`$, Fig. 5 shows that $`M_{\mathrm{core}}`$ is larger for the standard sequence at any given $`P3.5`$ h. The situation then reverses at the upper edge of the gap, which corresponds to $`M_{\mathrm{core}}=0`$ for the standard sequence. The effect of a lower central H abundance described above is similar to that reported by Laughlin et al. (1997) who followed the central H burning phase of very low–mass stars with $`M0.25\mathrm{M}_{}`$ (with constant mass). They note that stars which are fully convective on the ZAMS develop a radiative core towards the end of the main sequence because of the depletion of H and subsequent lowering of the central radiative opacities.

The $`P`$–SpT diagram suggests that at long orbital periods a fairly large fraction of systems has a significantly evolved donor. If these carried on evolving with a transfer rate $`1.5\times 10^9M_{}\mathrm{yr}^1`$ they would clearly overpopulate the period gap regime. This is made worse by the fact that some of these sequences reach their minimum period in the period gap (Fig. 3): at period bounce $`\dot{P}`$ = 0, and the detection probability $`d1/\dot{P}`$ increases sharply. We find a similar behaviour if mass transfer is treated self–consistently, with an angular momentum loss rate $`\dot{J}=`$ const. or $`\dot{J}J`$, calibrated such that the corresponding unevolved sequence reproduces the observed period gap (giving $`\dot{M}=1.5\times 10^9M_{}\mathrm{yr}^1`$ at $`P3`$ h).

An obvious way to prevent evolved sequences from evolving into the period gap region is to increase the mass transfer rate. This increases the deviation from thermal equilibrium and leads to period bounce above the upper edge of the observed gap. Test calculations with constant mass transfer rate show that period bounce occurs at $`P>3`$ h with $`\dot{M}3\times 10^9`$ and $`5\times 10^9`$ $`M_{}\mathrm{yr}^1`$ for sequences with respectively $`X_\mathrm{c}=0.16`$ and 0.05, for initial donor mass $`1\mathrm{M}_{}`$. The same behavior is found for donors with higher initial mass. In other words, for sequences with smaller $`X_\mathrm{c}`$ a higher $`\dot{M}`$ near the upper edge of the period gap is required than for less evolved donors, consistent with the fact that the thermal time is shorter if $`X_\mathrm{c}`$ is small. These evolved sequences with high transfer rates ($`\dot{M}>5\times 10^9M_{}\mathrm{yr}^1`$) also match the observed spectral types in the region $`P`$ = 3-6 h, just as high $`\dot{M}`$ sequences with ZAMS donors do (see §2 and Fig. 1).

Note however, that for long periods $`P`$ $`5`$ h a higher $`\dot{M}`$ has the same effect as for the unevolved sequences discussed in §2 (Fig. 1). The effective temperature is larger and therefore the spectral type earlier than for a sequence with smaller $`\dot{M}`$. Thus in order to account for the full range of rather cool spectral types at long periods, $`\dot{M}`$ cannot always be high along all evolved sequences.

To summarize: evolved sequences require high $`\dot{M}`$ ($`5\times 10^9`$) near the upper gap region to avoid crossing the gap, but a low mass transfer rate ($`\dot{M}5\times 10^9`$) for $`P6`$ h to explain the late spectral type CVs in this region. This suggests that the mean mass transfer rate increases during the secular evolution of nuclear evolved donors.

## 4 Discussion and conclusions

We performed calculations of the long–term evolution of CVs with up–to–date stellar models that fit observed properties of single low–mass stars exceptionally well. Our focus was on the orbital period– spectral type ($`P`$–SpT) diagram, and the observed location of CV secondaries which populate a band with spectral types significantly later than ZAMS stars at the same orbital period.

Our calculations tested if the observed spectral types can be reproduced by varying two main parameters, the mass transfer rate $`\dot{M}`$ (assumed constant), and the degree of nuclear evolution of the secondary before mass transfer starts (measured in terms of the initial central H fraction $`X_\mathrm{c}`$). Only CVs forming with massive ($`11.2\mathrm{M}_{}`$) donors, i.e. at long periods, can have donors where H is significantly depleted in the centre.

We summarize our results as follows: (i) In the framework of the standard discontinuous orbital braking model for the period gap, the observed gap between 2-3 h is well reproduced by sequences starting from ZAMS donors and proceeding at a mass transfer rate $`\dot{M}12\times 10^9M_{}\mathrm{yr}^1`$ near the upper edge of the gap. (ii) Higher transfer rates ($`\dot{M}5\times 10^9M_{}\mathrm{yr}^1`$) than assumed in (i) near the upper edge of the gap are required to explain the full range of late spectral types of secondaries in CVs with periods $`P=36`$ h. This is true whatever the evolutionary state of the donor at turn–on of mass transfer. (iii) A family of evolutionary sequences starting mass transfer from an evolved donor with varying $`X_\mathrm{c}<0.5`$ and low mass transfer rate ($`\dot{M}5\times 10^9M_{}\mathrm{yr}^1`$) covers the observed locations of CVs in the $`P`$SpT diagram for $`P6`$ h. (iv) For these sequences, higher mass transfer rates than in (i) are required near 3 h, otherwise they would evolve into the period gap and predict too early spectral types at shorter $`P`$.

If $`\dot{M}`$ is sufficiently large, the sequences reach their minimum period above the upper edge of the period gap. The donor becomes fully convective at a mass smaller than the canonical $`0.21\mathrm{M}_{}`$ of the standard (unevolved) sequence, but at a period longer than 3 h. If the system detaches at this point and the donor can re–establish thermal equilibrium, then mass transfer would resume at a period well below 2 h. This guarantees consistency with the standard discontinuous orbital braking model.

Points (ii) to (iv) show that if the mass transfer rate increases significantly along the secular evolution of evolved donors, these sequences can explain the observed scatter in the $`P`$SpT diagram both above 6 h, through the effect a lower $`X_\mathrm{c}`$ has on the evolutionary track (see §3 and Fig. 3), and below 6 h, through the effect of a higher $`\dot{M}`$ (see §2 and Fig. 1).

A large range of the secular mean mass transfer rate, as suggested by (ii), could be reconciled with the sharp boundaries of the observed period gap. In the standard model these are a result of the dominance of sequences with the canonical value $`12\times 10^9M_{}\mathrm{yr}^1`$ at periods close to 3 h. The apparent scatter of data in the $`P`$-SpT diagram does not support this dominance. Reasonably sharp boundaries could be maintained if the canonical value is the minimum $`\dot{M}`$ value of the range of transfer rates. The observed gap edges would then be defined by this minimum mass transfer rate sequence. Sequences with higher $`\dot{M}`$ would bounce before they reach the upper edge of the gap. They detach at longer periods and reattach at shorter periods than the standard sequence.

A more serious problem is that the large $`\dot{M}`$ range for otherwise similar system parameters, as suggested by (ii) for unevolved sequences, leaves the nature of the main control parameter determining the strength of orbital angular momentum losses completely undetermined. Our experiments show that $`X_\mathrm{c}`$ could be this parameter.

The full observed range of data in the $`P`$-SpT diagram could be explained if the orbital braking strength increases both with decreasing $`P`$ and decreasing $`X_\mathrm{c}`$. This ensures that period bounce prevents evolved sequences from overpopulating the period gap, while they still pass through rather late spectral types at long periods. The resulting $`\dot{J}`$ law must give the canonical value of $`\dot{M}`$ at 3 h for $`X_\mathrm{c}=0.7`$. A large $`\dot{M}`$ range for unevolved systems is not required.

With standard magnetic braking laws (e.g. by Verbunt & Zwaan 1981, Mestel & Spruit 1987) the transfer rate usually increases with period, and is smaller for sequences with more evolved donors (see e.g. Pylyser & Savonije 1989, Singer et al. 1993, Ritter 1994), i.e. has just the opposite differential behaviour than the one we propose. Given the significant uncertainty in these magnetic braking models and our poor knowledge of the underlying stellar magnetic dynamo, the suggestion that observations favour a non–standard law does not appear unrealistic. Further considerations of the theory of magnetic braking to assess our finding are clearly beyond the scope of this paper. We note, however, that there are magnetic braking scenarios in the literature that lead to an increase of $`\dot{M}`$ with progressing evolution. Kolb & Ritter (1992) obtained this for a Verbunt & Zwaan (1981)–type law when only the donor’s convective envelope is coupled to the orbital motion. Zangrilli et al. (1997) found this property for their boundary–layer dynamo.

At this point we re–emphasize that the input physics on which our stellar models are based provides an excellent description of observed properties of isolated low–mass stars. Unless there are effects peculiar to CV secondaries that are not taken into account by the models (e.g., effects due to the rapid rotation, or the irradiation by the white dwarf and accretion disc), one is forced to accept that a large fraction of CVs have a nuclearly evolved donor star. This is in stark contrast to predictions of standard models for the formation of CVs, where most CVs form with a secondary which is too low–mass ($`0.6\mathrm{M}_{}`$) to be evolved (Politano 1996, King et al. 1994). But this predominance of unevolved donors is mainly due to the neglect of systems forming with a secondary that is significantly more massive than the white dwarf (Ritter 2000). These undergo thermal–timescale mass transfer at a rate $`10^7M_{}\mathrm{yr}^1`$ and are usually associated with supersoft sources (e.g. di Stefano and Rappaport 1994). It is perfectly possible that these systems reappear as standard CVs once the mass ratio ($`q=`$ donor mass/WD mass) is sufficiently small. In descendants of this thermal–timescale evolution any degree of nuclear evolution prior to mass transfer is possible. De Kool (1992) finds that the impact of these survivor systems depends on the assumed initial mass ratio distribution which ultimately determines the secondary mass distribution in post–common envelope binaries. For de Kool’s model 1 (d$`N`$d$`q`$) survivor systems completely dominate the CV population. A full appraisal of the viability of the apparent predominance of nuclear evolved systems is beyond the scope of our study. This requires a full population synthesis with self–consistent treatment of $`X_\mathrm{c}`$–dependent angular momentum losses.

Finally, the accurate determination of donor masses could provide an observational test of our hypothesis that a significant fraction of CV donors is nuclearly evolved. Observational errors are still too large (e.g. Smith & Dhillon 1998) for this test. Systems with the most discrepant (late) spectral type should have the smallest mass at any given $`P`$.

### Acknowledgments

We are grateful to K. Beuermann, H. Ritter, H. Spruit, F. and E. Meyer for valuable discussions. We thank the referee, R. Smith, for a careful reading of the manuscript. I.B thanks the Max–Planck–Institut für Astrophysik (Garching) and the University of Leicester for hospitality during the realization of part of this work. The calculations were performed on the T3E at Centre d’Etudes Nucléaires de Grenoble. We thank A. Norton for improving the language of the manuscript.


text:
# Symmetry breaking and coarsening in spatially distributed evolutionary processes including sexual reproduction and disruptive selection

## Abstract

Sexual reproduction presents significant challenges to formal treatment of evolutionary processes. A starting point for systematic treatments of ecological and evolutionary phenomena has been provided by the gene centered view of evolution which assigns effective fitness to each allele instead of each organism. The gene centered view can be formalized as a dynamic mean field approximation applied to genes in reproduction / selection dynamics. We show that the gene centered view breaks down for symmetry breaking and pattern formation within a population; and show that spatial distributions of organisms with local mating neighborhoods in the presence of disruptive selection give rise to such symmetry breaking and pattern formation in the genetic composition of local populations. Global dynamics follows conventional coarsening of systems with nonconserved order parameters. The results have significant implications for the ecology of genetic diversity and species formation.

The dynamics of evolution can be studied by statistical models that reflect properties of general models of the statistical dynamics of interacting systems. Research on this topic can affect the conceptual foundations of evolutionary biology, and many applications in ecology, population biology, and conservation biology. Among the central problems is understanding the creation, persistence, and disappearance of genetic diversity. In this paper, we describe a model of sexual reproduction which illustrates mean field approaches (the gene-centered view of evolution) and the relevance of symmetry breaking and pattern formation in spatially distributed populations as an example of the breakdown of these approximations.

Pattern formation in genomic space has been of increasing interest in theoretical studies of sympatric speciation. These papers advance our understanding of the mechanisms of forming two species from one. However, they do not address the fundamental and practical problems of genetic diversity and spatial inhomogeneity within one species—a population whose evolution continues to be coupled by sexual reproduction. Moreover, and significantly, these papers do not address the implication of symmetry breaking and pattern formation for the gene centered view as a fundamental framework of evolutionary theory. In the following, we demonstrate that symmetry breaking and pattern formation invalidate the gene centered view (whether or not speciation occurs), and that they are important for the spatio-temporal behavior of the genetic composition of sexually reproducing populations. This has a wide range of implications for ecology, conservation biology, and evolutionary theory.

Before introducing the complications of sexual reproduction, we start with the simplest iterative model of exponential growth of asexually reproducing populations:

$$N_i(t+1)=\lambda _iN_i(t)$$
(1)
where $`N_i`$ is the population of type $`i`$ and $`\lambda _i`$ is their fitness. If the total population is considered to be normalized, the relevant dynamics is only of the proportion of each type, then we obtain

$$P_i(t+1)=\frac{\lambda _i}{_i\lambda _iP_i(t)}P_i(t)$$
(2)
where $`P_i`$ is the proportion of type $`i`$. The addition of mutations to the model, $`N_i(t+1)=_j\lambda _{ij}N_j(t)`$, gives rise to the quasi-species model which has attracted significant attention in the physics community. Recent research has focused on such questions as determining the rate of environmental change which can be followed by evolutionary change.

Sexual reproduction causes offspring to depend on the genetic makeup of two parents. This leads to conceptual problems (not just mathematical problems) in evolutionary theory because the offspring of an organism may be as different from the parent as organisms it is competing against. A partial solution to this problem is recognizing that it is sufficient for offspring traits to be correlated to parental traits for the principles of evolution to apply. However, the gene centered view is a simpler perspective in which the genes serve as indivisible units that are preserved from generation to generation. In effect, different versions of the gene, i.e. alleles, compete rather than organisms. This view simplifies the interplay of selection and heredity in sexually reproducing organisms.

We will show, formally, that the gene centered view corresponds to a mean field approximation. This clarifies the domain of its applicability and the conditions in which it should not be applied to understanding evolutionary processes in real biological systems. We will then describe the breakdown of the gene centered view in the case of symmetry breaking and pattern formation and its implications for the study of ecological systems.

It is helpful to explain the gene centered view using the “rowers analogy” introduced by Dawkins. In this analogy boats of mixed English- and German-speaking rowers are filled from a common rower pool. Boats compete in heats and it is assumed that a speed advantage exists for boats with more same-language rowers. The successful rowers are then returned to the rower pool for the next round. Over time, a predominantly and then totally same language rower pool will result. Thus, the selection of boats serves, in effect, to select rowers who therefore may be considered to be competing against each other. In order to make the competition between rowers precise, an effective fitness can be assigned to a rower. We will make explicit the rowers model (in the context of genes and sexual reproduction) and demonstrate the assignment of fitness to rowers (genes).

The rowers analogy can be directly realized by considering genes with selection in favor of a particular combination of alleles on genes. Specifically, for two genes, after selection, when allele $`A_1`$ appears in one gene, allele $`B_1`$ must appear on the second gene, and when allele $`A_1`$ appears on the first gene allele $`B_1`$ must appear on the second gene. We can write these high fitness organisms with the notation $`(1,1)`$ and $`(1,1)`$, and the organisms with lower fitness as $`(1,1)`$ and $`(1,1)`$. For simplicity, we assume below that the lower fitness organisms are non-reproducing. Models which allow them to reproduce, but with lower probabilities than the high fitness organisms, give similar results.

The assumption of placing rowers into the rower pool and taking them out at random is equivalent to assuming that there are no correlations in reproduction (i.e. no correlations in mate pairing) and that there is a sufficiently dense sampling of genomic combinations by the population (in this case only a few possibilities). Then the offspring genetic makeup can be written as a product of the probability of each allele in the parent population. This assumption describes a “panmictic population” which forms the core of the gene centered view often used in population biology. The assumption that the offspring genotype frequencies can be written as a product of the parent allele frequencies is a dynamic form of the usual mean field approximation neglect of correlations in interacting statistical systems. While the explicit dynamics of this system is not like the usual treatment of mean-field theory, e.g. in the Ising model, many of the implications are analogous.

In our case, the reproducing parents (either $`(1,1)`$ or $`(1,1)`$) must contain the same proportion of the correlated alleles ($`A_1`$ and $`B_1`$) so that $`p(t)`$ can represent the proportion of either $`A_1`$ or $`B_1`$ and $`1p(t)`$ can represent the proportion of either $`A_1`$ or $`B_1`$. The reproduction equation specifying the offspring (before selection) for the gene pool model are:

$`P_{1,1}(t+1)`$ $`=`$ $`p(t)^2`$ (3)
$`P_{1,1}(t+1)`$ $`=`$ $`P_{1,1}(t+1)=p(t)(1p(t))`$ (4)
$`P_{1,1}(t+1)`$ $`=`$ $`(1p(t))^2`$ (5)
where $`P_{1,1}`$ is the proportion of $`(1,1)`$ among the offspring, and similarly for the other cases.

The proportion of the alleles in generation $`t+1`$ is given by the selected organisms. Since the less fit organisms $`(1,1)`$ and $`(1,1)`$ do not survive this is given by $`p(t+1)=P_{1,1}^{}(t+1)+P_{1,1}^{}(t+1)=P_{1,1}^{}(t+1)`$, where primes indicate the proportion of the selected organisms. Thus

$$p(t+1)=\frac{P_{1,1}(t+1)}{P_{1,1}(t+1)+P_{1,1}(t+1)}$$
(6)
This gives the update equation:

$$p(t+1)=\frac{p(t)^2}{p(t)^2+(1p(t))^2}$$
(7)
There are two stable states of the population with all organisms $`(1,1)`$ or all organisms $`(1,1)`$. If we start with exactly 50% of each allele, then there is an unstable steady state. In every generation 50% of the organisms reproduce and 50% do not. Any small bias in the proportion of one or the other will cause there to be progressively more of one type over the other, and the population will eventually have only one set of alleles. This problem is reminiscent of an Ising ferromagnet at low temperature: A statistically biased initial condition leads to alignment.

This model can be reinterpreted by assigning a mean fitness (analogous to a mean field) to each allele as in Eq. (2). The fitness coefficient for allele $`A_1`$ or $`B_1`$ is $`\lambda _1=p(t)`$ with the corresponding $`\lambda _1=1\lambda _1`$. The assignment of a fitness to an allele reflects the gene centered view. The explicit dependence on the population composition (an Engligh-speaking rower in a predominantly English-speaking rower pool has higher fitness than one in a predominantly German-speaking rower pool) has been objected to on grounds of biological appropriateness . For our purposes, we recognize this dependence as the natural outcome of a mean field approximation.

We can describe more specifically the relationship between this picture and the mean field approximation by recognizing that the assumptions of no correlations in reproduction, a random mating pattern of parents, is the same as a long-range interaction in an Ising model. If there is a spatial distribution of organisms with mating correlated by spatial location and fluctuations so that the starting population has more of the alleles represented by 1 in one region and more of the alleles represented by $`1`$ in another region, then patches of organisms that have predominantly $`(1,1)`$ or $`(1,1)`$ form after several generations. This symmetry breaking, like in a ferromagnet, is the usual breakdown of the mean field approximation. Here, it creates correlations / patterns in the genetic makeup of the population. When correlations become significant then the species has two types, though they are still able to cross-mate and are doing so at the boundaries of the patches. Thus the gene centered view breaks down when multiple organism types form.

Understanding the spatial distribution of organism genotype is a central problem in ecology and conservation biology. The spatial patterns that can arise from spontaneous symmetry breaking through sexual reproduction, as implied by the analogy with other models, may be relevant. A systematic study of the relevance of symmetry breaking to ecological systems begins from a study of spatially distributed versions of the model just described. This model is a simplest model of disruptive selection, which corresponds to selection in favor of two genotypes whose hybrids are less viable. Assuming overlapping local reproduction neighborhoods, called demes, the relevant equations are:

$`p(x,t+1)`$ $`=`$ $`D(\overline{p}(x,t))`$ (8)
$`D(p)`$ $`=`$ $`{\displaystyle \frac{p^2}{p^2+(1p)^2}}`$ (9)
$`\overline{p}(x,t)`$ $`=`$ $`{\displaystyle \frac{1}{N_R}}{\displaystyle \underset{|x_j|R}{}}p(x+x_j,t)`$ (10)
$`N_R`$ $`=`$ $`\left|\{x_j||x_j|R\}\right|`$ (11)
where the organisms are distributed over a two-dimensional grid and the local genotype averaging is performed over a preselected range of grid cells around the central cell. Under these conditions the organisms locally tend to assume one or the other type. In contrast to conventional insights in ecology and population biology, there is no need for either complete separation of organisms or environmental variations to lead to spatially varying genotypes. However, because the organisms are not physically isolated from each other, the boundaries between neighboring domains will move, and the domains will follow conventional coarsening behavior for systems with non-conserved order parameters.

A simulation of this model starting from random initial conditions is shown in Fig. 1. This initial condition can arise when selection becomes disruptive after being non-disruptive due to environmental change. The formation of domains of the two different types that progressively coarsen over time can be seen. While the evolutionary dynamics describing the local process of organism selection is different, the spatial dynamics of domains is equivalent to the process of coarsening / pattern formation that occurs in many other systems such as an Ising model or similar cellular automata models. Fourier transformed power spectra (Figs. 24) confirm the correspondence to conventional coarsening by showing that the correlation length grows as $`t^{1/2}`$ after initial transients. In a finite sized system, it is possible for one type to completely eliminate the other type. However, the time scale over which this takes place is much longer than the results assuming complete reproductive mixing, i.e. the mean field approximation. Since flat boundaries do not move except by random perturbations, a non-uniform final state is possible. The addition of noise will cause slow relaxation of flat boundaries but they can also be trapped by quenched (frozen) inhomogeneity.

The results have significant implications for ecology of genetic diversity and species formation. The observation of harlequin distribution patterns of sister forms is generally attributed to nonhomogeneities in the environment, i.e. that these patterns reflect features of the underlying habitat (=selective) template. Our results show that disruptive selection can give rise to spontaneously self-organized patterns of spatial distribution that are independent of underlying habitat structure. At a particular time, the history of introduction of disruptive selection events would be apparent as a set of overlapping patterns of genetic diversity that exist on various spatial scales.

More specific relevance of these results to the theoretical understanding of genetic diversity can be seen in Fig. 5 where the population averaged time dependence of $`p`$ is shown. The gene centered view / mean field approximation predicts a rapid homogenization over the entire population. The persistence of diversity in simulations with symmetry breaking, as compared to its disappearance in mean field approximation, is significant. Implications for experimental tests and methods are also important. Symmetry breaking predicts that when population diversity is measured locally, rapid homogenization similar to the mean field prediction will apply, while when they are measured over areas significantly larger than the expected range of reproduction, extended persistence of diversity should be observed.

The divergence of population traits in space studied in our work can also couple to processes of speciation, i.e., processes that prevent interbreeding or doom the progeny of such breedings. These may include assortative mating, whereby organism traits inhibit interbreeding. Such divergences can potentially lead to the formation of multiple species from a single connected species (sympatric speciation). By contrast, allopatric speciation, where disconnected populations diverge, has traditionally been the more accepted process even though experimental observations suggest sympatric speciation is important.

Recent studies have begun to connect the process of symmetry breaking to sympatric speciation. Without considering pattern formation in physical space, we and other researchers have been investigating the role of pattern formation in genomic space as a mechanism or description of sympatric speciation. These studies include: a model of stochastic branching and fixation of subpopulations due to genetic drifts and local reproduction in genome space, general reaction-diffusion Turing pattern formation models in genomic space, and specific individual-based models of reproductive isolation involving assortative mating and disruptive selection (intrinsic disruptive selection, or disruptive selection arising from competition or sexual selection). Our work, presented here, is unique in discussing spatial inhomogeneity and genetic diversity within one species.

In conclusion, in formalizing sexual reproduction in evolutionary theory, we have found fundamental justification for rejecting the widespread application of the gene centered view. The formal mathematical analysis we presented to demonstrate the lack of applicability of the gene centered view is an essential step toward developing a sound conceptual foundation for evolution. We also showed that the gene centered view breaks down for species where local mating and disruptive selection give rise to symmetry breaking and pattern formation, which correspond to genetic inhomogeneity and trait divergence of subpopulations. The patterns formed undergo coarsening, following the usual universal spatio-temporal scaling behavior. The slow movement of boundaries between types cause long term persistence of genetic diversity through the local survival of (partially) incompatible types. This provides a new understanding of the development and persistence of spatio-temporal patterns of genetic diversity within a single species.

One should note that the context in which the gene centered view breaks down is of profound significance in applied aspects of modern ecology and conservation biology. The preservation of endangered species and ecosystems is currently at risk due to a dramatic decrease in their genetic diversity. We have described the implications of our results for the experimental observation of genetic diversity in endangered species. Our study of spatial patterns of genetic diversity in populations may also help guide the design of conservation areas and human directed breeding programs for endangered organisms.


text:
# Critical velocity in cylindrical Bose-Einstein condensates

## Abstract

We describe a dramatic decrease of the critical velocity in elongated cylindrical Bose-Einstein condensates which originates from the non-uniform character of the radial density profile. We discuss this mechanism with respect to recent measurements at MIT.

Superfluidity is one of the striking manifestations of quantum statistics, and the occurence of this phenomenon depends on the excitation spectrum of the quantum liquid. The key physical quantity is the critical velocity $`v_c`$, that is, the maximum velocity at which the flow of the liquid is still non-dissipative (superfluid). The well-known Landau criterion  gives the critical velocity as the minimum ratio of energy to momentum in the excitation spectrum:

$$v_c=\mathrm{min}\left(\frac{ϵ(k)}{k}\right)$$
(1)
(we put $`\mathrm{}=1`$ and the particle mass $`M=1`$). Actually, in liquid $`{}_{}{}^{4}\text{He}`$ the Landau criterion overestimates the critical velocity. The explanation of this fact was put forward by Feynman (see e.g. ) who suggested that superfluidity is destroyed by spontaneous creation of complex excitations (vortex lines, vortex rings, etc.). Extensive theoretical and experimental studies on this subject are reviewed in .

Bose-Einstein condensation of dilute trapped clouds of alkali atoms  offers new possibilities for the investigation of superfluidity . In the spatially homogeneous case, the spectrum of elementary excitations of a Bose-condensed gas is given by the Bogolyubov dispersion law :

$$ϵ(k)=\sqrt{\left(\frac{k^2}{2}\right)^2+2c_s^2\frac{k^2}{2}},$$
(2)
and the Landau critical velocity is equal to the speed of sound $`c_s=\sqrt{gn_0}`$ ($`n_0`$ is the condensate density, $`g=4\pi a`$, and $`a>0`$ is the scattering length). The first experimental observation of the critical velocity in trapped gaseous condensates, recently reported by the MIT group , gives a significantly smaller value of $`v_c`$.

The analyses in  and in recent theoretical publications (see e.g.  and references therein) employ the Feynman hypothesis and provide a qualitative explanation of the MIT experimental result. In this Letter we point out a simple geometrical effect which is characteristic for elongated cylindrical traps. Due to the non-uniform character of the radial density profile, the spectrum of axially propagating excitations in these traps is very different from the Bogolyubov dispersion law (2), and this difference leads to a strong decrease of the critical velocity. We show that this effect can at least partially explain the small critical velocity measured in the MIT experiment .

We consider an infinitely long cylindrical condensate which is harmonically trapped in the radial ($`\rho `$) direction. Then the condensate wave function $`\psi _0(\rho )`$ satisfies the Gross-Pitaevskii equation

$`\left({\displaystyle \frac{\mathrm{\Delta }_\rho }{2}}+V(\rho )\mu +g|\psi _0(\rho )|^2\right)\psi _0(\rho )=0,`$

where $`\mu `$ is the chemical potential, $`V(\rho )=\omega ^2\rho ^2/2`$ is the trapping potential, and $`\omega `$ the trap frequency. In the Thomas-Fermi regime, where the ratio $`\eta =\mu /\omega 1`$, the density profile is given by $`n_0(\rho )|\psi _0|^2=(\mu V(\rho ))/g`$ and the chemical potential is related to the maximum condensate density as $`\mu =n_{0max}g`$.

Elementary excitations can be regarded as quantized fluctuations of the condensate wavefunction . In our trapping geometry they are characterized by the axial ($`z`$) wavevector $`k`$ and radial angular momentum $`m`$. The corresponding part of the field operator reads

$`\delta \widehat{\psi }={\displaystyle \underset{m,k}{}}(u_{mk}\widehat{b}_{mk}v_{mk}^{}\widehat{b}_{mk}^{}),`$

where $`\widehat{b}_{mk}(\widehat{b}_{mk}^{})`$ are annihilation(creation) operators of the excitations. The excitation wave functions can be written in the form

$`(u,v)_{mk}(\rho ,z)=(u,v)_{mk}(\rho )\mathrm{exp}(im\varphi )\mathrm{exp}(ikz),`$

where $`\varphi `$ is the angle in the $`x,y`$ plane, and the radial functions $`u_{mk}(\rho )`$ and $`v_{mk}(\rho )`$ are solutions of the Bogolyubov-de Gennes equations (see e.g. )

$`ϵu`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Delta }_\rho }{2}}+{\displaystyle \frac{k^2}{2}}+V\mu +2n_0g\right)u+n_0gv,`$ (3)
$`ϵv`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Delta }_\rho }{2}}+{\displaystyle \frac{k^2}{2}}+V\mu +2n_0g\right)v+n_0gu.`$ (4)
Equations (3) and (4) constitute an eigenvalue problem. For given $`m`$ and $`k`$, they lead to a set of frequencies $`ϵ_{nm}(k)`$ characterized by the radial quantum number $`n`$ which takes integer values from zero to infinity. In the limit $`ϵ_{nm}(k)\mu `$, these modes were found for $`m=0`$ in the hydrodynamic approach in . For $`kR1`$, where $`R=(2\mu /\omega ^2)^{1/2}`$ is the Thomas-Fermi radial size of the condensate, the dispersion relation can be expanded in powers of $`k^2`$ :

$$ϵ_{n0}^2(k)=2\omega ^2n(n+1)+\frac{\omega ^2}{4}(kR)^2+O(k^4).$$
(5)
The lowest mode ($`n=0`$) represents axially propagating phonons: For this mode we have $`ϵ_{00}(k)=c_Zk`$, where the sound velocity $`c_Z=\sqrt{n_{0max}g/2}`$ is smaller by a factor of $`\sqrt{2}`$ than the Bogolyubov speed of sound at maximum condensate density, $`c_s`$ . The velocity $`c_s/\sqrt{2}`$ of axially propagating phonons has been measured in the MIT experiment .

The first correction to the linear behavior of the dispersion law $`ϵ_{00}(k)`$ at $`kR1`$ reveals its negative curvature: $`\delta ϵ_{00}=\omega (kR)^3/192`$ . In  the perturbative analysis leading to Eq.(5) was extended numerically to $`k1/R`$, still assuming that $`ϵ_{n0}(k)\mu `$. The calculations show that the group velocity of the first mode, $`dϵ_{00}/dk`$, monotonously decreases with increasing $`k`$ and can become significantly smaller than $`c_Z`$ characteristic for $`k1/R`$. This indicates that the critical velocity (1) associated with creating axial excitations ($`m=0`$, $`n=0`$) is smaller than $`c_s`$ and can also be reduced to below $`c_Z`$. The physical reason is that the decrease of the condensate density with increasing $`\rho `$ makes the axial superfluid flow less stable (see below).

The hydrodynamic approach used in  is not valid for $`ϵ\mu `$, where the excitation spectrum is no longer phonon-like and is dominated by the single particle dispersion relation $`ϵ(k)=k^2/2`$ (see Eq.(2)). The crossover between the two regimes occurs at $`kk_c=\mu ^{1/2}`$ and prevents the decrease of the group velocity with further increase in $`k`$. Obviously, the decrease of the critical velocity due to the radial inhomogeneity of the density profile can be dramatic only for $`k_c1/R`$, i.e. in large condensates with $`\eta 1`$.

The Bogolyubov-de Gennes equations (3) and (4) are valid for arbitrary $`k`$ and hence allow us to find the excitation spectrum in the crossover regime (see Fig.1) and establish the value of the critical velocity as a function of the Thomas-Fermi parameter $`\eta `$. Since the spectrum of excitations consists of a number of independent branches characterized by the quantum numbers $`n`$ and $`m`$, Eq.(1) gives a value $`v_c^{(nm)}`$ for the critical velocity corresponding to each mode. The results of our numerical calculations for the two lowest modes ($`n=0,m=0`$ and $`n=1,m=0`$) are presented in Fig.2. The breakdown of superfluidity occurs when the velocity of the flow matches the lowest of the velocities $`v_c^{(nm)}`$, which is proved to be $`v_c^{00}`$. The corresponding curve in Fig.2 indicates a significant decrease of the critical velocity compared to $`c_Z`$ already at $`\eta 10`$.

The decrease of the critical velocity with increasing ratio $`\mu /\omega `$ seems counterintuitive. One can increase this ratio by decreasing $`\omega `$ and keeping constant $`\mu `$. Then, at a constant density on the axis of the cylinder, for radially larger condensates one gets smaller $`v_c`$. However, this phenomenon has a clear physical explanation. The key point is the non-uniform character of the radial density profile. With increasing axial wavevector $`k`$, the wave functions of the excitations with $`kk_c`$ are more localized in the outer spatial region of the condensate and are thus more sensitive to the small value of $`n_0`$ in this region. It is this feature that provides a decrease of $`v_c`$ with increasing $`\eta `$, since for larger $`\eta `$ one has more possibilities to increase $`k`$ and still satisfy the condition $`kk_c`$. The described situation is not met in liquid helium where the density is practically constant and changes only in the region very close to the border of the sample.

We now turn to the discussion of the MIT experiment . The radial frequency in this experiment was $`\omega =2\pi \times 65`$ s<sup>-1</sup> and the chemical potential $`\mu =110`$ nK, so that $`\eta 35`$. The sound velocity $`c_s=6.2`$ mm/s was measured by observing a ballistic expansion. For these parameters we find $`v_c0.42c_s=2.6`$ mm/s, which is not far from the observed value $`v_c^{(exp)}1.6`$ mm/s . We do not pretend to explain the MIT data as the condensate in the experiment is not an infinite cylinder. The ratio of the radial to axial frequency for the MIT trap is $`3.3`$, which indicates that the effect of a finite axial size can be important. Nevertheless, our calculations may partially explain the small critical velocity observed in the experiment. Moreover, a smaller value of $`v_c`$ originating from our geometrical effect can facilitate the nucleation of vortex objects and thus provide a larger decrease of the critical velocity.

We would like to thank M. Lewenstein, A. Sanpera and A. Muryshev for fruitful discussions and numerical insights. The work was supported by the Austrian Science Foundation, by the Stichting voor Fundamenteel Onderzoek der Materie (FOM), by INTAS, and by the Russian Foundation for Basic Studies (grant 99-02-1802).