text:
# Testing intermediate-age stellar evolution models with VLT photometry of LMC clusters. II. Analysis with the Yale models1footnote 11footnote 1Based on observations collected at the European Southern Observatory, Chile (ESO No 64.L-0385).

## 1 Introduction

Stellar evolution theory is routinely tested using Galactic clusters, either globular (in general old, and both metal-poor and metal-rich) or open (intermediate-age or young, and generally more metal-rich than globular clusters). These clusters reflect the particular star formation and chemical enrichment history of the Milky Way, and certainly do not represent all possible stellar populations. For example, intermediate-age, metal-poor populations, which are not well represented in the Galaxy, are conspicuous in the Magellanic Clouds. This is an important part of the parameter space, which is relevant, for example, to the study of stellar populations of dwarf galaxies in the Local Group, or to that of galaxies at high redshifts, and therefore, in early evolutionary stages.

There are uncertainties that are critical to modeling intermediate-mass stellar models. Most notable of them are the treatments of convection (e.g., convective core overshoot) and mass loss. In this paper, we will address them, using intermediate-age Magellanic Cloud clusters as comparison templates.

Convective core overshoot (OS) is one of the physical parameters which must be taken into account during the core hydrogen burning phase of the evolution of intermediate-mass stars. Stellar evolution theory predicts that the CMD of young and intermediate-age star clusters is greatly affected by the amount of OS. In the early days of stellar structure theory, the convective core size was determined by the classic Schwarzschild (1906) criterion, ignoring the so-called OS, which is due to the inertial motion of materials beyond the formal core edge. Initially, the extent of core OS was thought to be negligible (Saslaw & Schwarzschild 1965). Subsequent theoretical studies, however, have emphasized the importance of convective core OS in stellar evolutionary models (Shaviv & Salpeter 1973).

Since then, many observational and theoretical studies have been devoted to investigating the effects of core OS (Prather & Demarque 1974; Bressan, Chiosi & Bertelli 1981; Bertelli, Bressan & Chiosi 1985; Stothers 1991; Meynet, Mermilliod & Maeder 1993). It is now generally believed that a moderate amount of core OS is necessary in order to reproduce the shape of the core contraction gap observed in CMDs of galactic open clusters (Maeder & Mermilliod 1981; Stothers & Chin 1991; Carraro et al. 1993; Daniel et al. 1994; Demarque, Sarajedini, & Guo 1994; Kozhurina-Platais et al. 1997; Rosvick & VandenBerg 1998), but the dependence of core OS on stellar mass is still poorly understood and the subject of on-going investigations (e.g. Ribas, Jordi, & Gimenez 2000).

The extent of OS is not currently derivable through direct theoretical calculations but can be determined empirically by comparing synthetic CMDs to observed ones, with OS in the models often being parameterized as a fraction of the pressure scale height at the formal edge of the convective core. In the CMDs of intermediate-age star clusters ($`2`$ Age $`5`$ Gyr), which have a distinctive gap near the turnoff (TO) caused by the rapid contraction of a hydrogen-exhausted convective core, the TO topology can be used to evaluate the extent of core OS (Demarque et al. 1994; Kozhurina-Platais et al. 1997). In the case of younger clusters, the importance of core OS can be estimated from the observed luminosity function of main sequence (MS) stars (Chiosi et al. 1989; Barmina, Girardi & Chiosi 2002) and from the ratio of main sequence to post-main sequence stars. Studies of detached eclipsing binaries can also provide empirical estimates of OS by direct comparison with stellar tracks in the log $`T_{\mathrm{eff}}`$ \- log $`g`$ plane (Schroder, Pols, & Eggleton 1997; Guinan et al. 2000; Ribas et al. 2000). However, the number of binary systems available for OS studies is very limited.

Many galactic open cluster studies have suggested a moderate amount of OS, i.e. $`20\%`$ of the local pressure scale height at the edge of a convective core. However, this has been tested only for near solar metallicity clusters in the Galaxy. So far, the presence of OS in metal-poor clusters in the Magellanic Clouds has only been tested in a limited sample of young (age $``$ 500 Myr) clusters (Lattanzio et al. 1991; Vallenari et al. 1991, 1994; Stothers & Chin 1992; Brocato, Castellani & Piersimoni 1994; Chiosi et al. 1995; Testa et al. 1999; Keller, Da Costa & Bessell 2001; Barmina et al. 2002). For Magellanic Clouds intermediate-age clusters, the necessary observations are more challenging due to the fainter magnitude of the TO, and HST or 8-m class telescopes under excellent seeing conditions are required.

The characteristics of mass loss during the RGB phase are still unclear in many aspects. Among the early papers on the subject, that of Reimers (1975) was most influential and his mass loss formula has been widely used to date:

$$\dot{M}=4\times 10^{13}\eta \frac{L}{gR}$$
(1)
where L, g, and R are luminosity, gravity and radius, respectively. The mass loss efficiency parameter $`\eta `$, which was later introduced, has been reported to vary from 0.25 to 2 – 3 for metal-rich stars Dupree (1986); Kudritzki & Reimers (1978); Renzini (1981). In contrast, the estimated range of $`\eta `$ for metal-poor stars appeared to be narrow, as most studies on metal-poor stars suggest $`\eta `$ = 0.3 – 0.7 Aaronson & Mould (1982); Mould & Aaronson (1982); Renzini (1981); Lee et al. (1994), and these values successfully reproduce the horizontal branch morphology of the old Galactic globular clusters Yi et al. (1999). While the application of Reimers’ formula for low mass stars predicts that the total amount of mass loss gets smaller as mass increases, it is not clearly demonstrated how the amount of mass loss depends on the mass of the stars in intermediate-age clusters (Tripicco, Dorman & Bell 1993; Liebert, Saffer & Green 1994; Carraro et al. 1996).

In this context, we started a project to investigate several aspects in the evolution of intermediate-age, metal-poor stellar populations, in particular, the amount of OS and the mass loss characteristics. To reach this goal, we observed three intermediate-age clusters using the Very Large Telescope (VLT) under excellent seeing conditions. The project and the observations are presented in Gallart et al. (2002; Paper I). In the current paper, we will compare the observations with the recent set of Yonsei-Yale ($`Y^2`$) stellar evolutionary models, and with models specifically calculated for this project using the same main input physics (but differing in the OS parameter), to constrain the amount of OS and the mass loss.

The three clusters under the study seem to present the signature of a core contraction gap in their TO region, which is predicted by the models as well for this age range. A possible obstacle in studies of OS using the shape of the core contraction gap is binary contamination. Castellani, Degl’Innocenti & Marconi (1999) pointed out that the contraction gap enhanced by core OS can be filled in by unresolved binary stars. In fact, some intermediate-age clusters, which are expected to have a core contraction gap, show a continuous sequence rather than a distinctive gap (Carraro et al. 1994; Paper I). Although this could be considered as a limiting factor for the topology method, we will show that it is possible, using synthetic CMDs, to put some constraints on the binary star population. It may be also expected that an increase of binary fraction in synthetic CMD simulations tends to decrease the value of OS required for a good fit. Testa et al. (1999) have claimed that a no-OS model with $`30\%`$ of unresolved binaries allows a fair fit to the observed MS luminosity function of NGC 1866 in the LMC, while Barmina et al. (2002), using a revised analysis of the same data, have concluded that the inclusion of unresolved binaries would not alter the need for OS (see also for other clusters, Vallenari et al. 1992; Carraro et al. 1994)

In this paper, we present a short description of the observational data (§ 2). Stellar models including OS treatment and synthetic CMD simulations are described in section 3. Detailed analysis on the steps of deriving ages, binary fraction, the amount of OS and mass loss is described in section 4, followed by a short discussion (§ 5) and conclusions (§ 6).

## 2 THE DATA

The three LMC clusters NGC 2173, SL 556 and NGC 2155 were observed with FORS1 at VLT-UT1, in service mode, through the filters $`V`$ and $`R`$, under excellent seeing conditions. PSF fitting photometry was obtained with the suite of programs DAOPHOTII/ALLSTAR/ALLFRAME (Stetson 1987, 1994), and calibrated using data of the same clusters obtained under photometric conditions using the MOSAIC camera at the CTIO-4m telescope. We estimate an accuracy of the zero points of the photometric transformation of $`\pm `$ 0.005. An assessment of the completeness and error characteristics of the data was performed through artificial star tests. Finally, we carefully subtracted the LMC field stars in the cluster field in a statistical fashion. All these aspects of the data treatment are discussed in detail in Paper I. In the following sections we will just use the photometrically calibrated, field-subtracted CMDs, and the results of the artificial star tests for completeness and errors.

