text:
# On the relative frequencies of core-collapse supernovae sub-types: the role of progenitor metallicity

## 1 Introduction

It has been known for sometime now that supernovae with different spectral signatures and lightcurves, such as IIP, IIL, IIdw, Ib, Ic etc., all explode by the same mechanism: the gravitational energy released by the collapse of the Fe core of a massive star is partially transferred to the stellar envelope and, under certain conditions, manages to expell it to space. Despite more than fourty years of intense theoretical investigations, the conditions for a succesfull explosion remain unclear yet (e.g. Janka et al. 2002 for a recent overview).

However, the appearance of the supernova to an external observer, i.e. its optical spectra and lightcurve, depend only little on the mechanism of the explosion. The main factor affecting the observed features of core collapse supernovae is the mass (and the size) of the H-rich envelope of the progenitor star (see e.g. Hamuy 2003 for a recent overview). This depends on the amount of mass loss suffered during the star’s life. This is, in its turn, mainly a function of the stellar mass and metallicity. \[ Note: in the case of binary stars, the interaction of the progenitor with its companion may play an even more important role in determining the mass loss.\] It is expected then that the relative numbers of the various core collapse supernovae subtypes should depend on the metallicity of the host galaxy (e.g. Maeder 1992, Mowlavi et al. 1998).

In this work we use the well-known metallicity - luminosity relationship for late type galaxies (e.g. Garnett 2002 and references therein) to show that the observed ratio NIb,c/NII between supernova of type Ib,c and II does indeed depend strongly on metallicity (Sect. 3). To our knoweledge, it is the first time that such a relation is shown to hold observationally. However, we find that the slope of that relation, as well as the value of the ratio at solar metallicity, are higher than those expected from non-rotating models of single stars with mass loss. We show that predictions of recently calculated stellar models with rotation and mass loss (Meynet & Maeder 2003) offer a quantitative agreement with observations, at least for progenitor stars with solar initial metallicity. Finally, in Sect. 4 we argue that systematic surveys of supernovae at redshifts up to $`z`$1, feasible with the SNAP satellite, would allow to probe the role of metallicity on supernova properties during earlier epochs in the history of the Universe.

## 2 Types of core collapse supernovae

In a recent comprehensive study, Heger et al. (2003a) suggested the following scheme for the progenitors of core collapse supernova types, on the basis of models of single non-rotating stars with mass loss (Fig. 1a):

1) For progenitor metallicities Z$``$Z, stars with initial mass MIN$`<`$MIb,c$``$34 M keep enough of their envelope at the end of their lives to explode as SNII (either IIP or IIL). More massive stars lose completely their H-rich envelope, become WR stars and explode as SNIb,c.

2) For Z$`>`$Z MIb,c decreases, down to $``$30 M at Z$``$3 Z(since higher metallicities favour higher mass loss rates through radiation pressure on the envelope).

3) For Z$`>`$Z$`>`$0.1 Z, weaker mass loss rates increase MIb,c, which reaches $``$40 M at Z$``$0.1 Z.

4) For Z$`<`$0.1 Z, mass losses are negligible and the upper mass limit for SNII is kept at 40 M (Fig. 1). That upper limit no more coincides with MIb,c, since stars in the mass range MIb,c$`>`$MIN$`>`$40 M and in that metallicity range fail to explode (at least under “standard” assumptions of current models, i.e. spherical symmetry) and collapse directly to form black holes111 According to Heger et al. (2003a) non-”normal” supernovae, powered by jets, may occur in that mass-metallicity region; but their optical signature differs from both SNII or SNIb,c.. Only the most massive stars (see dotted curve in upper panel of Fig. 1) produce a weak SNIb,c explosion in that metallicity range.

A direct test of that scheme was recently undertaken by Smartt et al. (2003). By using high quality images of nearby host galaxies of SN prior to the explosion, and by comparing to stellar evolution tracks with mass loss but no rotation, they were able to evaluate (or put upper limits on) the mass of the progenitor stars of eight core collapse supernovae. Seven of them were SNII (either P or L) and were found to have MIN$`<`$25 M, in agreement with theoretical expectations. The last one (SN2002ap), classified as Ic, was found to have MIN$`<`$40 M and Z$``$0.5 Z, in (marginal) agreement with the scheme proposed by Heger et al. (2003a).

In this work we propose an indirect test of the ideas of Heger et al. (2003a), concerning the number ratio of SNIb,c to SNII core collapse supernovae as a function of metallicity. For metallicities Z$`>`$0.1 Z, where the upper mass limit for SNII coincides with MIb,c this ratio is given by

$$\frac{N_{\mathrm{Ib},\mathrm{c}}}{N_{\mathrm{II}}}=\frac{_{M_{\mathrm{Ib},\mathrm{c}}}^{M_{\mathrm{SUP}}}\mathrm{\Phi }(M)𝑑M}{_{M_{\mathrm{INF}}}^{M_{\mathrm{Ib},\mathrm{c}}}\mathrm{\Phi }(M)𝑑M}$$
(1)
where $`M_{\mathrm{INF}}`$9 M is the lower mass limit for core collapse supernova (Heger et al. 2003a), $`M_{\mathrm{SUP}}`$100 M is the corresponding upper mass limit (on which the value of NIb,c/NII depends very little) and $`\mathrm{\Phi }(M)`$ is the stellar initial mass function. Eq. (1) is obviously correct for a “steady state” situation, i.e. for a star formation rate $``$constant during at least the last $``$25 Myr (the lifetime of a star with mass $`M_{\mathrm{INF}}`$); it does not apply e.g. in the case of a starburst more recent than 10 or 20 Myr.

For metallicities Z$`>`$0.1 Z, one sees that the ratio NIb,c/NII increases monotonically with Z, since the dividing mass MIb,c decreases with increasing Z. This is true as far as the stellar initial mass function (IMF) remains constant with metallicity, which appears to be the case according to most current studies (e.g. Kroupa 2002). Of course, the value of NIb,c/NII depends on the form of the adopted IMF in the massive star range, as can be seen on Fig. 1b. The two thick curves are constructed, respectively, with the Salpeter IMF (a single power law with slope X=1.35) and the Kroupa et al. (1993) IMF (a multi-slope power law with X=1.7 in the massive star range) 222The power-law stellar IMF of slope X is defined as $`dN/dMM^{(1+X)}`$.. In both cases we assume that the IMF extends to $`M_{SUP}`$=100 M. Indeed, although stars with MIN$`>`$100 M might have existed in $``$zero metallicity environments (Nakamura & Umemura 2001), their existence at non-zero metallicities is implausible, both on theoretical (e.g. Baraffe et al. 2002) and observational grounds (e.g. Heydari-Malayeri 2003). We assume then in the following that for Z$`>`$0.1 Z the IMF extends only up to $``$100 M.

An inspection of Fig. 1b shows that NIb,c/NII never exceeds $``$0.2, even at the highest metallicities calculated by Heger et al. (2003a) with their non-rotating models. This is in disagreement with observations. Indeed, Bressan et al. (2002) find that in normal galaxies NIb,c/NII$``$0.27, while Hamuy (2003) argues that one out of four core collapse supernovae belongs to type Ib,c, i.e. that the ratio NIb,c/NII should be $``$0.33.

In the next section we show that, by counting supernova types (and corresponding ratios) in galaxies of different metallicities, the prediction of Fig. 1b can be tested statistically. We find that the predictions of non-rotating models disagree not only with observed values at high metallicities, but also with the observed slope of the NIb,c/NII vs metallicity relation. We suggest that models including rotation offer a physically plausible quantitative solution.

## 3 The NIb,c/NII ratio in galaxies

The existence of a metallicity-luminosity relation is a well established fact about spiral galaxies (see Garnett 2002 and references therein). Oxygen abundances measured in the gaseous phase and at fixed fractional galactocentric radius (at the center, at one disk scalelength or at the disk half-radius) always increase with the blue luminosity of the galaxy. According to Garnett (2002) this relation can be expressed as

$$log\left(\frac{O}{H}\right)=0.16\mathrm{M}_\mathrm{B}6.4$$
(2)
where MB is the galaxy’s blue magnitude and the metallicity is taken at the galaxy’s effective radius Reff. The relation is valid over $``$2 orders of magnitude for O/H and over a range of $``$10 magnitudes, for spirals, irregulars and low surface brightness galaxies alike (see Fig. 2, top).

This important feature has not found a satisfactory explanation up to now. It could be due either to a greater efficiency of star formation in more massive galaxies (e.g. Ferreras and Silk 2001), to a younger age for less massive galaxies (e.g. Boissier et al. 2001), to a more important mass loss reducing the effective yield in less massive galaxies (e.g. Dekel and Silk 1986) or to some combination of these factors. Independently of its origin, this relation allows us to test statistically the theoretically predicted trend of NIb,c/NII with metallicity.

We use a recent version (january 2003) of the Asiago Supernova Catalogue, presented in Barbon et al. (1999)333Available at http://merlino.pd.astro.it/~supern/.. This catalogue provides the morphological type of the galaxy, the supernovae type and the blue magnitude of the parent galaxy when possible (usually from the RC3 or the LEDA database). From the entries of the catalogue, we keep only the galaxies of spiral or irregular morphologies, (for which the metallicity-luminosity relationship applies) and which show no signs of recent starbursts (so that Eq. (1) applies); we also take into account only the core collapse supernovae, i.e. those with a clear identification as one of the types Ib, Ib/c, Ic, or II. The parent galaxy’s magnitude is available for 280 of the core-collapse events. For that restricted sample, we find a value of NIb,c/NII=0.27, i.e. similar to the value given by Bressan et al. (2002). To study the effect of metallicity on the number of SN of each sub-type, we adopt magnitude bins containing each the same number of core collapse events (50). In the following we adopt errorbars $`\delta N=\sqrt{N}`$ for any number N, and the relative uncertainty of a ratio $`N/M`$ is $`\delta N/N+\delta M/M`$.

In Fig. 2 (middle) we plot the number of the core collapse SN sub-types in each bin. Despite the poor statistics, it can be seen that the total number of SNII increases in the low luminosity bins, while the one of SNIb,c decreases. This is translated into a clear increase of the NIb,c/NII ratio with galactic luminosity, presented in Fig. 2 (bottom).

The obvious physical explanation for the observed corelation between NIb,c/NII and MB involves the effect of metallicity, as explained in Sect. 2: more luminous and metal-rich galaxies are expected to have higher ratio of NIb,c/NII. Note that the data span the metallicity range 1/3 Z to 2 Z, i.e. the range where theoretical models predict indeed a variation of NIb,c/NII with metallicity.

However, before discussing further the metallicity dependence of NIb,c/NII, one should account for the fact that Eq. (2) is valid at a fixed radius Reff, while galaxies display important metallicity gradients and, in consequence, supernova inside a galaxy of a given magnitude are expected to occur in regions of different metallicities; this introduces, in principle, a dispersion of the NIb,c/NII ratio around some average value corresponding to the mean galaxian metallicity. The question is then whether the observed supernovae exploded in regions with metallicity close to the one given by Eq. (2), i.e. close to Reff. An inspection of Fig. 1 in Bressan et al. (2002) shows that this is indeed the case: the distribution of observed core collapse supernovae as a function of their galactocentric distance peaks at 0.6 R25 (R25 is the radius where the B band surface brightness is 25 mag arcsec-2); on the other hand, for exponential disks with “canonical” central surface brightness ($`\mu _0`$=21.65 mag/arcsec, the Freeman value), Reff corresponds to $``$0.56 R25444 For an exponential disk with a scalelength $`R_\mathrm{D}`$, the surface brightness profiles has the form: $`\mu (R)=\mu _0+1.085R/R_\mathrm{D}`$. For $`\mu _0`$=21.65, $`\mu `$=25 occurs at $`R=R_{25}3R_\mathrm{D}`$. The radius within which the integral of this profile is equal to half the integral up to an infinite radius is the half-light radius $`R_{\mathrm{eff}}`$, and is equal to 1.865 $`R_\mathrm{D}`$. Then, $`R_{\mathrm{eff}}/R_{25}`$ 0.6.. Thus, for most of the observed core collapse supernovae, the magnitude of the host galaxies can indeed be associated to the metallicities implied by Eq. (2).