## 3 MODELS

### 3.1 Stellar tracks & isochrones

We have constructed stellar tracks of metallicity $`Z=0.004`$ and $`0.007`$ with the Yale stellar evolution code YREC (Guenther et al. 1992), ignoring the effects of rotation. All input physics and model construction assumptions are the same as in the $`Y^2`$ Isochrones (Yi et al. 2001) except for the OS treatment. Three sets of tracks have been generated depending on the amount of core OS. The OS parameter is defined as a fraction of the pressure scale height at the edge of a convective core and prescribes the radial size of the OS region beyond the classical core convection boundary. We used 10, 20 and 30% of the local pressure scale height as the amount of OS, corresponding to “OS parameter” values of 0.1, 0.2 and 0.3, respectively.

From these tracks, we have constructed isochrones between 1 and 3 Gyr, with an interval in age of 0.1 Gyr. Color transformations are based on Lejeune, Cuisinier & Buser (1998), following the $`Y^2`$ scheme. It is generally found that a convective core develops when the stellar mass is larger than $`1.1`$$`M_{}`$, a result known to depend on chemical composition. For $`Y^2`$ isochrones, Yi et al. (2001) have adopted OS = 0.2 for young isochrones ($`<`$ 3 Gyr) and OS = 0.0 for old ones ($`3`$ Gyr) based on the observational studies of approximately solar abundance populations. Their isochrones have a 1 Gyr age grid for ages greater than 1 Gyr and do not provide the isochrones of ages greater than 2 Gyr but smaller than 3 Gyr. This is inconvenient to our study because our sample clusters appear to be in this age range. Besides, the $`Y^2`$ Isochrones group does not recommend any simplistic isochrone interpolation for the age range of 2 to 3 Gyr because it is in this age range that the transition between convective core and radiative core in the MS TO stars occurs (see also Woo & Demarque 2001). For this reason, we have computed a new $`Y^2`$ Isochrones in which the OS parameter has been kept constant for ages between 1 and 3 Gyr. Figure 1 illustrates the isochrones of two different OS parameters. Models with larger OS show brighter TOs and subgiant branches for each given age, implying that an age estimate for an observed CMD would be larger when using isochrones with larger OS.

### 3.2 Synthetic CMD code

Stellar mass is the fundamental parameter that determines the evolutionary path of each star in the CMD. By randomly distributing stellar mass with an initial mass function (IMF), and tracing each star’s evolutionary stage at a given time, one can construct synthetic CMDs for stellar systems with single age and metallicity. These simple stellar populations are relatively easy to model due to the simplifying assumption of instantaneous star formation rate and homogeneous chemical composition.

Several kinds of uncertainties still remain in computing synthetic CMDs for simple stellar populations. The first kind, which affects each star’s location on the CMD, comes from uncertainties in stellar evolutionary theory and stellar atmosphere models. Secondly, uncertainties in the number of stars in each location on the CMD come from the combinations of uncertainties in the evolutionary time scale of each stage, the slope of the IMF, and stochastic effects caused by small number statistics. Many CMD studies have aimed at testing physical parameters in stellar evolutionary theory by comparing theoretical shape and star counts with the observations while minimizing uncertainties from other sources. Depending on the nature of the targets and on the methods used in the analysis, these uncertainties can be minimized with rigorous effort. For example, plausible stochastic effects in number ratio analysis can be avoided by repeating the simulations and using large enough areas in the CMD to count stars in observed and synthetic CMDs.

A new CMD synthesis code has been created, including binary simulation, based on previous works (Woo 1998; Yi et al. 1999). For MS to RGB stars, Salpeter’s IMF with x = 1.35 has been used. For synthetic horizontal-branch construction, we have used the core helium-burning stellar models of Yi, Demarque & Kim (1997) with Lee, Demarque & Zinn (1990)’s prescription for mass distribution. Mass loss along the RGB is formally calculated with Reimers’ (1975) formula, and subtracted at the RGB tip to derive the mass of core helium burning stars. However, for this particular study, the amount of mass loss is empirically determined to reproduce the observed luminosity of the red clump.

Unresolved binary stars make a distinct sequence in the observed CMDs. It is critical to include binary stars for a detailed study of the CMD. We define the binary fraction (BF) as the number ratio of unresolved binary stars to the total number of stars in the CMD. The binary stars in the CMD are simulated in the following way. Once stellar mass is distributed according to the stellar IMF, some stars are randomly selected to be given a secondary star. The mass of the secondary star $`M_{se}`$ is assigned using the critical mass ratio, $`q_c`$ that sets the minimum mass ratio between secondary and primary:

$$M_{se}=(RAN(1q_c)+q_c)M_{pr}$$
(2)
where $`M_{pr}`$ is the primary star mass, and RAN is a random number with a flat distribution between 0 and 1. Such a flat distribution for the mass function of the secondary stars and a typical value of $`q_c`$=0.7 is adopted following Elson et al. (1998; see also Hurley & Tout 1999). Secondary stars with masses $`M_{se}0.7M_{pr}`$ would form a sequence indistinguishable from that of single stars. The BF that we will quote throughout the paper will be, therefore, the fraction of binary stars with $`M_{se}0.7M_{pr}`$, and a larger total BF should be expected.

Once the mass of the secondary star is determined, then the magnitudes $`M_V`$(pr) and $`M_V`$(se) of the primary and secondary, respectively, are combined to calculate the magnitude $`M_V`$(BI) for the unresolved binary system:

$$M_V(BI)=2.5log[10^{0.4M_V(pr)}+10^{0.4M_V(se)}]$$
(3)
### 3.3 Error & completeness simulation

The simulation of the observational errors (photometric shifts and completeness) in the synthetic CMDs has been performed on a star-by-star basis, using an empirical approach without any assumption on the nature of the errors or their propagation. The whole process has been performed in a way similar to that described in Aparicio & Gallart (1995) and Gallart et al. (1996a,b), and it is based on the results of the artificial star tests discussed in Paper I.

We have derived two tables from the artificial star tests. In the first one, the completeness information for each V, R magnitude interval is recorded. In the second table, for each artificial star that was recovered under the same quality criteria as the cluster data, the input and output magnitudes are listed. From this information, the simulation process for photometric errors and completeness in the synthetic CMD is the following: for each star in the synthetic CMD, with magnitudes $`V_s`$, $`R_s`$, we first use the completeness information in the corresponding magnitude interval to determine the probability of the star having been lost in the photometry, and keep it or remove it from the list of synthetic stars according to that probability. In a second step, for each synthetic star that has been kept, we select all artificial stars within a given distance in color and magnitude from the synthetic star. One of these stars is then picked up at random, and its $`\delta V`$ and $`\delta R`$, calculated as the differences between the recovered and injected magnitudes for that star, are added to the magnitudes $`V_s`$, $`R_s`$ in order to simulate the observational errors.

In Figure 2, synthetic CMDs before and after error simulations are compared as an example. In the left panel, single stars are placed on top of a given isochrone while the binary sequence has a dispersion according to the range of the mass ratio q ($`q_cq1`$). The error simulated CMD in the right panel shows many scattered stars in the CMD due to the photometric errors especially at faint magnitude levels. The result of the completeness simulation is that, globally, 10 to 15% of the total input stars are lost in the error simulated CMD.

## 4 Analysis: overshoot and mass loss

In this section, we compare the CMDs of the three LMC clusters, NGC 2173, SL 556 and NGC 2155, with isochrones and synthetic CMDs, in order to test the predictions of stellar evolution theory under different sets of assumed parameters regarding OS and mass loss. This will allow us not only to derive the cluster parameters (ages, metallicity, distance, and reddening) but also to constrain the input parameters in stellar model calculations (e.g., OS and mass loss).

### 4.1 The amount of overshoot in low-metallicity, intermediate-age populations

Cluster CMDs provide a powerful tool for deriving the OS parameter empirically. We have applied the following strategy on our sample clusters. First, we used the color difference $`\mathrm{\Delta }`$(V-R) between the base of the RGB and the TO to constrain the age range for each cluster under each OS assumption (Sarajedini & Demarque 1990). This has been complemented with visual isochrone fitting within that age range, to obtain a single best fitting age for each OS parameter. Second, since the presence of unresolved binary stars substantially affects the morphology of the TO and core contraction gap in the CMD, we attempt to find the BF using the ratio of the number of stars in two boxes around and below the core contraction gap. The last step takes advantage of the different spread in color near the core contraction gap depending on the assumed OS, to single out the best OS, and consequently, the best age for the cluster. Binary stars have substantial effects on the morphology of the TO and the core contraction gap, however, we will show that they do not greatly affect the color distribution of the stars, and therefore, the final conclusion on the best fitting OS parameter is almost independent of the assumed BF. These steps are explained in detail in the following three sections.

#### 4.1.1 Isochrone solutions

The major effect of core OS to stellar evolution is to supply hydrogen fuel from the OS-induced mixing region to the core, which lengthens the hydrogen core burning phase. This results in a brighter TO and subgiant branch since for a given age, more massive stars are going through the core contraction phase. Consequently, if the amount of OS assumed in the models is larger, an older age will be derived for a given observed CMD. Therefore, the OS parameter cannot be uniquely determined without an accurate age constraint and vice versa. Due to this degeneracy, we determined the best age for each OS parameter in the following steps.

As the TO gets redder with increasing age, the color difference, $`\mathrm{\Delta }`$(V-R) between the base of the RGB and the TO decreases. For a given age, a larger OS implies an increase of $`\mathrm{\Delta }`$(V-R) since more massive stars are located near the TO (Figure 3). Thus, for each OS parameter, an age can be estimated from the observed $`\mathrm{\Delta }`$(V-R) value. However, due to the uncertainties in determining the locations of the base of RGB and the TO in the observed CMD, this step only gives a plausible range of age for each OS. A typical error, 0.02 in $`\mathrm{\Delta }`$(V-R) gives age uncertainty $``$ 0.3 Gyr.

Within the age range, we have performed isochrone fitting to find the best fitting age. In this step, two further parameters are required: distance modulus and interstellar reddening. The distance modulus of the LMC is usually taken to be 18.5 $`\pm `$ 0.1, which is consistent with the recent compilation of Benedict et al. (2002). However, there is still some disagreement on the LMC distance modulus among experts in the field, and values ranging from 18.3 to 18.7 have been advanced. In addition, the LMC has a large angular distribution of clusters in the sky, among which the distance modulus may differ. For example, if we consider that the line of sight depth of the LMC is of the order of the LMC bar size, 3 degree ($``$ 2.6 kpc), the distance modulus can easily differ by up to 0.1 magnitude. Currently, the distance to individual clusters is unknown.

A mean reddening of E(B-V)$``$0.1 \[E(V-R)= 0.08\] or larger is tipically found for the LMC (Harris, Zaritsky, & Thompson 1997, and references therein). Most of these studies, however, focus on the central part of the LMC, or on regions containing young stars. Therefore, presumably, a larger amount of dust is present in these regions than in more external LMC regions, as is the case of the clusters in our study. Schlegel, Finkbeiner, & Davis (1998) calculate typical reddenings toward the LMC from the median dust emission in surrounding anuli, and find E(B-V)=0.075 (E(V-R)=0.06)), which should be close to the reddening found in the outer part of the LMC where internal extiction may be very low. No estimates of individual reddenings toward the clusters in our study are available.

Given these uncertainties, we performed isochrone fitting with limited freedom of distance modulus and reddening. We did not use a fiducial line obtained from the CMD in the fitting process since a fiducial line is difficult to determine around the convective core contraction gap, and binary stars can influence the fiducial line position. Since field star subtraction is not perfect, we rather put more weight on the isochrone shape near the TO and upper MS, which are relatively well populated.

With isochrone fitting, we single out the cluster age for each OS parameter, which is also complemented by comparing the observed CMD with the synthetic CMD simulated with each set of age and OS (Table 1). The true metallicities of our sample clusters are still poorly known but often suggested to be in the range Z=0.004 – 0.007 (Paper 1). We prefer to use a single value of metallicity for all three clusters in order to investigate systematic effects. Metallicity Z=0.007 seems to be inadequate as a representative value because the Z=0.007 isochrones do not yield reasonable matches (MS and RGB location, TO shape, and core contraction gap) to all three cluster CMDs due to the redder RGB. Therefore, we exclude the Z=0.007 solutions and only use Z=0.004 solutions for further analysis. However, it should be noted that NGC 2173 might be more metal-poor than Z=0.004 (Paper 1). For NGC 2155, OS=0.3 isochrones produce very poor matches to the observed CMDs and thus are excluded in the further analysis of this study. The final isochrone solutions are plotted with the observed CMD of each cluster (Figure 4). Note that the shapes of all three isochrones for each cluster are very similar except for the core contraction gap.

#### 4.1.2 Unresolved binary stars

i) The effect of binaries and overshoot on the shape of the gap

The effects of unresolved binary stars require special attention in cluster CMD studies since the binary sequence modifies the distribution of stars from the expected single star sequence. This contamination effect on the TO topology is more profound to intermediate-age clusters than to much younger or much older clusters since the core contraction gap is prominent only in intermediate-age clusters. We have performed CMD synthesis including binary stars in order to investigate how they affect the distribution of stars around the convective core contraction gap. For this test, we selected isochrone solutions of SL 556 with various BFs (0 to 30%).

Figure 5 shows how the shape of the CMD around the gap is affected by binaries and OS. For a given BF, larger OS results in a longer MS extending further toward the red due to the delayed core contraction phase. Consequently, the color distribution becomes much broader around the gap. Also, since the evolutionary time scale after core contraction phase is much faster with larger OS, the subgiant branch is less populous with increasing OS. For a given OS parameter, more stars are located around the gap with increasing BF. Binary stars make a rounded sequence, which resembles the shape of the single star sequence near the TO, which is $``$ 0.75 magnitude fainter. The stars above the gap are most likely binaries since single stars cannot be located above the gap even with much larger OS.

Binaries fill in the core contraction gap, which causes visual confusion near the TO and makes the derivations of age and OS parameter difficult. However, their sequence can still be differentiated from that of single stars using a detailed TO topology comparison. Furthermore, the number of stars above the core contraction gap, which are most likely unresolved binaries according to our analysis, can constrain the BF.

ii) Constraining the binary fraction

The binary frequency in a cluster CMD can be determined by artificial star tests using a binary sequence parallel to the MS (Rubenstein & Bailyn 1997). However, this method is time-consuming and still has a large uncertainty in estimating the BF. For our clusters, which are of intermediate-age and expected to display a core contraction gap, the distribution of stars around the core contraction gap can be utilized for estimating the BF.

We estimate the BF using the number ratio of the stars above the core contraction gap to the stars in the upper MS, respectively counted in two boxes in the CMD (Figure 6). The location of Box 1 is carefully selected to contain mostly binary stars and some subgiant stars (although for the oldest cluster NGC 2155 the number of subgiant stars increases). As mentioned in the previous section, we have compared a synthetic CMD without binaries to a CMD with 30% binaries, and determined the best location for the binary star box (Box 1). The location of Box 2 is selected to cover the TO region including a large fraction of MS stars and some binaries. If the total number of stars is fixed, a larger value of BF naturally results in an increase in number of stars in Box 1, while the number of stars in Box 2 is not significantly affected by the BF. Therefore, the number ratio of Box 1 to Box 2 can be a good indicator of the BF in the observed CMD.

The number ratio of stars in these two boxes for each cluster is shown in Figure 7 with a solid line. The uncertainties of the number ratios in the observed CMDs are estimated from the field star decontamination process and shown with dashed lines. In spite of the same box size, the observed number ratio of SL 556 is much greater than that of NGC 2173, which indicates a larger BF in SL 556. This result is consistent with the presence of a clear binary sequence parallel to the MS in the CMD of SL 556. In case of NGC 2155, a direct comparison with the other clusters is not possible since a Box 2 of much smaller size has been used due to the brighter limiting magnitude.