A best fit to the observed NIb,c/NII vs metallicity relation is given by

$$\frac{N_{\mathrm{Ib},\mathrm{c}}}{N_{\mathrm{II}}}=0.35log(Z/Z_{})+0.34$$
(3)
The theoretical NIb,c/NII ratio (i.e the curves of Fig. 1b) are also plotted in Fig. 2 (bottom), by assuming the luminosity-metallicity relation of Eq. (2). It can be seen that theory matches observations only at the bin of lowest luminosity, corresponding to a metallicity of $``$0.3 Z. At higher metallicities/luminosities there is an increasing difference between theory and data, by a factor of 2-3 at Z and up to 3.5 at $``$2 Z. The disagreement never exceeds 2 $`\sigma `$ but it is worrying that the observed slope is considerably more steep than the theoretical one.

Several effects might be at the origin of that discrepancy. A shallower stellar IMF than the one of Salpeter could certainly produce a larger NIb,c/NII ratio; but to match the NIb,c/NII ratio at Z, it would take a slope X$`<`$1 for the stellar IMF (see Fig. 3), i.e. much smaller than all current empirical evidence. On the other hand, an IMF with a metallicity dependent slope could obviously fit to the data. However, such a variability of the IMF is not supported by observations (see e.g. Kroupa 2002 and references therein).

As pointed out in Sect. 2, the theoretical predictions are based on the single star models of Heger et al. (2003a), calculated with mass loss. Other ingredients, not taken into account in those models, may affect the NIb,c/NII ratio. For instance, mass loss through Roche lobe overflow in a binary system is a well known channel to form WR stars (e.g. Chiosi and Maeder 1986), which are progenitors of SNIb,c. This channel was thought to be predominant at low metallicities: Maeder & Meynet (1994) found that all WR stars in the Small Magellanic Cloud are in binaries, but argued that the WR/O ratio in different metallicity environments suggests that the fraction of the WR stars due to that channel is small at high metallicities. Recently, Foellmi et al. (2003a) performed a thorough study of the WR population of the SMC, and found that, contrary to theoretical expectations, the binary fraction of WR in that galaxy is $``$40$`\%`$ rather than $``$100$`\%`$ theoretically expected. Also Foellmi et al. (2003b) found that the corresponding value for the LMC (which is $``$3 times more metal-rich than the SMC) is $``$30$`\%`$. These findings weaken even more the role of binarity in WR formation at high metallicities. In fact, assuming that the binary channel contributes by 20$`\%`$ to the Milky Way population of WR (Foellmi et al. 2003a,b), that binarity reduces the average WR mass to 19 M (for solar metallicity stars) and taking into account the results displayed in Fig. 3, one sees that the expected NIb,c/NII value at solar metallicity is 0.23, i.e. still lower than observed in solar metallicity environments. Moreover, the contribution of binaries to the WR fraction and to the NIb,c/NII ratio should be independent of metallicity and cannot therefore account for the mismatch obtained in Fig. 2 (bottom) between observed and theoretical slopes.

Models of rotating massive stars (combined with mass loss) have been recently developed by a few groups (Maeder & Meynet 2002, Heger et al. 2003b). Meynet & Maeder (2003) find that rotation produces stars with more massive convective cores, higher effective temperatures and smaller mass left at the end of the evolution, compared to non-rotating ones. In particular, they find that the minimal mass for a star to become WR is 37 M for non-rotating stars (i.e. close to the 34 M suggested by Heger et al. 2003a), but only 22 M for stars rotating initially at 300 km/s. In fact, such stars become only WNL stars, i.e. they keep a trace amount of their original hydrogen envelope. An inspection of Figs. 9 and 10 of Meynet & Maeder (2003) shows that stars rotating at 300 km/s with slightly higher mass, around 25 M, become WNE stars losing their entire H-envelope.

This reduction in MIb,c obviously affects the NIb,c/NII ratio. In Fig. 2 (bottom) we also plot the NIb,c/NII ratio assuming MIb,c=22 and 25 M, respectively, for the Salpeter and for the Kroupa et al. (1993) IMFs and for solar metallicity only (open and filled squares, respectively). It can be seen that there is a much better agreement with the observations. As argued by Meynet & Maeder (2003), rotation provides a better fit to various observations concerning WR stars, like statistics of WR subtypes and surface abundances. According to the results displayed in Fig. 2, it also provides a nice explanation to the observed NIb,c/NII ratio, at least for solar metallicity stars. Indeed, if it is assumed that MIb,c is as high as 34 M, it requires unnaturaly low values of the slope X of a power-law IMF in order to obtain the observed NIb,c/NII ratio, as shown in Fig. 3. Only for values of MIb,c at least as low as 22-25 M, the observed NIb,c/NII ratio can be reproduced with “reasonable” values of X. We suggest that this is a further argument in favour of the rotating star models of Meynet & Maeder (2003).

Assuming that the fit of Eq. (3) represents well the metallicity dependence of the NIb,c/NII ratio, and that the only factor affecting that dependence is the variation of MIb,c with metallicity, one may use Eq. (1) to determine that variation. This is done in Fig. 4, for two IMFs. The thick portions of the curves indicate the region of galaxian metallicity (or blue magnitude) where NIb,c/NII has been observationally determined (i.e. between MB values of -23 and -14 in Fig. 2c). It can be seen that MIb,c varies from 40-50 M at Z$``$0.3 Z to slightly below 20 M at Z$``$2 Z. This observationally determined variation of MIb,c may offer potentially important constraints to models of massive star evolution.

We note that Bressan et al. (2002) find that the ratio NIb,c/NII in Seyfert galaxies is $``$1, i.e. about three to four times larger than in normal galaxies. They explore several possible reasons for such a high ratio, including metallicity effects on MIb,c, and they conclude that the most probable one is a recent starburst, of age smaller than the lifetime of the smallest core collapse supernova (in which case Eq. (1) is invalid). They base their quantitative estimates on non-rotating models with mass loss. We think their explanation certainly holds, but rotating models may alleviate the constraints on the age of the starburst they obtain in their Fig. 3.

## 4 Relative supernova frequencies in galaxies

Supernova frequencies are usually given for galaxies of different morphological types, in units of SNu (i.e. in number of SN per 100 yr and per 1010 LB,⊙). The rational behind this definition of frequency units is that the galaxy’s blue luminosity measures its total mass, as also suggested by the Tully-Fisher relation in the blue. Note that in starburst galaxies the blue luminosity traces rather the young stellar population, while in normal galaxies there is a significant contribution from the old, metal poor, population.

As expected then, the SN frequencies expressed in SNu depend little (if any at all) on the galaxy’s morphological type (e.g. Cappellaro et al. 1999, Turatto 2000), as shown on Fig. 4a. For the same reason, the ratio NIb,c/NII is also insensitive to the galaxy’s morphological type (Fig. 4b) and it is of little utility in galaxy studies. We only plot those ratios in Fig. 4b as a consistency check of our results with those of previous surveys. On the other hand, Cappellaro et al. (1999) have shown that the frequencies of core collapse supernovae (II+Ib,c) depend significantly on the galaxy’s U-V colour, in the sense that higher frequencies correspond to smaller U-V values, i.e. to younger galaxies; they suggested that this dependence reflects directly the dependence of the core collapse SN rate on the current star formation rate. If confirmed, this finding would allow then to probe the star formation rates of late type galaxies, i.e. their current status (provided that the relevant uncertainties, that is the effect of extinction on U-V indices, are adequatly taken into account).

We suggest here that instead of SN frequencies, ratios between subtypes of core collapse supernovae coud be used as a probe of the metallicity of the host galaxies, which is a measure of the integrated star formation activity of the galaxy, i.e. of its past history. Measurements of numbers of SN subtypes at higher redshifts would then allow to probe in much greater detail the history of metallicity evolution in the Universe. For that purpose, it would not be sufficient to simply plot the NIb,c/NII ratio as a function of redshift (as in e.g. Heger et al. 2003a). A detailed survey, equivalent to the one presented in Fig. 2, would be necessary: plotting the NIb,c/NII ratio as a function of the (rest frame) blue magnitude of galaxies at different redshifts, would enable us to check whether the variation of that ratio with blue magnitude was the same in the past; since this relationship is an indirect one (i.e via the metallicity of the host galaxy), such a survey would offer also clues about the metallicity-luminosity relation in the past.