We have performed 10 Monte Carlo simulations for each given BF and averaged the number of stars in each box. Dots in Figure 7 show the mean number ratio of the two boxes in the synthetic CMDs, with error bars indicating 1 $`\sigma `$ dispersion. Simulations show a general trend of the binary effect on the number ratio, which increases with increasing BF as expected. According to our simulations, the number ratio varies slightly with different box locations and sizes, however, the trend appears to remain the same if we keep the same boxes in the observed CMD and the synthetic CMD. A comparison between the number ratios from the empirical data and from the simulations suggests 10-20% of unresolved binaries (with $`q0.7`$, the BF could be larger for any $`q`$) are required for all three clusters. These estimates are similar to those of Elson et al. (1998) for the LMC cluster NGC 1818 and of Rubenstein & Bailyn (1997) for the Galactic globular cluster NGC 6752.

#### 4.1.3 Overshoot from the color distribution around the core contraction gap

As mentioned in § 3.1, the empirical TO morphology may constrain the amount of OS and find the most probable isochrone solution in Table 1 and Figure 4. For a quantitative analysis, we have drawn Box 3 around the core contraction gap and analyzed the color distribution of stars in the box. The location of the box for each cluster, which contains most stars in core contraction phase, is shown in Figure 8.

For all synthetic CMD models with each isochrone solution (OS and age), the color distributions of stars in Box 3 are plotted in Figure 9. Each panel contains the distribution of stars from the observed CMD (histogram) and from simulations with various values of BF (10%-30%), which are suggested in the previous section. Note that once scaled to the total number of stars in Box 3, the color distribution is almost independent of BF, while OS affects the color distribution significantly. In general, the distribution is broadened and shifted to the red with increasing OS.

The best fits for NGC 2173 and SL 556 are models with OS=0.2, as clearly indicated by the median of the distribution. The various BF does not significantly affect the color distribution and, thus, OS estimation. The best OS value for NGC 2155, which is older than the other clusters, seems to lie between 0.1 and 0.2. We conclude that for all 3 LMC clusters the OS parameter is around 0.2. The fact that OS seems to be between 0.1 and 0.2 for NGC 2155 may indicate that OS slightly decreases with cluster age. A larger number of clusters in this age range would be necessary to confirm the trend.

### 4.2 Mass loss and the location of the red clump

The distinct feature of the red clump (RC) stars in the CMDs of intermediate-age clusters deserves considerable attention for its potential use as a distance indicator. However, it can be useful only when its characteristics (i.e., luminosity dependence on age and metallicity) are properly understood (Girardi & Salaris 2001; Grocholski & Sarajedini 2002).

An age dependence of the RC luminosity is naturally expected for the younger clusters with TO mass $`1.8M_{}`$, for which helium ignites in a non-degenerate core, since the core mass of helium burning stars increases as the age of a cluster is decreasing. In contrast, in the case of stars undergoing the He-flash, the degenerate core mass is nearly constant and independent of age, but the envelope mass can play a role in controlling the RC luminosity since the energy output from the hydrogen burning shell is decreasing with a smaller envelope mass. Thus, if mass loss during the RGB phase depends on age, there can be an age dependence in the RC luminosity of intermediate-age clusters.

In this section, we derive and discuss mass loss for each cluster using an empirically obtained relation between RC luminosity and RC mass. First, we have performed synthetic CMD simulations, using He-burning tracks (Yi et al. 1997), and measured the RC luminosity, $`M_V`$, as a function of the mean RC mass (Figure 10 top panel). It is shown that the median RC magnitude increases with larger RC mass, which is a generic feature of He-burning stellar tracks at this mass range (Vandenberg 1985; Seidel, Demarque & Weinberg 1987). Second, RC mass for each observed CMD is empirically derived from this luminosity–mass relation (points in Figure 10 top panel). The observed RC magnitude is converted to $`M_V`$ with the distance modulus, and a typical error (0.1 mag.) is assumed considering the uncertainties in distance modulus. It is interesting to note that NGC 2173 and SL 556 have a fainter RC magnitude than the older cluster NGC 2155, which is also shown in Figure 4, and consequently a smaller mean RC mass is inferred. Third, the amount of mass loss along the RGB has been estimated by subtracting the mean RC mass from the RGB tip mass of each isochrone. The mass loss errors are calculated from the assumed uncertainties in the RC magnitude.

We note that the estimated mass loss value for NGC 2155 is similar to that obtained using mass loss efficiency, $`\eta `$= 0.5 - 0.7 in Reimers’ formula, which successfully reproduces the horizontal branch morphology in old globular clusters Yi et al. (1999). However, the mass loss derived for SL 556 is too large to be compatible with $`\eta `$= 0.5 - 0.7 and may indicate a much higher mass loss efficiency in the context of Reimers’ formula. To illustrate this, we derived RC luminosity, $`M_V`$, for clusters with various age assuming a fixed mass loss value, 0.2 and 0.6 $`M_{}`$, respectively (Figure 10 bottom panel).

The origin of such a large estimate of mass loss in SSL 556 and disagreement with $`\eta `$=0.5-0.7 remains unclear to us and may have an important implication on the uncertainty/variability of the mass loss efficiency even in similar environments. The validity of mass loss estimates are subject to the uncertainties in distance modulus, age and metallicity although the color difference between the TO and the RGB is clear enough to indicate a reasonably accurate age estimate and the assumed metallicity provides a good match on the RGB slope. Further investigations are essential to draw more confident conclusions, especially with a better distance modulus constraint.

In the case of NGC 2173, it is unclear whether the RC stars experienced He-flash since the RGB tip mass of this cluster lies near the critical mass at which stars develop a degenerate core, although it is known that the critical mass becomes smaller with increasing OS (Chiosi et al. 1989). If not, then the lower RC luminosity of NGC 2173 may result from the fact that RC luminosity decreases as RGB tip mass becomes much larger than 1.5 $`M_{}`$ (see Girardi, Mermilliod & Carraro 2000). However, this cannot be the case for SL 556 and NGC 2155 since their RGB tip mass is similar or less than $`1.5`$ $`M_{}`$.

## 5 Discussion: Age gap among LMC clusters

Jensen, Mould & Reid (1988) first noticed a gap in the age distribution of LMC clusters, with none of them with ages between 4 and 10 Gyr. This age gap would be present also in metallicity, with younger LMC clusters having \[Fe/H\]$`0.7`$ and older ones having \[Fe/H\]$`2.0`$ (Olszewski et al. 1991). Geisler et al. (1997) performed a survey of candidate old clusters and confirm the existence of a gap in the LMC cluster formation, which they situate similarly between 3 and 8 Gyr ago. Shortly after, Sarajedini (1998) claimed to have found three 4 Gyr old LMC clusters, namely NGC 2155, SL 556 and NGC 2121, which would be in the age gap. With deeper observations, Rich, Shara & Zurek (2001) have estimated the age of these clusters to be 3.2 Gyr old, and claimed that they are probably the first clusters to have formed in the younger group of the LMC clusters. Our age estimate of NGC 2155 is slightly younger (by $`0.3`$ Gyr) but consistent with that of Rich et al. (2001), supporting the extension of the LMC age gap up to 3 Gyr on its intermediate-age side.

## 6 Conclusion

In this paper, we have performed detailed synthetic CMD simulations on the effect of OS and binary stars, compared with observed CMDs of three LMC clusters, NGC 2173, SL 556, and NGC 2155. Unresolved binary stars fill in the core contraction gap and make a unique sequence near the gap, which cannot be reproduced by single stars alone, even with a larger OS amount. We utilize the number ratios of stars around the core contraction gap to derive a first order estimate of the BF. Based on our star count analysis, all three clusters have a BF of 10% to 20% (with $`q0.7`$).

Using the color distribution of stars around the core contraction gap, which does not significantly depend on the BF, we find that the best OS parameter for each cluster is close to 20% of the local pressure scale height.

With the best solutions of the BF and OS from star count and color distribution analysis, we present synthetic CMDs for each cluster with the observed CMD in Figure 11, 12, and 13, respectively.

From the overall analysis, such as isochrone fitting, star counts, color distributions, and synthetic CMD comparison, we conclude that a moderate amount of OS, $`20\%`$ of the local pressure scale height, is essential to reproduce the observed shape around the core contraction gap in the CMDs of all three clusters, which are more metal-poor than the intermediate-age open clusters observed in the Galaxy, implying OS does not depend on metallicity at least for this metallicity range. With OS=0.2, our best age estimates are 1.5, 2.1 and 2.9 Gyr for NGC 2173, SL 556 and NGC 2155, respectively.