Such a detailed survey would be possible with future instruments, such as the SNAP satellite (see web page at: http://snap.lbl.gov/). SNAP is expected to detect $``$2000 SNIa up to redshifts of 1.7. Taking into account that core collapse supernovae are fainter (by $``$ two magnitudes) than SNIa and 2-3 times as frequent, one sees that SNAP will be able to construct large enough samples for the purposes of the proposed study at several redshift ranges up to z$``$1.

## 5 Summary

This work explores (some of) the consequences of metallicity dependent mass loss of massive stars on the statistics of core collapse supernovae of various subtypes, and in particular on the ratio of type Ib,c vs. type II supernovae. According to current models of massive stars, this ratio should increase with increasing metallicity (Heger et al. 2003a). The exact form of that relationship depends (a little) on the adopted stellar IMF and mostly on the value of MIb,c, the minimal value of initial mass for a star to explode as SNIb,c.

We show first that the observed NIb,c/NII does indeed display a metallicity dependence. For that purpose we construct a statistically significant sample of Ib,c and II supernovae as a function of their host galaxy’s blue magnitude and we use the well-known metallicity-luminosity relation for late type galaxies (Garnett 2002).

However, the observed NIb,c/NII vs Z relationship is steeper than theoretically expected, if one adopts non-rotating star models (from e. g. Heger et al. 2003a or Meynet & Maeder 203). We argue that neither the assumptions of a metallicity dependent IMF or of an enhanced fraction of Ib,c supernovae originating in binary systems offer viable solutions to the problem. We suggest that rotating stellar models with mass loss, such as those recently calculated by Meynet & Maeder (2003) offer a much better quantitative explanation. We base our argument on models of solar initial metallicity and we urge calculations at lower and higher than solar metallicities in order to substantiate our conclusions.

Furthermore, we suggest that surveys of core collapse supernovae with the SNAP satellite will allow to probe the effect of metallicity on supernovae properties during a large fraction of the past history of the Universe, at least up to redshift $`z`$1.

###### Acknowledgements.

We are grateful to the referee, G. Meynet, for his critical comments and suggestions, which improved the paper substantially.


format: text

text:
# Multiple Folding Pathways of the SH3 domain

## 1 INTRODUCTION

Recent experimental studies indicate that several proteins exhibit simultaneously a variety of intermediates and folding pathways. Kiefhaber identified at low denaturant concentration a fast pathway (50 ms) in the folding of lyzosyme with no intermediates and a slow phase (420 ms) with well-populated intermediates. Choe et al. observed the formation of a kinetic intermediate in the folding of villin 14T upon decreasing the temperature, and Silverman et al. observed the extinction of a slow phase in the folding of the P4-P6 domain upon changes in ion concentration. Kitahara et al. studied a pressure-stabilized intermediate of ubiquitin, identified as an off-pathway intermediate in previous kinetics experiments at basic conditions. All these studies suggest that environmental conditions favor some folding pathways over others.

Major theoretical efforts in the study of protein folding have focused on small, single domain proteins . It is found in experiments that these proteins undergo folding transition with no accumulation of kinetic intermediates in the accessible range of experimental conditions. However, other kinetics studies of two-state proteins suggest the presence of short-lived intermediates that cannot be directly detected experimentally. Recently, Sánchez et al. explained the curved Chevron plots — the nonlinear dependence of folding and unfolding rates on denaturant concentration — of 17 selected proteins by assuming the presence of an intermediate state. Led by these studies, we hypothesize that single domain proteins may exhibit intermediates in the folding transition under suitable environmental conditions.

To test our hypothesis, we perform a molecular dynamics study of the folding pathways of the c-Crk SH3 domain (PDB access code 1cka). The SH3 domain is a family of small globular proteins which has been extensively studied in kinetics and thermodynamics experiments. We select the c-Crk SH3 domain ($`57`$ residues) as the SH3 domain representative and perform molecular dynamics simulations over a broad range of temperatures. We determine the kinetic partition temperature $`T_{KP}`$ below which the model protein exhibits slow folding pathways and above which the protein undergoes a cooperative folding transition with no accumulation of intermediates. Below $`T_{KP}`$, we study the presence of intermediates in the slow folding pathways and resolve their structure. We find that one of the intermediates populates the folding transition for temperatures as high as $`T_{KP}`$. We discuss the relevance of our results in light of recent experimental evidence.

## 2 RESULTS

The SH3 domain is a $`\beta `$-sheet protein (upper triangle of Fig. 1a). Our previous thermodynamic studies of the c–Crk SH3 domain revealed only two stable states at equilibrium conditions: folded and unfolded. Both states coexist with equal probability at the folding transition temperature, $`T_F=0.626`$, at which the temperature dependence of the potential energy has a sharp change, and the specific heat has a maximum (experimentally, this temperature corresponds to $`67^{}`$C). Thus, our model reproduces the experimentally-determined thermodynamics of the SH3 domain.

### 2.1 Initial Unfolded Ensemble

Our initially unfolded ensemble consists of 1100 protein conformations that we sample from a long equilibrium simulation at a high temperature $`T_0=1.0`$ at equal time intervals of $`10^4`$ time units (t.u). This time separation is long enough to ensure that the sampled conformations have low structural similarity among themselves. We calculate the frequency map — the plot of the probability of any two amino acids forming a contact — of this unfolded ensemble (lower triangle of Fig. 1a). At $`T=1.0`$, only nearest and next nearest contacts have high frequency, and the frequency decreases rapidly with the sequence separation between the amino acids.

When we quench the system from $`T=1.0`$ to a target temperature, $`T_{\mathrm{target}}`$ (see Methods), the system relaxes in approximately $`1500`$ t.u. Due to the finite size of our heat bath, the heat released by the protein upon folding increases the final temperature of the system by $`0.03`$ units above $`T_{\mathrm{target}}`$. After relaxation, the protein stays for a certain time in the unfolded state, then undergoes a folding transition. During this time interval, the protein explores unfolded conformations, and we calculate the frequency map of the unfolded state for different target temperatures.

At $`T_{\mathrm{target}}=T_F`$, the secondary structure is unstable (Fig. 1b), with average frequency $`\overline{f}=0.24`$ (see Methods). Successful folding requires the cooperative formation of contacts throughout the protein in a nucleation process. At $`T_{\mathrm{target}}=0.54`$, the secondary structure is more stable (Fig. 1c, $`\overline{f}=0.50`$). Thus, the conformational search for the native state (NS) is optimized by limiting the search to the formation of a sufficient number of long range contacts. At $`T_{\mathrm{target}}=0.33`$, the lowest temperature studied, secondary structure elements form during the rapid collapse of the model protein in the first $`1500`$ t.u. (Fig. 1d, $`\overline{f}=0.73`$). During collapse, some tertiary contacts — contacts between secondary elements — may also form. The formation of these contacts prior to the proper arrangement of secondary structure elements may lead the protein model to a kinetic trap. Finally, folding proceeds at this temperature through a thermally activated search for the NS.

### 2.2 Kinetic Partition Temperature

In order to determine the temperature below which we can distinguish fast and slow folding pathways, we compute the distribution of folding times $`p(t_F,T)`$ (Fig. 2a–e), as well as the average $`t_F`$ (Fig. 2f) and standard deviation $`\sigma _F`$. The ratio $`r(T)t_F/\sigma _F`$ measures the average folding time in units of the standard deviation $`\sigma _F`$. This quantity is particularly useful when the value of the standard deviation correlates with the value of the average as we change $`T_{\mathrm{target}}`$. For instance, single-exponential distributions $`e^{t_F/t_F}/t_F`$ have $`r1`$.

We expect $`r1`$ for $`T_{\mathrm{target}}>T_F`$, because at these high temperatures the folding transitions become rare events and are single-exponential distributed. As we decrease $`T_{\mathrm{target}}`$, we expect $`r>1`$ just below $`T_F`$, because the folded state becomes more stable than the unfolded state, and the folding transitions are favored. Distributions with $`r>1`$ indicate a narrow distribution centered in $`t_F`$, so that most of the simulations undergo a folding transition for times of the order of the average folding time. However, if we continue decreasing $`T_{\mathrm{target}}`$, we expect some folding transitions to be kinetically trapped, and the folding time distribution will spread over several orders of magnitudes. Such distributions have $`r<1`$. Thus, there is a temperature below $`T_F`$ where the maximum of $`r(T)`$ occurs, and which signals the onset of slow folding pathways. We use the maximum of $`r(T)`$ to calculate $`T_{KP}`$.

Fig. 2g suggests that $`T_{KP}=0.54`$, which corresponds to a maximally compact distribution of folding times***Assuming a linear relation between experimental and simulated temperatures and taking into account that $`T_F=67^{}`$C, we estimate $`T_{KP}20^{}`$C (Fig. 2d). We find that the ratio approaches one as we increase the temperature above $`T_{KP}`$, and the distribution of folding times approximates a single-exponential distribution. In particular, the distribution of folding times fits the single-exponential distribution $`e^{t_F/t_F}/t_F`$ for $`T=0.64`$, the closest temperature to $`T_F`$ that we study. The ratio $`r(T)`$ decreases monotonically below $`T_{KP}`$, indicating that the distribution of folding times spreads over several orders of magnitude. This is the consequence of an increasing fraction of folding simulations kinetically trapped (Fig. 2a-b). The average folding time $`t_F`$ is minimal not at $`T_{KP}`$, but at a lower temperature $`T_{t_F}=0.49`$ (Fig. 2f). At this temperature, we find that the protein becomes temporarily trapped in approximately $`7\%`$ of the folding transitions. On the other hand, the remaining simulations undergo a folding transition much faster, thus minimizing $`t_F`$. Interestingly, $`r(T_{t_F})1.0`$, even though the distribution of folding times at this temperature is non-exponential.

### 2.3 Folding Pathways

Below $`T_{KP}`$, an increasing fraction of the simulations undergo folding transitions that take a time up to three orders of magnitude above the minimal $`t_F`$. In addition, $`t_F`$ increases dramatically (Fig. 2f). At the lowest temperatures studied, we distinguish between the majority of simulations that undergo a fast folding transition (the fast pathway) and the rest of the simulations that undergo folding transitions with folding times spanning three orders of magnitude (the slow pathways). At the low temperature $`T=0.33`$, the potential energy of the fast pathway has on average the same time evolution of all the simulations at $`T_{KP}=0.54`$, indicating that there are no kinetic traps in the fast pathway.

For each folding simulation that belongs to the slow pathways, we sample the potential energy at equal time intervals of $`100`$ t.u. until folding is finished (see Methods). Then, we collect all potential energy values and construct a distribution of potential energies. We find that below $`T=0.43`$, the distribution is markedly bimodal (Fig. 3a). The positions of the two peaks along the energy coordinate do not correspond to the equilibrium potential energy value of the folded state (Fig. 3b). Therefore we hypothesize the existence of two intermediates in the slow pathways. We denote the two putative intermediates as $`I_1`$ and $`I_2`$ for the high energy and low energy peaks, respectively. As temperature decreases, the peaks shift to lower energies, but the energy difference between the two peaks, approximately six energy units, remains constant (Fig. 3b). A constant energy difference implies that the two putative intermediates differ by a specific set of native contacts. As temperature decreases, other contacts not belonging to this set become more stable and are responsible for the overall energy decrease. At $`T=0.33`$, we record the distribution of survival times for both intermediates and find that they fit a single-exponential distribution, supporting the hypothesis that each intermediate is a local free energy minima and has a major free energy barrier (Fig. 3c).

To further test the single free energy barrier hypothesis, we select a typical conformation representing intermediate $`I_2`$ and perform 200 folding simulations, each with a different set of initial velocities for a set of temperatures in the range $`0.33T0.52`$. For each simulation, we record the time that the protein stays in the intermediate and find that the average survival time fits the Arrhenius law for temperatures below $`T=0.44`$ (Fig. 3d). This upper bound temperature roughly coincides with the temperature $`T=0.43`$ below which $`I_2`$ becomes noticeable in the histogram of potential energies (Fig. 3a). This result indicates that the free energy barrier to overcome intermediate $`I_2`$ becomes independent of temperature for low temperatures, or analogously, that the same set of native contacts must form (or break) to overcome the intermediate.

Next, we determine the structure of the two intermediates. For each intermediate, we randomly select three conformations and find that they are structurally similar. Conformations belonging to intermediate $`I_1`$ have a set of long-range contacts (C1) with a high occupancy and a set of long-range contacts (C2) with no occupancy at all (Fig. 3e). Contacts in C1 represent a $`\beta `$-sheet made by three strands: the two termini and the strand following the RT-loop, which we name strand “A” (see $`I_1`$ in Fig. 4). Contacts in C2 represent the base of the n-Src loop and the contacts between the RT-loop and the distal hairpin (see $`I_2`$ in Fig. 4). In addition, $`I_1`$ has a set of medium-range contacts (C3) with high occupancy (Fig. 3e) representing the distal hairpin and a part of the n-Src loop. For a slow folding transition, the $`\beta `$-sheet (C1 contacts) forms in the early events and strand “A” can no longer move freely. This constrained motion prevents strand “A” from forming contacts with one of the strands of the distal hairpin, which we name as strand “B” (see $`I_2`$ in Fig. 4). Similarly, strand “B” cannot move freely because it is a part of C3. The missing contacts between strand “A” and strand “B” are the contacts that form the base of the n-Src loop.

Conformational changes leading the protein away from intermediate $`I_1`$ involve either dissociation of the $`\beta `$-sheet, thus breaking some contacts of C1, or dissociation of the distal hairpin, thus breaking some contacts of C3. We find that the latter dissociation may lead the protein conformation to intermediate $`I_2`$. Intermediate $`I_2`$ has contacts of C1, but lacks the set of contacts (C4) that form the base of the the distal hairpin (see NS in Fig. 4).

Once we identify the structure of the intermediates, we investigate whether intermediate $`I_1`$ is present at larger temperatures when no distinction can be made concerning fast and slow folding pathways. To test this hypothesis, we sample the protein conformation during the folding transition at equal time intervals of $`60`$ t.u. for each of the $`1100`$ simulations, and compare these conformations to intermediate $`I_1`$ with a similarity score function (see Methods). For each folding transition, we record only the highest value of the similarity score, thus obtaining $`1100`$ highest score values. At $`T_{KP}`$, the histogram of the highest scores is bimodal, with $`25\%`$ of the folding simulations passing through intermediate $`I_1`$ (Fig. 3f). We find that at $`T_{KP}`$, simulations that undergo the folding transition through $`I_1`$ show kinetics of folding no different than those of the rest of simulations.

## 3 DISCUSSION

It was shown that the simplified protein model and interaction potentials that we use here reproduced in a certain range of temperatures the
experimentally-determined two-state thermodynamics of the SH3 domain. The qualitative predictive power of the model encouraged us to study the folding kinetics in a wider range of temperatures. From our relaxation studies of the initial unfolded ensemble, we observe that the structure of the unfolded state is highly sensitive to $`T_{\mathrm{target}}`$. The role of the unfolded state in determining the folding kinetics has already been pointed out in recent experimental and theoretical studies. We observe nucleation, folding with minimal kinetic barriers, and thermally activated mechanisms for the different observed unfolded states.

In previous studies, various methods have been developed to determine the temperature that signals the onset of slow folding pathways. Socci et al. determined a glass transition temperature, $`T_g`$, at which the average folding time is half way between $`t_{\mathrm{min}}`$ and $`t_{\mathrm{max}}`$, where $`t_{\mathrm{min}}`$ is the minimun average folding time and $`t_{\mathrm{max}}`$ is the total simulation time. This method is sensitive to the a priori selected $`t_{\mathrm{max}}`$. The authors varied $`t_{\mathrm{max}}`$ in the range $`0.27\times 10^9<t_{\mathrm{max}}<0.960\times 10^9`$, and they found a $`10\%`$ error in the calculation of $`T_g`$. Also, Gutin et al. estimated a critical temperature, $`T_c`$, at which the temperature dependence of the equilibrium potential energy leveled off. From their results, one can evaluate a $`20\%`$ error in their calculation of $`T_c`$. Both $`T_g`$ and $`T_c`$ are temperatures that authors use to characterize the onset of multiple folding pathways. In our study we use $`T_{KP}`$, which signals the breaking of time translational invariance of equilibrium measurements for temperatures below this value. We estimate a $`2\%`$ error in our calculation of $`T_{KP}`$ from uncertainties in the location of $`T_{KP}`$ in Fig. 2g.

At $`T_{KP}`$, secondary structure elements are partially stable, and the search for the NS reduces to the formation of tertiary contacts. Furthermore, $`T_{KP}`$ is a relatively high temperature that prevents the stabilization of improper arrangements of the protein conformation, thus minimizing the occurrence of kinetic traps. Below $`T_{KP}`$, the model protein exhibits two intermediates with well-defined structural characteristics. This modest number of intermediates is a direct consequence of the prevention of non-native contacts. This prevention reduces dramatically the number of protein conformations. Furthermore, since a low energy value implies that most of the native interactions have formed, there are few conformations having both low energy and structural differences with the NS.

It is found experimentally  that proteins exhibit only a discrete set of intermediates. Even though in real proteins amino acids that do not form a native contact may still attract each other, experimental and theoretical studies confirm that native contacts have a leading role in the folding transition. Protein engineering experiments show that transition states in two-state globular proteins are mostly stabilized by native interactions. To quantitatively determine the importance of native interactions in the folding transition, Paci et al. studied the transition states of three two-state proteins with a full-atom model. They found that on average, native interactions accounted for approximately $`83\%`$ of the total energy of the transition states. Of relevance to our studies of the SH3 domain are the full-atom study and the protein engineering experiments showing that the transition state of the src-SH3 domain protein is determined by the NS. On the other hand, evidence exists that in some proteins, non-native contacts are responsible for the presence of intermediates. In their study of the homologous Im7 and Im9 proteins, Capaldi et al. identified a set of non-native interactions responsible for a intermediate state in the folding transition of Im7 protein. Mirny et al. performed Monte Carlo simulations of two different sequences with the same NS in the $`3\times 3`$ lattice. One sequence presented a series of pathways with misfolded states due to non-native interactions.

We investigate the kinetics of formation of the two intermediates in a wide range of temperatures. At low temperatures, simulations that undergo folding through intermediate $`I_1`$ reveal that contacts between the two termini form earlier than the contacts belonging to the folding nucleus. This result coincides with an off-lattice study of a 36-monomer protein by Abkevich et al. In this study, the authors found an intermediate in the folding transition of their model protein. Inspection of the intermediate revealed no nucleus contacts, but a different set of long-range contacts already formed. In addition, we learned of the work by Viguera et al. after completion of our study. They reported that a mutant of the $`\alpha `$-spectrin SH3 domain undergoes a folding transition through one intermediate. The authors observed that the newly-introduced long-range contacts had already been formed in the denaturated state, preceding the formation of the transition state of the protein. Thus, environmental conditions that favor stabilization of long-range contacts other than the nucleus contacts may induce intermediates in the folding transition.

Alternatively, short-range contacts in key positions of the protein structure may also be responsible for slow folding pathways. After completion of our study, Karanicolas et al. reported their studies on the Gō model of the forming binding protein WW domain. The authors found a slow folding pathway in the model protein, and a cluster of four short-range native contacts that are responsible for this pathway. However, the authors observed that it was the absence, not the presence, of these native contacts in the unfolded state that generated bi-phasic folding kinetics. Thus, environmental conditions that favor destabilization of short range contacts may promote the formation of intermediate states in the folding transition.

We also investigate the survival time of intermediate $`I_2`$, and find that the free energy barrier separating $`I_2`$ from the NS is independent of temperature. Thus, the average survival time follows Arrhenius kinetics. The value of the free energy barrier is approximately $`5.85`$ energy units, indicating that about six native contacts break when the protein conformation reaches the transition state that separates $`I_2`$ from the NS. At the low temperatures where intermediates $`I_1`$ and $`I_2`$ are noticeable, thermal fluctuations are still large enough so that the observed survival times of $`I_2`$ should be much smaller, if only any six native contacts were to break. Thus we hypothesize that it is allways the same set of native contact that must break in the transition $`I_2NS`$. Our observations of the transition $`I_1I_2`$ support this hypothesis. In this transition, we find that the set of contacts C4 allways breaks.

At $`T_{KP}`$, we do not detect the intermediates from kinetics measurements of the average folding time, or analogously, from the folding rate. Thus we analyze the folding transition with the similarity score function that tests the presence of intermediate $`I_1`$. Then we find that this intermediate is populated in $`25\%`$ of the folding transitions. In a recent study, Gorski et al. reported the existence of an intermediate in the folding transition of protein Im9 under acidic conditions ($`pH=5.5`$). This finding led authors to formulate the hypothesis that Im9 has an intermediate at normal conditions ($`pH=7.0`$), but it is too unstable to be detected with current kinetic experimental techniques. Interestingly, the homologous protein Im7 ($`60\%`$ sequence identity) undergoes folding transition through an intermediate in all tested experimental conditions, supporting the authors’ hypothesis. Thus, changes in both the environmental conditions and the amino acid sequence may uncover hidden intermediates in the folding transition of a two-state protein. In addition, a detailed study at $`T_{KP}`$ may reveal the intermediates. This is particularly useful for computer simulations, because simulations at low temperatures when intermediates are easily identifiable may require several orders of magnitude longer than simulations at $`T_{KP}`$.

### 3.1 Conclusion

We perform molecular dynamics analysis of the folding transition of the Gō model of the c-Crk SH3 domain in a broad range of temperatures. At the folding transition temperature, we observe that only the folded and unfolded states are populated. However, as we decrease the temperature, parameters monitoring the folding process such as potential energy and root mean square distance with respect to the native state, $`rmsd`$, suggest the presence of intermediates. We determine the kinetic partition temperature $`T_{KP}`$ below which we observe two folding intermediates, $`I_1`$ and $`I_2`$, and above which we do not observe accumulation of intermediates. Below $`T_{KP}`$, intermediate $`I_1`$ forms when the two termini and the strand following the RT-loop form a $`\beta `$-sheet, prior to the formation of the folding nucleus. This intermediate effectively splits the folding transition into fast and slow folding pathways. Dissociation of part of the $`\beta `$-sheet leads the protein to the native state. We also find that stabilization of this $`\beta `$-sheet and subsequent dissociation of the distal hairpin may lead to intermediate $`I_2`$.

The key structural characteristics of intermediate $`I_1`$ allow us to define a similarity score function that probes the presence of the intermediate in a folding transition. We find that $`I_1`$ is populated even at $`T_{KP}`$. This result suggests that one can obtain information regarding the existence of putative intermediates by studying the folding trajectories at $`T_{KP}`$. However, at this temperature, no intermediates are noticeable if one limits the analysis only to the distribution of folding times.

We observe that the folding pathways of the model SH3 domain are highly sensitive to temperature, suggesting the important role of the environmental conditions in determining the folding mechanism. Our findings suggest that the SH3 domain, a two-state folder, may exhibit stable intermediates under extreme experimental conditions, such as very low temperatures.

## 4 MATERIALS AND METHODS

### 4.1 Model Protein and Interactions

We adopt a coarse-grained description of the protein by which each amino acid is reduced to its Cβ atom (Cα in case of Gly). Details of the model, the surrounding heat bath, and the selection of structural parameters are discussed in detail in a previous study. The selection of the set of interaction parameters among amino acids is of crucial importance for the resulting folding kinetics of the model protein . Experimental and theoretical studies of globular proteins  suggest that native topology is the principal determinant of the folding mechanism. Thus, we employ a variant of the Gō model of interactions — a model based solely on the native topology — in which we prevent formation of non-native interactions, since we are solely interested in the role that native topology and native interactions may have in the formation of intermediates. We perform simulations and monitor the time evolution of the protein and the heat bath with the discrete molecular dynamics algorithm . The higher performance of this algorithm over conventional molecular dynamics allows one to increase the computational speed up to three orders of magnitude.

### 4.2 Frequencies and Folding Simulations

To calculate the frequency map at $`T=1.0`$, we probe the presence of the native contacts in each of the $`1100`$ initially unfolded conformations. Then, we compute the probability of each native contact to be present. To calculate the frequency map at $`T_{\mathrm{target}}`$, we select one particular folding transition and we probe the presence of the native contacts during the time interval that spans after the initial relaxation and before the simulation reaches the folding time $`t_F`$. To compute $`t_F`$, we stop the folding simulation when $`90\%`$ of the native contacts form. Then, we trace back the folding trajectory and record $`t_F`$ when the root mean square distance with respect to the native state, $`rmsd`$, becomes smaller than $`3`$Å. We consider all protein conformations occurring for $`t>t_F`$ as belonging to the folded state and of no relevance to the folding transition.

### 4.3 Similarity Score Function

We introduce the similarity score function, $`S=(a/23)(15b)/15`$, where $`a`$ is the number of native contacts belonging to the set of contacts C1, and $`b`$ is the number of native contacts belonging to set C2 (Fig. 3e). C1 has $`23`$ contacts and C2 has $`15`$ contacts. If the protein is unfolded, then $`ab0`$, thus $`S0`$. Similarly, if the protein is folded, then $`a23`$ and $`b15`$, thus $`S0`$ again. Finally, if the protein adopts the intermediate $`I_1`$ structure, then $`a23`$ and $`b0`$, thus $`S1`$.

We acknowledge E. I. Shakhnovich for insightful discussions, R. Badzey for careful reading of the manuscript, and the Petroleum Research Fund for support.

Figure 1: (a) Upper triangle: c–Crk SH3 domain contact map with $`160`$ native contacts. The secondary structure elements are the clusters of contacts that are organized perpendicularly to the map diagonal. Long-range contacts between the two termini are enclosed in the circle, and long-range contacts between the RT-loop and the distal hairpin are enclosed in the square. Lower triangle: the frequency map for the initial set of 1100 unfolded conformations at $`T=1.0`$. (b) Frequency map of the unfolded state at $`T=T_F=0.626`$. We compute the frequencies for a particular folding transition, whose potential energy trajectory we show in the inset (see Methods). Same for (c) $`T=0.54`$ and (d) $`T=0.33`$, the lowest temperature studied.

Figure 2: (a–e) Histograms of folding times for selected temperatures. At $`T=0.33`$ and $`T=0.36`$, the two lowest temperatures studied, histograms have a maximum for long folding times ($`^{\mathrm{}}`$), which suggests the existence of putative intermediates. At $`T=0.33`$, a maximum in the histogram ($`^+`$), not present at $`T=0.36`$, corresponds to short lived kinetic traps. The distributions of folding times are unimodal at higher temperatures. At $`T=0.54`$, the histogram is compact, and has no tail of long folding times. At $`T=0.64`$, the histogram fits a single-exponential distribution $`e^{t_F/t_F}/t_F`$ for times larger than the relaxation time of $`1500`$ t.u. (dashed line). We estimate the errors of the histogram bars as the square root of each bar. (f) Average folding time versus temperature. Each dot represents the folding time for a particular folding transition. (g) Ratio $`r`$ of the average and the standard deviation, $`r=t_F/\sigma _F`$, for the distribution of folding times. The ratio approaches one above $`T_{KP}`$ and zero below $`T_{KP}`$. The ratio is maximal at $`T_{KP}`$, indicating a compact distribution of folding times at this temperature.

Figure 3: (a) Distributions of the potential energies of the slow folding pathways for temperatures below $`T=0.43`$. The distributions are bimodal, suggesting two putative intermediates $`I_1`$ and $`I_2`$. (b) The potential energy of the distribution peaks (✳ and $`\mathrm{}`$) increases with temperature, but the energy difference between peaks remains constant. The energy of the peaks is significantly larger than the equilibrium energy of the folded state ($``$). (c) Distributions of survival times at $`T=0.33`$ for the high energy intermediate $`I_1`$ (✳), $`t_F=1.81\times 10^6`$ and $`\sigma =1.85\times 10^6`$, and the low energy intermediate $`I_2`$ ($`\mathrm{}`$), $`t_F=2.47\times 10^6`$ and $`\sigma =2.43\times 10^6`$, fit to single-exponential distributions. (d) Arrhenius-fit of the average survival time of intermediate $`I_2`$ below $`T=0.44`$. This upper bound temperature coincides with the temperature below which the distribution of the potential energies (Fig. 3a) of the slow folding pathways becomes bimodal. (e) Upper triangle: Absent contacts (filled squares) and present contacts (crosses, “C4”) in intermediate $`I_1`$. Upon the transition $`I_1I_2`$, these contacts reverse their presence (so that the filled squares are the present contacts and the crosses are the absent contacts). There are five more squares than crosses, which roughly accounts for the difference of six energy units between the two intermediates. Lower triangle: long-range contacts “C1” are present in intermediate $`I_1`$, and long-range contacts “C2” are absent. There are $`23`$ contacts in “C1” and $`15`$ contacts in “C2”. (f) Probability that a folding transition at $`T=T_{KP}=0.54`$ contains a protein conformation with similarity $`S`$ to intermediate $`I_1`$ (see Methods).

Figure 4: Schematic diagram of fast and slow folding pathways. At $`T=0.33`$, approximately in $`15\%`$ of the simulations, the model protein undergoes a folding transition through the slow folding pathways. We show the protein structure in $`I_1`$ and $`I_2`$ using the secondary structural elements of the native state, although some of these elements are not formed. In intermediate $`I_1`$, both termini and the strand A form a $`\beta `$-sheet (in the ellipse). The corresponding set of native contacts is C1. Dissociation of the $`\beta `$-sheet leads to rearrangements of the protein conformation and successful folding to the native state (NS). However, dissociation of the distal hairpin (in red) leads to more localized rearrangements that may lead the protein to intermediate $`I_2`$. Upon $`I_1I_2`$ transition, contacts of C2 (the two ellipses in $`I_2`$) form, but contacts of C4 (the ellipse in NS) break.


text:
# (Vanishing) Twist in the Saddle-Centre and Period-Doubling Bifurcation

## 1 Introduction

The dynamics near a periodic orbit of a Hamiltonian system can be studied in terms of a local Poincaré section transversal to the orbit. In two degrees of freedom the first return map restricted to the surface of constant energy is an area preserving map with a fixed (or periodic) point. If the multipliers $`\mu _i`$, $`i=1,2`$ of the fixed point have modulus 1 (but are not equal to $`\pm 1`$) the fixed point is called elliptic. Then $`\mu _1=\overline{\mu _2}`$ and we can write $`\mu =\mathrm{exp}(2\pi i\omega )`$ with the rotation number $`0<\omega <1/2`$ of the periodic orbit. If the periodic orbit is elliptic and the rotation number is irrational the map can (formally) be transformed to Birkhoff normal form which in action-angle variables $`(\phi ,I)`$ reads

$$(\phi ,I)(\phi +2\pi \mathrm{\Omega }(I),I).$$
(1)
The action $`I`$ is like the radial coordinate in polar coordinates, hence the map in normal form maps circles to circles by rotating them $`\mathrm{\Omega }(I)`$ times. The rotation number (or winding number) $`\mathrm{\Omega }(I)`$ near the periodic orbit can be expanded as

$$\mathrm{\Omega }(I)=\omega +\tau _0I+\frac{1}{2}\tau _1I^2+\mathrm{}$$

The twist (or torsion) $`\tau (I)`$ is the derivative of the rotation number with respect to the action,

$$\tau (I)=\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}I}(I)=\tau _0+\tau _1I+\mathrm{}$$

When the rotation number is a strictly monotone function of the action in some interval the map (1) restricted to the corresponding annulus is called a monotone twist map. Moser’s KAM theorem  states that the invariant torus $`I=I_0`$ of (1) persist under perturbation when its frequency $`\mathrm{\Omega }(I_0)`$ is diophantine and its twist $`\tau (I_0)`$ does not vanish. Arnold’s KAM theorem  is the same statement for flows where the nonvanishing of the twist corresponds to the isoeneregetic nondegeneracy conditions. A well known corollary of the KAM theorem is the stability of an elliptic fixed point in two degrees of freedom when $`\mathrm{\Omega }(0)=\omega 1/3,1/4`$ and the twist at the origin is non-vanishing, $`\tau (0)=\tau _00`$.

When the twist vanishes the perturbed dynamics can be more complicated. The stability of an elliptic point can be lost when its twist vanishes, see  for an example of an unstable elliptic point with $`\omega =1/5`$. The effects of vanishing twist away from the origin was first described by Howard , and the resulting effects have been observed in many examples . The probably most spectacular effect is the appearance of so-called meandering curves . The properties of non-twist maps also show interesting behaviour under renormalisation  and recently it has been shown  that an extension of the KAM theorem can also be proved in this context. In  it was finally shown that the vanishing of twist at the fixed point generically occurs in a one parameter family when the rotation number of a fixed point passes through the interval $`[1/4,1/3]`$. When the twist vanishes at the fixed point a twistless torus is created in a twistless bifurcation . After creation the twistless torus passes through resonances and in this way non-twist maps generically appear in one parameter families of area preserving maps. The truncated resonant Birkhoff normal form shows that this twistless torus eventually collides with a saddle-centre bifurcation that gives rise to the period 3 orbits that collide with the fixed point when $`\omega =1/3`$ . Such a connection between resonance and vanishing twist can also be found in 4 dimensional symplectic maps .