We constrainted mass loss along the RGB phase from the observed RC luminosity of each cluster, utilizing the RC luminosity – mass relation derived from synthetic CMDs. We found that the mas loss estimate for NGC 2155 is consistent with the typical mass loss parameterization that works for the Galactic globular clusters. However, the SL 556 data, and perhaps the data for the other young cluster NGC 2173 as well, indicate a substantially larger amount of mass loss than suggested by Reimers’ formula. This may indicate once again the complexity of mass loss processes.

This research is part of a Joint Project between Universidad de Chile and Yale University, partially funded by the Fundación Andes. Our data was collected as part of an ESO Service Mode run. C.G. acknowledges partial support from Chilean CONICYT through FONDECYT grant number 1990638. This research has been supported in part by NASA grant NAG5-8406 (P.D.). We thank H. Yong for his contribution in the early stage of this research.


format: text

text:
# Growth Algorithms for Lattice Heteropolymers at Low Temperatures

## I Introduction

Lattice polymers have been studied intensively to understand phenomena like the globule-coil transition of polymers, protein folding, etc. Protein folding (or, more precisely, protein fold prediction), one of the central problems of computational biology, refers to the determination of the ground state of protein molecules – which grosso modo is also its native state – from their amino acid sequence. Due to rapid advances in DNA analysis the number of known sequences has increased enormously, but progress in understanding their 3-dimensional structure and their functions has lagged behind owing to the difficulty of solving the folding problem.

Simplifying the description of a protein by replacing each amino acid by a simple point particle on a site of a regular lattice implies of course a great reduction of complexity, and one might wonder how much one can learn by this for real proteins. But even if this simplification is too strong, searching for lowest energy states of such models represents a paradigmatic example of combinatorial optimization. This will indeed be our main motivation: Finding algorithms that explore efficiently the low-energy states of a complicated energy landscape with many local minima. In addition to finding the ground state we want these algorithms also to sample excited states correctly, so that they provide a complete thermodynamic description – though we shall restrict ourselves in this paper to presenting results on ground states only.

A popular model used in these studies is the so-called HP model Dill85 ; Lau89  where only two types of monomers, H (hydrophobic) and P (polar) ones, are considered. Hydrophobic monomers tend to avoid water which they can only by mutually attracting themselves. The polymer is modeled as a self-avoiding chain on a regular (square or simple cubic) lattice with repulsive or attractive interactions between neighboring non-bonded monomers. Although also other interaction parameters have been used in the literature, almost all examples treated in this paper use energies $`ϵ_{HH}=1,ϵ_{HP}=ϵ_{PP}=0`$. The only other model studied here has also two types of monomers, for simplicity also called H and P (although they have identical hydrophobicities), but with $`ϵ_{HH}=ϵ_{PP}=1,ϵ_{HP}=0`$ tp92 . Chain lengths considered in the literature typically are between $`N=30`$ and $`N=100`$. Shorter chains do not present any problem, longer ones are too difficult.

A wide variety of computational strategies have been employed to simulate and analyze these models, including conventional (Metropolis) Monte Carlo schemes with various types of moves yue95 ; rrp97 ; deutsch97 , chain growth algorithms without tp92  and with re-sampling g98 ; bast98  (see also bd96 ), genetic algorithms um93 ; kd01 , parallel tempering irbaeck , and generalizations thereof chikenji ; chikenji2 , an ‘evolutionary Monte Carlo’ algorithm lw01 , and others tt96 . In addition, Yue and Dill yd93 ; yd95  also devised an exact branch-and-bound algorithm specific for HP sequences on cubic lattices which gives all low energy states by exact enumeration, and typically works for $`N\genfrac{}{}{0pt}{}{<}{}7080`$. If the chain is too long, it does not give wrong output but no output at all.

It is the purpose of the present letter to present two new variants of the Pruned-Enriched Rosenbluth Method (PERM) g97  and to apply them to lattice proteins. PERM is a biased chain growth algorithm with re-sampling (“population control”) with depth-first implementation. It has a certain resemblance to genetic algorithms, except that the latter are usually implemented breadth-first and do not allow to obtain correct Gibbs-Boltzmann statistics.

The original version of PERM was used for lattice protein folding in g98 ; bast98  and did extremely well. With one exception, it could find all known lowest energy configurations for all sequences tested in g98 ; bast98 , and found a number of new lowest energy states. The one case where it could not find the ground state in an unbiased and blind search was a 64-mer designed in um93  (see Fig. 1), but this is not surprising: Any chain growth algorithm should have problems in finding this configuration, since it has to grow a long arc which at first seems very unnatural and which is stabilized only much later. Indeed, at that time no other Monte Carlo method had been able to find this state either. But a very efficient algorithm, the Core-directed Growth method (CG) bd96  was overlooked in g98 ; bast98 . Thus PERM was not tested on the most difficult example known at that time, a 88-mer forming a $`\beta /\alpha `$-barrel whose ground state energy was known exactly. In the meantime, also two other improved Monte Carlo algorithms were published irbaeck ; chikenji . All this motivated us to take up the problem again.

## II The Algorithm

PERM is built on the old idea of Rosenbluth and Rosenbluth (RR) rr55  to use a biased growth algorithm for polymers, where the bias is corrected by means of giving a weight to each sample configuration. While the chain grows by adding monomers, this weight (which also includes the Boltzmann weight if the system is thermal) will fluctuate. PERM suppresses these fluctuations by pruning configurations with too low weight, and by “enriching” the sample with copies of high-weight configurations g97 . These copies are made while the chain is growing, and continue to grow independently of each other. PERM has been applied successfully to a wide class of problems, including e.g. the $`\mathrm{\Theta }`$ transition in homopolymers g97 , trapping of random walkers on absorbing lattices mg02 , and stretching collapsed polymers in a poor solvent gh02 . It can be viewed as a special realization of a “go with the winners” strategy aldous  which indeed dates back to the beginning of the Monte Carlo simulation era, when it was called “Russian roulette and splitting” kahn . Among statisticians, this approach is also known as sequential importance sampling (SIS) with re-sampling liu .

Pruning and enrichment were done in g97 ; g98 ; bast98  by choosing thresholds $`W_n^<`$ and $`W_n^>`$ depending on the estimate of the partition sums of $`n`$-monomer chains. These thresholds are continuously updated as the simulation progresses. If the current weight $`W_n`$ of an $`n`$-monomer chain is less than $`W_n^<`$, a random number $`r`$ is chosen uniformly in $`[0,1]`$. If $`r<1/2`$, the chain is discarded, otherwise it is kept and its weight is doubled. Thus low-weight chains are pruned with probability 1/2. Many alternatives to this simple choice are discussed in liu , but we found that more sophisticated strategies had little influence on the efficiency, and thus we kept the above in the present work. The determination of $`W_n^<`$ and $`W_n^>`$ will be discussed later. In principle we could use the same as in g98 ; bast98 , but we simplified it since the new variants are more robust, and some of the tricks employed in g98 ; bast98  are not needed.

On the contrary, we found that different strategies in biasing and, most of all, in enrichment had a big effect, and it is here the present variants differ from those in g98 ; bast98 . There, high-weight configurations were simply cloned (with the number of clones determined from the ratio of the actual weight to $`W_n^>`$), and the weight was uniformly shared between the clones. For relatively high temperatures this is very efficient g97 , since each clone has so many possibilities to continue that different clones very quickly become independent from each other. This is no longer the case for very low temperatures. There we found that clones often evolved in the same direction, since one continuation has a much higher Boltzmann weight than all others. Thus, cloning is no longer efficient in creating configurational diversity, which was the main reason why it was introduced.

The main modification made in the present paper is thus that we no longer make identical clones. Rather, when we have a configuration with $`n1`$ monomers, we first estimate a predicted weight $`W_n^{\mathrm{pred}}`$ for the next step, and we count the number $`k_{\mathrm{free}}`$ of free sites where the $`n`$-th monomer can be placed. If $`W_n^{\mathrm{pred}}>W_n^>`$ and $`k_{\mathrm{free}}>1`$, we choose $`2kk_{\mathrm{free}}`$ different sites among the free ones and continue with $`k`$ configurations which are forced to be different. Thus we avoid the loss of diversity which limited the success of old PERM. We tried several strategies for selecting $`k`$ which all gave similar results. Typically, we used $`k=\mathrm{min}\{k_{\mathrm{free}},W_n^{\mathrm{pred}}/W_n^>\}`$.