The techniques of  can also be applied to the higher order resonances, in particular for $`\omega =1/4`$. In this paper we study the two remaining generic bifurcations at even stronger resonance $`\omega =0,1/2`$. By definition the corresponding fixed point is not elliptic. The two cases will be denoted as the $`1:1`$ and $`1:2`$ resonance, or as the saddle-centre and period doubling bifurcation, respectively. The main result is that near the saddle-centre bifurcation the twist always vanishes, while it does not vanish near the period doubling bifurcation.

The method is based on the analysis of the resonant normal form, in which the Poincaré map is approximated by the time $`1`$ map of a one degree of freedom system, see e.g. . This normal form is an approximation, that is local near the bifurcation in parameter space and in phase space. At first we will completely ignore this, and just analyse the normal forms in the following two sections. In section 4 we will address the problem of non-locality in phase space, and also comment on the effect of higher order perturbations on the twistless tori. We will talk of invariant tori even though the invariant curves $`H(u,v)=h`$ of the normal form may not be compact. This will also be justified in Sec. 4. Finally we treat a saddle-centre bifurcation in the Hénon map as an example.

## 2 Saddle-Centre Bifurcation

The normal form of a Hamiltonian system with two degrees of freedom near the $`1:1`$ resonance has the form

$$H(u,v)=\frac{v^2}{2}+\frac{u^3}{3}+\epsilon u.$$
(2)
The coefficient of $`u^3`$ has been scaled so that it equals $`1/3`$. The variables $`u`$ and $`v`$ are canonically conjugate variables on a local transversal Poincaré section and $`\epsilon `$ is a parameter, typically the energy of the original system. The Poincaré map is given by the time $`1`$ map of the flow of $`H(u,v)`$. If the period $`T`$ is large, the time $`1`$ map advances little. The rotation number of the full system is the period $`1`$ divided by the period of the reduced one degree of freedom flow, $`\mathrm{\Omega }=1/T`$. The more familiar $`\mathrm{\Omega }=2\pi /T`$ is obtained when the time $`2\pi `$ map is taken instead of the time $`1`$ map, but the time $`1`$ map is more natural at least for the example of the Hénon map we are going to discuss. Since $`\mathrm{\Omega }`$ is determined by $`T`$, we now study in detail the period $`T`$ of the one degree of freedom system given by $`H`$.

The critical points and critical values of the energy map $`H:^2`$ and their dependence on $`\epsilon `$ give the main structure to the bifurcation. Instead of a one degree of freedom system $`H(u,v;\epsilon )`$ depending on the parameter $`\epsilon `$ one may consider $`H(u,v,\epsilon ,\theta )`$ as a Hamiltonian in $`^3\times S^1`$, with action $`\epsilon `$ and conjugate angle $`\theta `$. The set of critical values of the energy-momentum map $`(H,\epsilon ):^3\times S^1^2`$ is called the bifurcation diagram. A simple way to compute it is to find the critical values of the energy map of $`H(u,v)`$ and consider their parameter dependence on $`\epsilon `$. The Hamiltonian has critical points $`(u,v)=(\pm \sqrt{\epsilon },0)`$ and corresponding critical values $`h=2(\epsilon )^{3/2}/3`$. They exist when $`\epsilon <0`$ and the upper sign corresponds to a local minimum of $`H`$, while the lower sign gives a saddle. The corresponding phase portraits are shown in Fig. 1.

The dynamics is given by Hamiltons equation $`\dot{u}=v`$ and eliminating $`v`$ using the Hamiltonian gives a first order equation for $`u`$. After separation of variables the period of motion with energy $`h`$ is given by the elliptic integral

$$T(h,\epsilon )=\frac{\mathrm{d}u}{\sqrt{2h\frac{2}{3}u^32\epsilon u}},$$
(3)
where the integration path is encircling the interval on the real axis where the argument of the square root is positive. If there are two positive intervals either one can be taken, the result is the same. By scaling $`u=z(\sigma \epsilon )^{1/2}`$, where $`\sigma =\mathrm{sign}(\epsilon )`$ and introducing the one essential parameter

$$\gamma =\frac{3h}{2(\sigma \epsilon )^{3/2}},$$
(4)
the period $`T`$ is an elliptic integral on the curve

$$:w^2=P_3(z)=2\gamma z^33\sigma z.$$
(5)
The case $`\epsilon =0`$ has to be excluded in this scaling, but it is simple to treat it separately. We are mostly interested in the case where $`\sigma =1`$. The essential integral now reads

$$S(\gamma )=\frac{1}{w}dz$$
(6)
and it is related to the period by

$$T(h,\epsilon )=\frac{\sqrt{3/2}}{(\sigma \epsilon )^{1/4}}S(\gamma ).$$
(7)
The polynomial $`P_3`$ has one or three real roots. The collision of two real roots corresponds to the unstable equilibrium and its separatrix. It occurs when the discriminant

$$\mathrm{\Delta }=108(\sigma +\gamma ^2)$$

vanishes. $`\mathrm{\Delta }=0`$ is only possible for $`\sigma =1`$ and hence the critical parameters for which a double root occurs are given by $`\gamma =\pm 1`$, hence

$$9h^2=4\epsilon ^3,$$
(8)
which has a cusp at the origin. At the origin $`h=\epsilon =0`$ all three roots collide in the saddle-centre bifurcation. The discriminant (8) is shown in Fig. 2.

In the case of the saddle-centre bifurcation the bifurcation diagram is given by the discriminant of $`P_3`$. The bifurcation diagram divides the parameter plane $`(\epsilon ,h)`$ into two regions: one with 3 real roots and one with 1 real and two complex roots. The latter has positive $`\epsilon `$ everywhere, while the former is the wedge shaped region contained in the negative half-plane. For $`\epsilon <0`$ the phase portraits contain a pair of stable/unstable fixed points. The critical value of the energy of the stable fixed point is a local minimum given by the lower branch of the bifurcation diagram, while that of the unstable fixed point is a saddle given by the branch with positive $`h`$. For this range of energies there are 3 real roots. It will turn out that the upper branch of the bifurcation diagram corresponding to the unstable fixed point is crucial for the existence of vanishing twist. For $`\epsilon >0`$ the phase portrait is without fixed points. Even though the topology is trivial in this case we will now show that the rotation number has a maximum on a certain invariant torus containing points near the origin in the phase space. At this maximum of the rotation number the twist vanishes. The vanishing twist occurs at the vertical tangents of the contours of the rotation number shown in Fig. 3. The fact that the invariant curves are all unbounded for $`\epsilon >0`$ will be dealt with in section 4. For now observe that the integral (3) is finite, even though the invariant curves are unbounded in $`v`$.

The main feature of the level lines of the rotation number as shown in Fig. 3 is that it diverges when the unstable periodic orbit is approached. This occurs for the positive critical value of $`h`$ when $`\epsilon <0`$. Everywhere else the rotation number is a well defined, smooth and bounded function of $`h`$ and $`\epsilon `$. Accordingly the level lines “hug” the curve of critical values that correspond to the unstable orbit. Already from this property one can deduce the existence of a curve with vanishing twist using topological arguments. Here we proceed along the analytical route, because it will give us more detailed information. Note that the curve of critical values with negative $`h`$, see Fig. 2, does not appear in Fig. 3. One reason for this is that we chose to plot the rotation number for the (non-compact) invariant tori with motion between $`\mathrm{}`$ and the smallest real zero of $`P_3`$. These invariant tori do not contain critical points of the energy map, even though the corresponding energy might be a critical value. Accordingly the rotation number is smooth across this line of critical values. The critical point corresponding to the critical values is the stable fixed point at the local minimum of $`H`$. But even if we would plot the rotation number of the bounded invariant tori near the local minimum of the potential the picture is unchanged. The reason is that for a cubic elliptic curve the integrals of first kind over either one of the two real intervals (if they exist) are equal.

Since $`\mathrm{\Omega }=1/T`$ and $`2\pi T=J/h`$ the derivative of the rotation function is

$$\frac{\mathrm{\Omega }}{J}=\frac{\mathrm{\Omega }/h}{J/h}=\frac{T/h}{2\pi T^3}.$$

Hence the twist vanishes when

$$\frac{T}{h}=\frac{(3/2)^{3/2}}{\epsilon ^{7/4}}\frac{S}{\gamma }=0$$

and this is only possible for finite $`\epsilon `$ when

$$\frac{S}{\gamma }=\frac{1}{w^3}dz=0.$$

This complete elliptic integral can be written as a linear combination of Legendre’s standard integrals. In this way a condition for the vanishing of the twist is now obtained. The relevant case for this purpose is that of one real root, for which the phase portrait has no fixed point. The integrand $`w`$ is positive for $`u(\mathrm{},z_0)`$, where $`z_0`$ is the single root of $`P_3(z)`$. Let the factorized polynomial be given by

$$P_3(z)=(zz_0)P_2(z),P_2(z)=(z\zeta _1)^2+\zeta _2^2,$$
(9)
so that the complex roots are $`\zeta _1\pm i\zeta _2`$. Denote the distance between the real and complex roots by $`r`$, hence $`r^2=P_2(z_0)`$, so that the discriminant is given by $`\mathrm{\Delta }=4r^4\zeta _2^2`$. Legendre’s standard integral of the first kind $`K(k)`$ has differential

$$\omega _1=\frac{1}{\sqrt{P_3(z)}}\mathrm{d}z$$
(10)
up to a constant factor, where the modulus $`k`$ is given by

$$k^2=\frac{1}{2}\left(1+\frac{z_0\zeta _1}{r}\right)=\frac{1}{2}\left(1+\frac{\mathrm{sign}(\gamma )}{\sqrt{1+\alpha ^2}}\right),$$
(11)
and in the last equality the parameter $`\alpha =\zeta _2/(z_0\zeta _1)`$ has been introduced. In the second equality in addition $`\mathrm{sign}(z_0\zeta _1)=\mathrm{sign}(\gamma )`$ is used, which is true because $`z_0=\zeta _1`$ in (9) together with the vanishing of the quadratic coeffcient in (5) implies $`z_0=\zeta _1=0`$, and therefoe the polynomial has no constant term and $`h=\gamma =0`$.

A non-standard form of the differential of Legendre’s standard integral of second kind $`E(k)`$ is

$$\omega _2=\frac{P_2(z)\mathrm{d}z}{(z(z_0\pm r))^2}\omega _1,$$

up to the same constant factor as in (10). The differential $`\mathrm{d}z/w^3`$ we are interested in is of the second kind, and can therefore be written as a linear combination of $`\omega _1`$ and $`\omega _2`$ with constant coefficients, up to a total differential:

$$\frac{\mathrm{d}z}{w^3}=A\omega _1+B\omega _2+\mathrm{d}F$$

where $`F=Q_2(z)/((z(z_0\pm r))w)`$. Together with the undetermined coefficients of the quadratic polynomial $`Q_2`$ this gives a system of 5 linear equations for the 5 unknown coefficients. Solving these equations gives

$$A=\frac{2r}{\mathrm{\Delta }}((z_0\zeta _1)rr^2+4\zeta _2^2),B=\frac{4r}{\mathrm{\Delta }}(r^24\zeta _2^2)$$

Since the quadratic coefficient of $`P_3`$ is zero, the roots of $`P_3`$ add up to zero. Therefore the real parts satisfy $`z_0+2\zeta _1=0`$, hence

$$2k^2=13\zeta _1/r\text{and}r^2=9\zeta _1^2+\zeta _2^2.$$

With these equations the coefficients $`A`$ and $`B`$ can be expressed in terms of $`k`$ alone, up to the factor $`(r^2\mathrm{\Delta })^1`$. The condition of vanishing twist, $`S/\gamma =0`$, finally reads

$$(8k^49k^2+1)K(k)=(16k^416k^2+1)E(k),$$

where $`K(k)`$ and $`E(k)`$ stand for elliptic integrals of the first and the second kind, respectively.

For $`k=0`$ we have equality since both elliptic integrals equal $`\pi /2`$. The first derivatives of either side vanishes, but the second derivatives are $`35\pi /4`$ and $`65\pi /4`$, respectively, so that the left hand side is larger for small $`k`$. For $`k=1`$ the prefactor of $`K`$ vanishes, while that of $`E`$ gives 1 and $`E(1)=1`$. Hence for $`k1`$ the right hand side dominates. This proves that there exists a solution of this equation for $`k(0,1)`$. Numerically we find $`k_0^20.7097215`$. In order to calculate the corresponding $`\gamma _0`$ we observe that $`\alpha `$ as introduced in (11) is related to $`\gamma `$ by

$$\alpha =\frac{\zeta _2}{z_0\zeta _1}=\frac{\zeta _2}{3\zeta _1}=\frac{1}{\sqrt{3}}\frac{\mathrm{\Gamma }^2+\sigma }{\mathrm{\Gamma }^2\sigma },$$
(12)
where

$$\mathrm{\Gamma }=(\gamma +\sqrt{\sigma +\gamma ^2})^{1/3}.$$
(13)
Using (11) and $`k_0`$ the corresponding value of $`\alpha `$ is $`\alpha _02.164255`$, and from $`\alpha _0`$ using (12) we find $`\gamma _00.9152203`$. Therefore the curve of vanishing twist in the parameter plane occurs for positive $`\epsilon `$ when

$$3h=2\gamma _0\epsilon ^{3/2}.$$
(14)
Since $`\gamma _0<1`$ the curve of vanishing twist is bent downward as compared to the bifurcation curve (8) for $`\epsilon <0`$ and $`h>0`$. See Fig. 3 for a graph of this curve together with the numerically computed lines of constant rotation number. The lines of constant rotation number have vertical slope at their intersection with the critical curve, as must be the case.

The scaling of $`u`$ reduces the number of parameters to one. The essential parameter $`\gamma `$, in its dependence on $`h`$ and $`\epsilon `$, organizes the bifurcation. It allows to compute explicitly all the important characteristics of the twistless torus. Combining (11) and (12) shows that $`k`$ is a function of $`\gamma `$. The curves in the parameter plane that have the same value of $`\gamma `$ are given by (4). The most prominent ones are the curve of critical values $`\gamma =\pm 1`$, as shown in Fig. 2, and $`\gamma =\gamma _0`$, the curve of twistless tori, see Fig. 3. They all have the same shape of a semicubical parabola, except when $`\gamma =0`$, hence $`h=0`$, or $`\gamma =\pm \mathrm{}`$, hence $`\epsilon =0`$ with $`\pm h>0`$.

$$\begin{array}{cccc}k^2& \gamma & \alpha & \text{curve}\\ & & & \\ 1& +1& 0& 3h=2(\epsilon )^{3/2}\\ \frac{1}{4}(2+\sqrt{3})& +\mathrm{}& 1/\sqrt{3}& \epsilon =0,h>0\\ 0.709721497& 0.91522& 2.164255& 3h=\gamma _0(\epsilon )^{3/2}\\ \frac{1}{2}& 0& \mathrm{}& \epsilon =0,h>0\\ \frac{1}{4}(2\sqrt{3})& \mathrm{}& 1/\sqrt{3}& \epsilon =0,h<0\\ 0& 1& 0& 3h=2(\epsilon )^{3/2}\end{array}$$

Approaching the bifurcation point along the curve (4) gives $`k(\gamma )`$ in the limit. The function $`k(h,\epsilon )`$ is therefore not continuous at the origin. Moving on curves (4) in the parameter plane for any $`\gamma 1`$ the change in the period $`T(h,\epsilon )`$ given by (3) is elementary. From (7) it follows that $`T`$ is proportional to $`|\epsilon |^{1/4}`$, and the constant of proportionality is determined from (7) and (6). The divergence of the period upon approaching the bifurcation is therefore not caused by the elliptic integral, but merely by the algebraic dependence $`|\epsilon |^{1/4}`$. In this way the rotation number of the stable periodic orbit is given by

$$\mathrm{\Omega }(\gamma =1,\epsilon )=\frac{1}{\sqrt{2}\pi }(\epsilon )^{1/4}0.225079|\epsilon |^{1/4}.$$

The integral $`S(1)`$ contained in this expression can be easily calculated using residue calculus because for $`\gamma =1`$ the curve has a double root, $`P_3(z)=(z+2)(z1)^2`$. For general values of $`\gamma `$ and in particular for the twistless torus the elliptic integral $`S(\gamma )`$ (6) needs to be calculated. The single real root is given by

$$z_0=\mathrm{\Gamma }\frac{\sigma }{\mathrm{\Gamma }},\text{hence}z_0(\gamma _0)0.5535942.$$
(15)
Finally $`S(\gamma )`$ is obtained as

$$S(\gamma )=\frac{4K(k)}{\sqrt{3z_0\sqrt{1+\alpha ^2}/2}}.$$

In particular when $`\gamma =\gamma _0`$ the rotation number of the twistless torus is