When selecting a $`k`$-tuple $`A=\{\alpha _1,\mathrm{}\alpha _k\}`$ of mutually different continuations $`\alpha _j`$ with probability $`p_A`$, the corresponding weights $`W_{n,\alpha _1}\mathrm{},W_{n,\alpha _k}`$ are (see Appendix)

$$W_{n,\alpha _j}=\frac{W_{n1}q_{\alpha _j}k_{\mathrm{free}}}{k\left(\genfrac{}{}{0pt}{}{k_{\mathrm{free}}}{k}\right)p_A}.$$
(1)
Here, the importance

$$q_{\alpha _j}=\mathrm{exp}(\beta E_{n,\alpha _j})$$
(2)
of choice $`\alpha _j`$ is the Boltzmann-Gibbs factor associated with the energy $`E_{n,\alpha _j}`$ of the newly placed monomer in the potential created by all previous monomers, and the terms in the denominator of Eq. (1) arise from correcting bias and normalization.

For the choice of continuations among the $`k_{\mathrm{free}}`$ candidates, we used two different strategies:

(1) In the first, called nPERMss for “new PERM with simple sampling”, we chose them randomly and uniformly, with the only restriction that they are mutually different. Accordingly, $`W_n^{\mathrm{pred}}=W_{n1}k_{\mathrm{free}}`$ footnote1 , and

$$W_{n,\alpha }=W_{n1}q_\alpha k_{\mathrm{free}}/k.$$
(3)
This has the advantage of simplicity, but it might at first appear to be inefficient. A priori, we would expect that some bias in favour of continuations with high Boltzmann weights or against continuations which run into dead ends might be necessary for efficiency.

(2) In the second, called nPERMis for “new PERM with importance sampling”, we did just that. For each possible placement $`\alpha [1,k_{\mathrm{free}}]`$ of the $`n`$-th monomer we calculated its energy $`E_{n,\alpha }`$ and its number $`k_{\mathrm{free}}^{(\alpha )}`$ of free neighbours, and used modified importances defined by

$$\stackrel{~}{q}_\alpha =(k_{\mathrm{free}}^{(\alpha )}+1/2)\mathrm{exp}(\beta E_{n,\alpha })$$
(4)
to choose among them. The predicted weight is now $`W_n^{\mathrm{pred}}=W_{n1}_\alpha \stackrel{~}{q}_\alpha `$. The replacement of $`q_\alpha `$ by $`\stackrel{~}{q}_\alpha `$ is made since we anticipate that continuations with less free neighbours will contribute less on the long run than continuations with more free neighbours. This is similar to “Markovian anticipation” cylinder  within the framework of old PERM, where a bias different from the short-sighted optimal importance sampling was found to be preferable.

The actual choice was made such that, for a given $`k`$ (remember that $`k`$ was already fixed by the ratio $`W_n^{\mathrm{pred}}/W_n^>`$), the variance of the weights $`W_n`$ is minimal. For $`k=1`$ this is standard importance sampling, $`p_\alpha =\stackrel{~}{q}_\alpha /_\alpha ^{}\stackrel{~}{q}_\alpha ^{}`$, and the variance of $`W_n`$ for fixed $`W_{n1}`$ would be zero if we had not replaced $`q_\alpha `$ by $`\stackrel{~}{q}_\alpha `$: $`W_{n,\alpha }=W_{n1}q_\alpha /p_\alpha =W_{n1}q_\alpha /\stackrel{~}{q}_\alpha _\alpha ^{}\stackrel{~}{q}_\alpha ^{}`$. For $`k>1`$, the probability to select a tuple $`A=\{\alpha _1,\mathrm{}\alpha _k\}`$ is found to be

$$p_A=\frac{\underset{\alpha A}{}\stackrel{~}{q}_\alpha }{_A^{}_{\alpha ^{}A^{}}\stackrel{~}{q}_\alpha ^{}}.$$
(5)
The corresponding weights are determined according to Eq. (1). The variance of the weight increase $`W_{n,\alpha }/W_{n1}`$, summed over all $`k`$ continuations within the tuple, would again be zero if $`q_\alpha `$ were not replaced by $`\stackrel{~}{q}_\alpha `$.

nPERMis is more time consuming than nPERMss, but it should also be more efficient. While Eq. (5) with $`\stackrel{~}{q}_\alpha `$ replaced by $`q_\alpha `$ would be optimal if the chain growth were a Markov process, it is not guaranteed to be so in the actual (non-Markovian) situation. We tried some alternatives for $`p_A`$, but none gave a clear improvement.

A noteworthy feature of both nPERMss and nPERMis is that they cross over to complete enumeration when $`W_n^<`$ and $`W_n^>`$ tend to zero. In this limit, all possible branches are followed and none is pruned as long as its weight is not strictly zero. In contrast to this, old PERM would have made exponentially many copies of the same configuration. This suggests already that we can be more lenient in choosing $`W_n^<`$ and $`W_n^>`$. For the first configuration hitting length $`n`$ we used $`W_n^<=0`$ and $`W_n^>=\mathrm{}`$, i.e. we neither pruned nor branched. For the following configurations we used $`W_n^>=CZ_n/Z_0(c_n/c_0)^2`$ and $`W_n^<=0.2W_n^>`$. Here, $`c_n`$ is the total number of configurations of length $`n`$ already created during the run, $`Z_n`$ is the partition sum estimated from these configurations, and $`C`$ is some positive number $`1`$. The following results were all obtained with $`C=1`$, though substantial speed-ups (up to a factor 2) could be obtained by choosing $`C`$ much smaller, typically as small as $`10^{15}`$ to $`10^{24}`$. The latter is easily understandable: with such small $`C`$, the algorithm performs essentially exact enumeration for short chains, giving thus maximal diversity, and becomes stochastic only later when following all possible configurations would become unfeasible. We do not quote the optimal results since they are obtained only for narrow ranges of $`C`$ which depend on the specific amino acid sequence, and finding them in each case would require an extensive search.

Since both nPERMss and nPERMis turned out to be much more efficient and robust than old PERM, we did not use special tricks employed in g98  like growing chains from the middle rather than one of the ends, or forbidding contacts between polar monomers.

In the following, when we quote numbers of ground state hits or CPU times between such hits, these are always independent hits. In PERM we work at a fixed temperature (no annealing), and successive “tours” g97  are independent except for the thresholds $`W_n^{<,>}`$ which use partially the same partition sum estimates. The actual numbers of (dependent) hits are much larger.

For both versions, results are less sensitive to the precise choice of temperature than they were for old PERM. As a rule, optimal results were obtained at somewhat lower temperatures, but in general all temperatures in the range $`0.25<T<0.35`$ gave good results for ground state search.

## III Results

(a) We first tested the ten 48-mers from yue95 . As with old PERM, we could reach lowest energy states for all of them, but within much shorter CPU times. As seen from Table 1, nPERMis did slightly better than nPERMss, and both were about one order of magnitude faster than old PERM. For all 10 chains we used the same temperature, $`\mathrm{exp}(1/T)=18`$, although we could have optimized CPU times by using different temperatures for each chain. In the following we quote in general only results for nPERMis, but results for nPERMss were nearly as good.

The CPU times for nPERMis in Table 1 are typically one order of magnitude smaller than those in bd96 , except for sequence #9 whose lowest energy was not hit in bd96 . Since in bd96  a SPARC 1 machine was used which is slower by a factor $`10`$ than the 167 MHz Sun ULTRA I used here, this means that our algorithms have comparable speeds.

(b) Next we studied the two 2-d HP-sequences of length $`N=100`$ of Ref. rrp97 . They were originally thought to have ground states fitting into a $`10\times 10`$ square with energies -44 and -46 rrp97 , but in bast98  configurations fitting into this square were found with lower energies, and moreover it was found that the configurations with lowest energies ($`E=47`$ resp. $`E=49`$) did not fit into this square. In the present work we studied only configurations of the latter type.

For the second of these sequences, new lowest energy configurations with $`E=50`$ were found later in chikenji , within 50 h CPU time on a 500 MHz DEC 21164A. We now hit this energy 7 times, with an average CPU time of 5.8 h on a 667 MHz DEC 21264 between any two hits.