$$\mathrm{\Omega }(\gamma =\gamma _0,\epsilon )=\frac{\sqrt{z_0(\gamma _0)\sqrt{1+\alpha _0^2}}}{4K(k_0)}\epsilon ^{1/4}0.1374244\epsilon ^{1/4}.$$
(16)
The constant of proportionality is close to $`\sqrt{3}/4\pi `$. Using the above value of $`z_0`$ the position of the rightmost point of the twistless torus in phase space is located at $`u_0=z_0\sqrt{\epsilon }`$. This means that for $`\epsilon =\epsilon _0`$ this point is on the same side as the stable periodic orbit was for $`\epsilon =\epsilon _0`$ before the bifurcation, but by a factor of $`1/z_0(\gamma _0)2`$ closer to the origin. See Fig. 4 for an illustration in the form of a standard bifurcation diagram showing position in phase space $`u`$ versus bifurcation parameter $`\epsilon `$.

## 3 Period-Doubling Bifurcation

The normal form of a Hamiltonian system with two degrees of freedom near a periodic orbit in $`1:2`$ resonance is

$$H(u,v)=\frac{v^2}{2}+D\frac{u^4}{4}+\epsilon u^2,$$
(17)
where $`D=\pm 1`$. As in the case of the saddle-centre bifurcation the variables $`u`$ and $`v`$ are canonically conjugate variables on a local transversal Poincaré section and $`\epsilon `$ is a parameter, typically corresponding to the energy of the original system. Contourplots of this Hamiltonian show the intersection of invariant tori with the Poincaré section transversal to the bifurcating orbit. For $`D=1`$ they are shown in Fig. 5, for $`D=1`$ in Fig. 6, for $`\epsilon <0`$, $`\epsilon =0`$, and $`\epsilon >0`$, respectively.

The critical points and critical values of the energy map $`H:^2`$ and their dependence on the parameter $`\epsilon `$ describe the structure of the bifurcation. The Hamiltonian has critical points at $`(0,0)`$ and at $`(u,v)=(\pm \sqrt{2\epsilon D},0)`$ with critical values $`h=0`$ and $`h=D\epsilon ^2`$. The origin is a local minimum for $`\epsilon >0`$, a saddle otherwise. The second critical point exists when $`\epsilon D<0`$, and is a local minimum for $`\epsilon <0`$, a saddle otherwise. The set of critical values, see Fig. 7, therefore is the union of the line $`h=0`$ with the half of the parabola $`h=D\epsilon ^2`$ for which $`D\epsilon <0`$. For $`D>0`$ the only unstable branch in the bifurcation diagram is $`h=0`$ for $`\epsilon <0`$. It divides the two regions of real motion. A third region $`\{h<\epsilon ^2\}\{\epsilon >0,\epsilon ^2<h<0\}`$ is not in the image of $`H`$, hence there is no real motion corresponding to $`(h,\epsilon )`$ from this region.

For $`D<0`$ the line of critical values $`h=0`$ again has a saddle as critical point when $`\epsilon <0`$. In addition the half parabola $`h=\epsilon ^2`$, $`\epsilon >0`$ also corresponds to a saddle of $`H`$ that is not at the origin. When $`D<0`$ the whole plane $`(h,\epsilon )`$ is in the image of $`H`$ considering all $`\epsilon `$. The critical values divide the plane into three regions with 0, 2, and 4 real roots.

The dynamics of the reduced one degree of freedom system is given by $`\dot{u}=v`$ and eliminating $`v`$ using the Hamiltonian gives a first order equation for $`u`$. Separation of variables then gives the period of motion in the reduced one degree of freedom system as the elliptic integral

$$T(h,\epsilon )=\frac{\mathrm{d}u}{\sqrt{2h\frac{1}{2}Du^42\epsilon u^2}}.$$

The number of parameters could be reduced by introducing the ratio $`h/\epsilon ^2`$, but for clarity we do not introduce this scaling. The period $`T`$ is an elliptic integral on the curve

$$:w^2=P_4(z)=P_2(z^2)=2h\frac{1}{2}Dz^42\epsilon z^2.$$
(18)
The discriminant of $`P_4`$ is

$$\mathrm{\Delta }=256Dh(Dh+\epsilon ^2)^2.$$

It vanishes for $`h=0`$ and $`Dh=\epsilon ^2`$, corresponding to the critical values $`h=0`$ and $`h=D\epsilon ^2`$ already found above. As usual the set of critical values is contained in the discriminant, however, the discriminant vanishes on a larger set. Namely the branch of the parabola with $`D\epsilon >0`$ is not part of the critical values. The three regions already found correspond to regions with 4, 2, and 0 real roots of $`P_4`$, respectively, see Fig. 7. The part of the discriminant that is not part of the critical values is dashed. The polynomial can be factored as

$$P_4(z)=\frac{1}{2}D(\xi _{}z^2)(\xi _+z^2),\frac{1}{2}D\xi _\pm =\epsilon \pm \sqrt{\epsilon ^2+Dh}.$$

Comparing coefficients gives $`\xi _{}\xi _+=4hD`$ and $`\xi _{}+\xi _+=4\epsilon D`$. The factorization of $`P_4`$ has real factors, i.e. $`\xi _\pm `$, with one exception. It occurs when the quadratic equation $`P_2(\xi )=0`$ has complex roots. The position of the roots in the regions of the parameter plane is as follows. The two real cases for $`D=1`$ are

$`+2:`$ $`h>0`$ $`\xi _{}<0<\xi _+`$ $`2[\sqrt{\xi _+},\sqrt{\xi _+}]`$
$`+4:`$ $`\epsilon <0,\epsilon ^2<h<0`$ $`0<\xi _{}<\xi _+`$ $`2[\sqrt{\xi _{}},\sqrt{\xi _+}]`$
In the first column a label is given containing the number of real roots. All 4 cases appear when $`D=1`$ for

$`2:`$ $`h<0`$ $`\xi _+<0<\xi _{}`$ $`2[\mathrm{},\sqrt{\xi _{}}]`$
$`4:`$ $`\epsilon >0,0<h<\epsilon ^2`$ $`0<\xi _+<\xi _{}`$ $`2[\mathrm{},\sqrt{\xi _{}}]`$
$`0i:`$ $`\epsilon <0,0<h<\epsilon ^2`$ $`\xi _+<\xi _{}<0`$ $`[\mathrm{},\mathrm{}]`$
$`0c:`$ $`h>\epsilon ^2`$ $`\xi _\pm \text{ complex}`$ $`[\mathrm{},\mathrm{}]`$
The last column gives the interval of real motion along which the period integral is taken. The region without real roots contains two parts separated by the branch of the parabola $`h=D\epsilon ^2`$ with $`D\epsilon >0`$. On this branch there occurs a collision of complex roots at $`z^2=2\epsilon D`$, and they move from the imaginary axis into the complex plane.

The intervals of integration and multiplication factors as given in the last column of the previous table can be read off from the phase diagrams, see Fig. 5, 6. In the case $`D=1`$ with 4 real roots there are different orbits for the same $`(h,\epsilon )`$. One is the compact orbit near the stable fixed point with extent $`2[\sqrt{\xi _+},\sqrt{\xi _+}]`$. The interval given above is for the non-compact orbit, which has half of this period.

In order to calculate the derivative of $`T`$ with respect to $`h`$ the integral is first written in standard form and then differentiated. This more traditional approach (as compared to the previous section) is preferable in this case because the roots of the quartic are easily written down.

Denote the ratio of the roots $`\xi _\pm `$ by

$$r=\frac{\xi _{}}{\xi _+}=\frac{\epsilon \sqrt{\epsilon ^2+Dh}}{\epsilon +\sqrt{\epsilon ^2+Dh}}.$$

Then the period $`T(h,\epsilon )`$ in the 6 cases is given by

$`+2,2:`$ $`T=`$ $`\left({\displaystyle \frac{8(12k^2)}{\epsilon }}\right)^{1/2}K(k),`$ $`k^2={\displaystyle \frac{1}{1r}},`$ (19)
$`+4,0i:`$ $`T=`$ $`\left({\displaystyle \frac{8(k^22)}{\epsilon }}\right)^{1/2}K(k),`$ $`k^2=1r,`$ (20)
$`4:`$ $`T=`$ $`\left({\displaystyle \frac{8(1+k^2)}{\epsilon }}\right)^{1/2}K(k),`$ $`k^2={\displaystyle \frac{1}{r}},`$ (21)
$`0c:`$ $`T=`$ $`\left({\displaystyle \frac{8(2k^21)}{\epsilon }}\right)^{1/2}K(k),`$ $`k^2={\displaystyle \frac{1}{2}}+{\displaystyle \frac{\epsilon }{2\sqrt{h}}}.`$ (22)
The level lines of the period $`T`$ (and hence the rotation number $`\mathrm{\Omega }=1/T`$) are shown in Fig. 8 for $`D=1`$ and in Fig. 9 for $`D=1`$. These numerically computed pictures show that there are no vertical tangents, hence the twist does not vanish. This is now proved by differentiating the period in Legendre normal form.

The twist vanishes when

$$0=\frac{T}{h}=\frac{T}{k}\frac{k}{h}.$$

The last factor $`k/h`$ does not vanish for $`k(0,1)`$. In the cases $`+2`$, $`2`$, and $`0c`$ it seems to vanish when $`k=1/2`$. But this implies that $`\epsilon =0`$ and the singularity cancels. Therefore it is enough to consider $`T/k=0`$. After removing common non-vanishing factors the conditions for vanishing twist are $`V_i(k)=0`$, $`i=1,2,3`$ where

$`+2,2,0c:`$ $`V_1(k)=`$ $`(12k^2)E(k)(1k^2)K(k),`$ (23)
$`+4,0i:`$ $`V_2(k)=`$ $`\left(1{\displaystyle \frac{k^2}{2}}\right)E(k)(1k^2)K(k),`$ (24)
$`4:`$ $`V_3(k)=`$ $`(1+k^2)E(k)(1k^2)K(k).`$ (25)
The three equalities $`V_i(k)=0`$, $`i=1,2,3`$ are never satisfied on the range $`k(0,1)`$. Obviously $`V_i(0)=0`$, while $`V_i(1)=1,1/2,2`$ for $`i=1,2,3`$, respectively. So we need to show that $`V_1`$ is negative and $`V_2`$, $`V_3`$ are positive for $`k(0,1)`$. Differentiating $`V_1`$ and $`V_2`$ gives the simple results

$`V_2^{}(k)=`$ $`{\displaystyle \frac{3}{2}}k(K(k)E(k))`$ (26)
$`V_3^{}(k)=`$ $`3kE(k)`$ (27)
In the last case $`V_3^{}(k)`$ is obviously non-positive, so that the twist in the case $`4`$ is a monotone funtion rising from 0, hence it is nonzero. Similarly also $`V_2^{}(k)`$ is non-positive, which follows from the well known inequality $`K(k)>E(k)`$. The first function $`V_1(k)`$ is negative on $`(0,1)`$, but not monotone. Rewriting it as

$$V_1(k)=(1k^2)(E(k)K(k))k^2E(k)<0$$

the inequality is clear because both terms are negative for $`k(0,1)`$. Therefore the twist never vanishes in a neighborhood of the period doubling bifurcation.

The relation between the cases $`D=\pm 1`$ is interesting. The main observation is that changing the sign of $`D`$ and $`\epsilon `$ inverts the overall sign of the potential $`V(x;D,\epsilon )=Dx^4/4+2\epsilon x^2`$. Therefore changing the signs of $`D`$, $`\epsilon `$, and $`h`$ leaves the roots of $`P_4`$ invariant. By this mapping the regions in parameter space with the same numbers of real roots are mapped into each other. The integration path does change in a less trivial way. The integration needs to be taken over the positive intervals of $`P_4`$ on the real axis. Changing the sings of $`D`$, $`\epsilon `$, and $`h`$ does change the sign of $`P_4`$. As a result the periods for $`D=1`$ are the $`\alpha `$ cycles of the elliptic curve, while those for $`D=1`$ are the $`\beta `$ cycles. Hence the period for the case $`4`$ can be obtained from that of $`+4`$ by replacing $`K(k)`$ by $`K^{}(k)=K(\sqrt{1k^2})`$. Similarly, the cases $`+2`$ and $`2`$ are mapped into each other.