For the first sequence of Ref. rrp97  we now hit several hundred times states with $`E=48`$, with ca. 2.6 min CPU time between successive hits. One of these configurations is shown in Fig. 2.

(c) Several 2-d HP-sequences were introduced in um93 , where the authors tried to fold them using a genetic algorithm. Except for the shortest chains they were not successful, but putative ground states for all of them were found in bast98 ; irbaeck ; chikenji . But for the longest of these chains (with $`N=64`$, see Fig. 1), the ground state energy $`E_{\mathrm{min}}=42`$ was found in bast98  only by means of special tricks which amount to non-blind search. With blind search, the lowest energy reached by PERM was -39. We should stress that PERM as used in bast98  was blind for all cases except this 64-mer, in contrast to wrong statements made in lw01 .

We now found putative ground states for all chains of um93  with blind search. For the 64-mer the average CPU time per hit was ca. 30h on the DEC 21264, which seems to be roughly comparable to the CPU times needed in chikenji ; irbaeck , but considerably slower than bd96 . As we already said in the introduction, this sequence is particularly difficult for any growth algorithm, and the fact that we now found it easily is particularly noteworthy.

On the other hand, nPERMis was much faster than bd96  for the sequence with $`N=60`$ of um93 . It needed $`10`$ seconds on the DEC 21264 to hit $`E_{\mathrm{min}}=36`$, and $`0.1`$ second to hit $`E=35`$. In contrast, $`E=36`$ was never hit in bd96 , while it took 97 minutes to hit $`E=35`$.

(d) A 85-mer 2-d HP sequence was given in konig , where it was claimed to have $`E_{\mathrm{min}}=52`$. Using a genetic algorithm, the authors could find only conformations with $`E47`$. In Ref. lw01 , using a newly developed evolutionary Monte Carlo (EMC) method, the authors found the putative ground state when assuming large parts of its known structure as constraints. This amounts of course to non-blind search. Without these constraints, the putative ground state was not hit in lw01  either, although the authors claimed their algorithm to be more efficient than all previous ones.

Both with nPERMss and with nPERMis we easily found states with $`E=52`$, but we also found many conformations with $`E=53`$. For nPERMis at $`\mathrm{exp}(1/T)=90`$ it took ca. 10 min CPU time between successive hits on the Sun ULTRA 1. One of those conformations is shown in Fig. 3.

(e) As two easy cases we studied the two longest sequences from kd01 , since we can compare there with CPU times given in kd01  for three versions of a supposedly very efficient genetic algorithm. These 2-d HP sequences with lengths $`N=33`$ and 48 have ground state energies -14 and -23, respectively. In kd01 , the most efficient version needed on average $`45`$min CPU (on an unspecified machine) to reach a ground state of the 33-mer. For the 48-mer only energy -22 could be reached, within $`2.5`$h per hit. Using $`\mathrm{exp}(1/T)=40`$, it took the Sun ULTRA 1 just 0.4 sec to hit one ground state of the 33-mer, 7 sec to hit $`E=22`$ for the 48-mer, and 16 min to hit a ground state of the 48-mer. Thus the present algorithm is roughly 1000 times faster than that of kd01 .

(f) Four 3D HP sequences with $`N=58`$, $`103`$, $`124`$, and $`136`$ were proposed in dfc93 ; lfd94  as models for actual proteins or protein fragments. Low energy states for these sequences were searched in tt96  using a newly developed and supposedly very efficient algorithm. The energies reached in tt96  were $`E=42`$, $`49`$, $`58`$ and $`65`$, respectively. With nPERMis, we now found lower energy states after only few minutes CPU time, for all four chains. For the longer ones, the true ground state energies are indeed much lower than those found in tt96 , see Table 2. Examples of the putative ground state configurations are shown in Figs. 4 to Fig. 7.

Note the very low temperatures needed to fold the very longest chains in an optimal time. If we would be interested in excited states, higher temperatures would be better. For instance, to find $`E=66`$ for the 136-mer (which is one unit below the lowest energy reached in tt96 ), it took just 2.7 seconds/hit on the DEC 21264 when using $`\mathrm{exp}(1/T)=40`$.

(g) Several 3-d HP sequences were studied in yd95 , where also their exact ground state energies were calculated using the ‘constrained hydrophobic core construction’ (CHCC) which is essentially an exact enumeration method tailored specifically to HP sequences on the cubic lattice. According to yd95 , CHCC can be used to find all exact ground state configurations for chains of length $`N70`$ to 88, depending on their degeneracies.

The longest chains given explicitly in yd95  together with their native configurations are a four helix bundle with $`N=64`$ and $`E_{\mathrm{min}}=56`$, and a chain with $`N=67`$ folding into a configuration resembling an $`\alpha /\beta `$ barrel with $`E_{\mathrm{min}}=56`$, too. Both have low degeneracy.