The obvious symmetry in Figures 8 and 9 with respect to changing the sign of $`h`$ is related to the fact that changing the sign of $`D`$ and $`h`$ leaves the ratio $`r`$, and therefore also the corresponding modulus $`k^2`$, invariant. This means that the level lines in region $`+2`$ and $`2`$ can be obtained from each other by reflection on the $`h=0`$ axis. In a similar way the regions $`+4`$ and $`0i`$ have the same rotation number.

The period doubling could have been treated in a scaled version with only one essential parameter $`\delta =h/\epsilon ^2`$. However, the presentation seems more transparent in the unscaled version. Similar to the case of the saddle-centre bifurcation the modulus $`k^2`$ of the elliptic integral for the period is constant on the parabolas $`h=\delta \epsilon ^2`$. Again dependence on the parameters on these curves is simply algebraic, as before through $`|\epsilon |^{1/4}`$. The value of the modulus is not defined at the origin, but depends on the parabola on which it is approached. But in any case, there are no twistless tori near the origin in this bifurcation.

## 4 Universality

A major problem in our approach seems to be that the invariant twistless tori in the normal form are not compact. The normal form is obtained from an expansion near the bifurcation point, and is therefore local in phase space and local in the parameter. How can the rotation number of a non-local invariant torus be determined from this local normal form? To answer this question higher order terms need to be considered in the Hamiltonian. The Poincaré map near the bifurcation can be described by the time $`1`$ map of the Hamiltonian

$$\stackrel{~}{H}(u,v,t;\epsilon )=H(u,v;\epsilon )+G(u;\epsilon )+R(u,v,t;\epsilon ).$$

The remainder terms $`R`$ containing the periodic time dependence can be made arbitrarily small, but they cannot in general be removed while retaining a non-zero radius of convergence of the normal form. The first term $`H`$ is the normal form analysed in the previous chapters. We will concentrate on the saddle-centre case (2), but similar remarks apply to (17). The invariant tori near the origin for $`\epsilon >0`$ of $`H`$ are not compact. The higher order terms in $`G`$ can compactify them. The results about vanishing twist can be applied when $`G`$ does compactify these curves. However, the precise form of $`G`$ does not matter. Under the compactness assumption KAM theory can be applied to $`H+G`$, where $`R`$ is the perturbation. Many of the invariant tori of $`H+G`$ will persist. In particular a twistless invariant torus of $`H+G`$ will persist if it is sufficiently irrational. The curve (14) in the parameter plane therefore does not have invariant twistless tori in its perimage for every point, instead just for a cantorset of points. This is well understood. The main issue in this section is to understand the effect of adding $`G`$ to $`H`$. For small $`\epsilon `$ the essential contribution to the diverging period comes from the dynamics near the origin, while the dynamics on the invariant torus away from the origin has finite period. This is why the local normal form can give a statement about the dynamics on a non-local invariant torus near the bifurcation point. In particular we will now show that the curve of vanishing twist that was found to be emanating from the cusp of the saddle-centre bifurcation has a universal shape sufficiently close to the cusp singularity. In particular this means that the constant $`\gamma _0`$ that determines the shape of the curve of twistless tori in relation to the curve of critical values of the unstable orbits has the universal value $`\gamma _00.91522`$. Quantities derived from $`\gamma _0`$, like $`z_0`$ and the coefficients in (16), are accordingly also universal. Our calculation will show that the value of $`\gamma _0`$ is not influenced by the higher order terms $`G`$ in the Hamiltonian. Moreover, the following argument also shows that integrating the non-compact invariant torus up to $`v=\pm \mathrm{}`$ does not introduce and additional error.

Let the high-order truncated normal form Hamiltonian be

$$H(u,v;\epsilon )=\frac{1}{2}v^2+\frac{1}{3}u^3+\epsilon u+G(u,\epsilon )$$

where $`G(0,\epsilon )=0`$ and $`G(u,0)=0`$ is an analytic function. It is not necessary to assume that the higher order terms in $`G`$ depend on $`v`$ also, see e.g. , but even with such a dependence a slightly modified argument would work. As already explained, it is now assumed that $`G(u,\epsilon )`$ is such that the invariant curves near the origin for $`\epsilon >0`$ are compact.

Hamiltons equation for $`u`$ reads $`\dot{u}=v`$ as before. The period is obtained by solving $`H(u,v;\epsilon )=h`$ for $`v=v(u;h,\epsilon )`$ and then by integrating

$$T(h,\epsilon )=\frac{\mathrm{d}u}{v(u;h,\epsilon )}.$$

The idea is to split the integral into two parts; one part near the origin, where the main contribution originates, and the rest, which is called $`T_3`$. In addition the singular integral near the origin is split again into two parts, $`T_1`$ which has the same integrand as in the previous calculation with $`G=0`$, and a correction $`T_2`$ which contains the contribution from $`G`$. The integral $`T_1`$ will be the most singular, $`T_2`$ is mildly singular, and $`T_3`$ is regular. For sufficiently small $`\epsilon `$ then $`T_1`$ dominates and the previous result is recovered.

To achieve the splitting into $`T_1`$ and $`T_2`$ the multiplicative structure of the original inegrand has to be recovered. Solving $`H=h`$ for $`v`$ and inserting into $`\dot{u}`$ gives

$$\dot{u}^2=v(u;h,\epsilon )^2=Q_0(u;h,\epsilon )\widehat{Q}(u;h,\epsilon )$$

where $`\widehat{Q}=1+O(u)`$ when $`|u|C1`$ for some fixed constant $`C`$, and $`Q_0`$ is a polynomial of degree 3 in $`u`$ whose zeroes approach those of the original case with $`G=0`$ when $`\epsilon 0`$. Denote by $`\gamma _C`$ the part of $`H(u,v;\epsilon )=h`$ for which $`|u|<C`$, and $`\overline{\gamma }_C`$ the rest of the invariant torus. Then the period integral can be split as

$`T(h,\epsilon )=`$ $`{\displaystyle _{\gamma _C}}{\displaystyle \frac{\mathrm{d}u}{\sqrt{Q(u,\epsilon ,h)}}}+{\displaystyle _{\overline{\gamma }_C}}{\displaystyle \frac{\mathrm{d}u}{\sqrt{Q(u,\epsilon ,h)}}}`$ (28)
$`=`$ $`{\displaystyle _{\gamma _C}}{\displaystyle \frac{\mathrm{d}u}{\sqrt{Q_0(u,\epsilon ,h)}}}+{\displaystyle _{\gamma _C}}{\displaystyle \frac{\widehat{Q}(u,h,\epsilon )^{1/2}1}{\sqrt{Q_0(u,\epsilon ,h)}}}du+{\displaystyle _{\overline{\gamma }_C}}{\displaystyle \frac{\mathrm{d}u}{\sqrt{Q(u,\epsilon ,h)}}}`$ (29)
$`=`$ $`T_1+T_2+T_3.`$ (30)
Now $`T_3`$ is regular, and gives a bounded contribution, so we can ignore it. The integral $`T_1`$ is singular in the limit $`h,\epsilon 0`$. $`T_2`$ is less singular, and in particular bounded, because the numerator goes to zero in this limit. So we only need to show that $`T_1`$ approaches the complete integral when $`\epsilon 0`$. The integral is

$$T_1(h,\epsilon )=\frac{4\sqrt{3/2}}{r^{1/2}}F(\varphi ,k),\varphi =2\mathrm{arctan}\sqrt{\frac{u_0+C}{r}},$$

where $`u_0`$ is the real root, $`\zeta `$ is the complex roots of $`Q_0`$, and $`r`$ and $`k`$ are as before,

$$r=\sqrt{2u_0^2+|\zeta |^2},2k^2=1+\frac{3u_0}{2r}.$$

The integral $`T_1`$ does not diverge because of the modulus $`k^21`$, but because in the prefactor $`r0`$. In fact, we already observed that the modulus $`k^2(\sqrt{3}+2)/4<1`$ inside the first quadrant $`\epsilon 0,h>0`$. The behaviour of $`r`$ for small $`\epsilon `$ is obtained from

$$r^2=\frac{3}{2}u_0\sqrt{1+\alpha ^2}=D(\gamma )\epsilon ^{1/4}.$$

This shows that $`\varphi \pi `$, and the integral approaches the complete integral with the same modulus $`k(\gamma )`$ as before. Therefore even though $`T_1`$ is an incomplete integral, in the limit of small $`\epsilon `$ its value approaches that of the complete integral $`T(h,\epsilon )`$, which diverges. Since the other integrals $`T_2`$ and $`T_3`$ are finite the analysis obtained from the complete integral over the non-compact curve of $`H(u,v;\epsilon )`$ with $`G=0`$ therefore correctly describes the behaviour of the rotation number of the compact invariant curve obtained when $`G0`$.

## 5 Example: Hénon Map

The Hénon map in the area preserving case,

$$(x^{},y^{})=(yk+x^2,x),$$

illustrates the above. There are saddle-centre bifurcations in the Hénon map for which the invariant tori for $`\epsilon >0`$ are not compact, and hence the vanishing twist cannot be observed. This applies to the initial bifurcation at $`k=1`$ that creates the pair of fixed points, and also to many saddle-centre bifurcations that occur for $`k>4`$. However for $`k=1`$ a pair of period 3 orbits is created at the origin in a saddle-centre bifurcation of the third iterate of the map, and the corresponding invariant tori are compact. One of the period three points is located on the symmetry line $`x=y`$. In new coordinates $`(u,v)=(x,y+x)`$ the third iterate of the map expanded near the origin with parameter $`k=1\epsilon `$ is

$$(u^{},v^{})=((uv)(1+4v),v+\epsilon +(u^22uv+3v^2))+O(u^3,v^3,\sqrt{|\epsilon |}^3).$$

The (exact) location of the fixed points is $`(u,v)=(\pm \sqrt{\epsilon },0)`$, with trace of the Jacobian $`22\sqrt{\epsilon }+8\epsilon 8(\epsilon )^{3/2}`$ so that the (approximate) multiplier is $`\mu 1+\mathrm{i}\sqrt{2}(\epsilon )^{1/4}`$ and $`\overline{\mu }`$ for the fixed point at $`u=\sqrt{\epsilon }`$. The corresponding rotation number $`\omega `$ of the elliptic fixed point is obtained from $`\mu =\mathrm{exp}(2\pi \mathrm{i}\omega )`$, so that $`\omega =(\epsilon )^{1/4}/(\sqrt{2}\pi )`$. For small positive $`\epsilon `$ the Hénon map possesses compact invariant curves near the origin. A higher order normal form would give a $`G`$ that describes these invariant curves. The Hénon map is non-integrable, so that only sufficiently irrational invariant curves of the (high order) normal form will exist in the Hénon map. For the situation under consideration numerical experiments show that many of these invariant curves do exist.

Considering the third iterate of the Hénon map turns the pair of period three orbits into three pairs of fixed points with heteroclinic connections. In the normal form there is only one pair of fixed points, and the unstable fixed point has a homoclinic connection. To match the prediction in the case of more than one unstable fixed point the period and hence rotation number must be calculated for the heteroclinic connection. For $`\epsilon >0`$ there are invariant tori on which the dynamics becomes slow near the three points that are close to the three bifurcation points of the third iterate of the map. Hence the rotation number for the third iterate of the map between two successive such points gives the correct rotation number. The results together with the position of the twistless curve are shown in table 1. The values shown converge to the predicted values given in (15) and (16), however, fairly small $`\epsilon `$ are needed to see this. The convergence to the true value $`0.1374244\mathrm{}`$ occurs approximately as $`\frac{1}{2}\epsilon ^{5/4}`$.