Finding ground states for the 64-mer was no problem for nPERMis. For $`\mathrm{exp}(1/T)=50`$, the DEC ALPHA 21264 machine needed on average 26.8 min CPU time to hit one of them. Things are a bit more interesting for the 67-mer. One of its ground states is shown in Fig. 8. Assume we want to let this grow, starting from the $`\beta `$ sheet end (monomer #67). Then we see that we always can form immediately stabilizing H-H bonds, and that we would be never seriously misled if we would place monomers greedily, at positions where they have low energies. Indeed, starting from this end we had no problems with nPERMis: It took on average 67 min to hit a native state on the DEC ALPHA 21264.

On the other hand, when starting with monomer #1, we were unsuccessful and the lowest energy reached was $`E=53`$, even after much longer CPU times. This is easily understood from Fig. 8: starting from this end we have to go repeatedly into directions which seem very unnatural at first sight, and which get stabilized much later.

Notice that the difference between the two growth directions is not that there is a folding nucleus when starting from #67, and no folding nucleus when starting from #1. After the first quarter is built up, both give the same $`\alpha /\beta `$ pair. Building this first quarter is no problem even when starting from #1, at least when we use $`C1`$ (in which case it is built essentially by complete enumeration). Thus the existence of a nucleus in the traditional sense is not sufficient. Instead it is crucial that further growth from this nucleus does not lead into false minima of the energy landscape.

(h) Next we studied the two-species 80-mer with interactions (-1,0,-1) that was introduced in tp92 . It was constructed in tp92  such as to fold into a four helix bundle with $`E=95`$, but two configurations with $`E=98`$ were found in g98  which essentially are $`\beta `$ sheet dominated. These configurations were hit on average once every 80 hours on a 167 MHz Sun ULTRA 1. Later they were also found in chikenji2 , with similar CPU time as far as we can tell. With nPERMis we needed only 5.3 hours / hit, on the same Sun ULTRA 1 (and for $`8\mathrm{exp}(1/T)12`$).

(i) Finally we also studied the 3-d HP sequence of length 88 given in bd96 . As shown there, it folds into an irregular $`\beta /\alpha `$-barrel with $`E_{\mathrm{min}}=72`$. This is the only chain whose ground state we could not find by our method, instead we only reached $`E=69`$. This is in contrast to the CG method which could find the lowest energy easily bd96 . The difficulties of PERM with this sequence are easily understood by looking at one of the ground states, see Fig. 9. The nucleus of the hydrophobic core is formed by amino acids #36-53. Before its formation, a growth algorithm starting at either end has to form very unstable and seemingly unnatural structures which are stabilized only by this nucleus, a situation similar to that in Fig. 1. In order to fold also this chain, we would have either to start from the middle of the chain (as done in bast98  for some sequences) or use some other heuristics which help formation of the hydrophobic core. Since we wanted our algorithm to be as general and “blind” as possible, we did not incorporate such tricks. The CG method, in contrast, is based on constructing an estimate of the hydrophobic core and the hydrophilic shell, and letting the chain grow to fill both in an optimal way, using a heuristic cost function.

Before leaving this section we should say that for all chains studied in this paper we found also states with $`E=E_{\mathrm{min}}+1,E_{\mathrm{min}}+2,\mathrm{}`$. Thus none of the sequences showed an energy gap above the (putative or exact) ground state. If such a gap is indeed typical for good folders, then none of the above sequences should be considered as good folders.

A list containing all sequences for which we found new lowest energy configurations is given in Table 3.

## IV Discussion

In the present paper we presented two new versions of PERM which is a depth-first implementation of the ‘go-with-the-winners’ strategy (or sequential importance sampling with re-sampling). The main improvement is that we now do not make identical clones of high weight (partial) configurations, but we branch such that each continuation is forced to be different. We do not expect this to have much influence for systems at high temperatures, but as we showed, it leads to substantial improvement at very low temperatures. The two versions differ in using simple sampling (nPERMss) resp. importance sampling (nPERMis) when choosing among possible branches.

Although the method could be used for a much wider range of applications (see chemnitz  for applications of PERM), we applied it here only to lattice heteropolymers with two types of monomers. These represent toy models of proteins, and we hope that our results will also foster applications to more realistic protein models. We showed only results for lowest energy configurations, but we should stress that PERM and its new variants are not only optimization algorithms. They also give information on the full thermodynamic behaviour. We skipped this here since finding ground states is the most difficult problem in general, and sampling excited states is easy compared to it.

Comparing our results to previous work, we see that we found the known lowest energy states in all cases but one. Moreover, whenever we could compare with previous CPU times, the comparison was favourable for our new algorithms, except for the CG method of Beutler and Dill bd96 . But we should stress that the latter is very specific to HP chains, uses strong heuristics regarding the formation of a hydrophobic core, and does not give correct Boltzmann weights for excited states. All that is not true for our method. In general nPERMis did slightly better than nPERMss, although the difference was much less than a priori expected.

In principle, essentially the same algorithms can also be used for off-lattice systems. This was already true for the original version of PERM which performed well for Lennard-Jones polymers at temperatures around the $`\mathrm{\Theta }`$-transition hg95 , but rather badly for collapsed heteropolymers at temperatures much below the $`\mathrm{\Theta }`$ temperature fg-unpub . Work is presently in progress to see whether the new versions of PERM perform better, and whether they can be used efficiently to study protein folding with realistic interactions.

Acknowledgements: We are indebted to Ralph Andrzejak for useful discussions and for critically reading the manuscript. WN was supported by the DFG (SFB 237).

## Appendix

In this appendix we shall collect some basic facts about random sampling, when tuples of instances are selected instead of individual instances. The discussion will be very general. On the other hand, we will not deal with problems specific to sequential sampling, i.e. we will assume that we sample only for the choice of a single item (e.g. for the placement of a single monomer).

Our central aim is thus to estimate a partition sum

$$Z=\underset{i=1}{\overset{N}{}}q_i$$
(6)
where the importances $`q_i`$ might e.g. be Boltzmann-Gibbs factors, and where $`N`$ is assumed to be finite (the generalization to infinite $`N`$ and to integrals instead of sums is straightforward). A conventional Monte Carlo (MC) procedure consists in choosing ‘instances’ $`i(\alpha ),\alpha =1,2,\mathrm{}`$ with probabilities $`p_{i(\alpha )}`$ such that each instance gives an unbiased estimate $`\widehat{Z}_1(\alpha )`$ (the index “1” will be explained in a minute). Thus, given $`M`$ such instances and letting $`M`$ tend to infinity, we have

$$Z=\underset{M\mathrm{}}{lim}\frac{1}{M}\underset{\alpha =1}{\overset{M}{}}\widehat{Z}_1(\alpha ).$$
(7)
One easily sees that

$$\widehat{Z}_1(\alpha )=\frac{q_{i(\alpha )}}{p_{i(\alpha )}}$$
(8)
does the job. Indeed,

$$\underset{M\mathrm{}}{lim}\frac{1}{M}\underset{\alpha =1}{\overset{M}{}}\frac{q_{i(\alpha )}}{p_{i(\alpha )}}=\underset{i=1}{\overset{N}{}}p_i\frac{q_i}{p_i}=Z.$$
(9)
At the same time we can also estimate the variance of $`\widehat{Z}_1`$. We have

$`\mathrm{Var}\widehat{Z}_1`$ $`=`$ $`\widehat{Z}_1^2\widehat{Z}_1^2`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}p_i\left({\displaystyle \frac{q_i}{p_i}}\right)^2Z^2`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{q_i^2}{p_i}}Z^2.`$
Up to now everything is correct for any choice of the probabilities $`p_i`$. They get fixed e.g. by $`p_i=1/N`$ (uniform sampling) or by demanding $`\mathrm{Var}\widehat{Z}_1`$ to be minimal, under the constraint $`_ip_i=1`$. This simple variational problem gives $`p_i^{\mathrm{opt}}q_i`$ which is known as importance sampling. For perfect importance sampling one finds furthermore that $`\mathrm{Var}\widehat{Z}_1=0`$.

Let us now assume that we select each time not one instance but $`K`$ instances, all of which are different. This requires of course $`KN`$. Moreover we will assume $`K<N`$, since otherwise this would amount to an exact summation of $`Z`$. An advantage of such a strategy should be that we obtain a more widely and uniformly spread sample. When $`NK`$, this should not have a big effect, but in our applications both $`N`$ and $`K`$ are small and the effect is substantial.

Thus each event consists in choosing a K-tuple $`\{i_1,i_2,\mathrm{},i_K\}`$, with the $`i_j`$ mutually different, from some probability distribution $`p_{i_1,i_2,\mathrm{},i_K}`$. We consider tuples related by permutations as identical, i.e. without loss of generality we can assume that $`i_1<i_2<\mathrm{}<i_K`$. Each choice $`\alpha `$ of a tuple $`\{i_1(\alpha ),i_2(\alpha ),\mathrm{},i_K(\alpha )\}`$ will lead to an estimate $`\widehat{Z}_K(\alpha )`$. Instead of Eq.(8) we have now

$$\widehat{Z}_K(\alpha )=\frac{N\underset{k=1}{\overset{K}{}}q_{i_k(\alpha )}}{K\left(\genfrac{}{}{0pt}{}{N}{K}\right)p_{i_1(\alpha )\mathrm{}i_K(\alpha )}},$$
(10)
since one verifies easily that $`\widehat{Z}_K(\alpha )=Z`$.

The variance of $`\widehat{Z}_K(\alpha )`$ is calculated just like that of $`\widehat{Z}_1(\alpha )`$,

$$\mathrm{Var}\widehat{Z}_K=\left(\genfrac{}{}{0pt}{}{N1}{K1}\right)^2\underset{i_1<\mathrm{}<i_K}{}\frac{(\underset{k=1}{\overset{K}{}}q_{i_k})^2}{p_{i_1\mathrm{}i_K}}Z^2.$$
(11)
Importance sampling is again obtained by minimizing it with respect to $`p_{i_1\mathrm{}i_K}`$, giving the result

$$p_{i_1\mathrm{}i_K}^{\mathrm{opt}}=\left(\genfrac{}{}{0pt}{}{N1}{K1}\right)\frac{\underset{k=1}{\overset{K}{}}q_{i_k}}{_{j=1}^Nq_j}.$$
(12)
The variance of $`\widehat{Z}_K`$ vanishes again for this choice.

On the other hand, for uniform (or “simple”) sampling, with

$$p_{i_1\mathrm{}i_K}^{\mathrm{ss}}=\left(\genfrac{}{}{0pt}{}{N}{K}\right)^1,$$
(13)
we obtain

$$\mathrm{Var}\widehat{Z}_K=\frac{(NK)N^2}{K(N1)}\mathrm{Var}q.(\mathrm{simple}\mathrm{sampling})$$
(14)
For $`K=1`$ this is the obvious result $`\mathrm{Var}\widehat{Z}_1=N^2\mathrm{Var}q`$, while for $`K=N`$ it gives $`\mathrm{Var}\widehat{Z}_N=0`$ as it should. For general $`1<K<N`$ the factor $`1/K`$ is trivial and results from the fact that each event corresponds to $`K`$ instances, while the factor $`(NK)/(N1)`$ gives the non-trivial improvement due to the fact that only different instances are chosen in each event.

Finally, when using Eq.(10) for sequential sampling, one has to attribute weights to each individual instance, instead of giving a weight only to the entire tuple. The obvious solution is

$$W_{i_k(\alpha )}=\frac{q_{i_k(\alpha )}N}{K\left(\genfrac{}{}{0pt}{}{N}{K}\right)p_{i_1(\alpha )\mathrm{}i_K(\alpha )}}.$$
(15)