text:
# Coarse-grained simulation of polymer translocation through an artificial nanopore

## 1 Introduction

In 1993, Branton and Deamer demonstrated that DNA could be threaded through a biological pore and that, by measuring the variations in the accompanying ionic current, information about DNA secondary structure could be obtained (Kasianowicz et al., 1996). Although these results were very promising, little progress has been made towards actual sequencing of DNA due to the limited voltage range that can be applied across a biological pore and the difficulty in measuring the current variations, because the shot noise is comparable to the expected signal (Bokhari et al., 2002). An alternative approach currently under investigation (Bokhari et al., 2002) consists in using the new capabilities in nanoscale semiconductor technology to create a more optimized pore geometry. A device consisting of a nanometer-sized pore in a thin silicon membrane with vertical transistors positioned along the wall of the pore may, in principle, be able to detect charge passing through the pore by measuring image charges in the transistor. By reducing the noise and maximizing the charge and current sensitivity, the proposed geometry may achieve DNA sequencing with order of magnitude improvements in speed and cost and minimal pre- and post-measurement procedures compared to current sequencing techniques. The expected scientific and therapeutic benefits will be dramatic, for example enabling rapid routine screening for treatable genetic disorders.

Although an empirical approach based on trial and error can lead to some progress in the fabrication of such devices, a fundamental physical understanding of the delicate interplay between the physico-chemical behavior of the macromolecule in solution and the interaction with the confining geometry induced by the solid surface of the pore is highly desirable. Molecular modeling is a powerful tool to investigate the collective behavior of systems resulting from the complex interactions between their individual components. Such an approach may allow us to gain insight into the phenomena of driven polymer translocation through a nanometer-scale pore from both an applied and a fundamental point of view. From an applied point of view, we can assess the feasibility of sensing and identifying individual nucleotides, to aid in the design of semiconductor-based devices and in the interpretation of the experimental measurements. In order to achieve this goal, atomistic modeling using high accuracy molecular models must be used to represent the properties of single stranded DNA, aqueous solution and substrate quantitatively. With this pre-requisite, it is possible to study the local structure, conformation, and short time dynamics of single-stranded DNA in confined geometries, and to investigate the translocation of ions and counterions across the channel as a function of the pore diameter and shape and of the magnitude of the potential difference applied. From a more fundamental point of view, it is of interest to study the mechanism of polymer translocation and the behavior of polyelectrolytes in confined geometries. This is of relevance to the process of pore entry, including studies of the role of counterion charge density gradients and of the role played by hydrodynamic effects. A better understanding of this process can guide the design of pore geometries that effectively steer DNA molecules toward the pore.

Polymer translocation involves traversal of a free energy barrier arising from several distinct physical effects, including loss of conformational entropy of the macromolecule upon traversal (an entropic barrier) (Muthukumar, 1999, 2001), mismatch in dielectric constant between the aqueous solution and the silicon pore (an electrostatic barrier) (Parsegian, 1969) and specific interactions of the macromolecule with the pore surface (an enthalpic barrier) (Zimm et al., 1992). In this study, we use a coarse-grained approach, in which groups of atoms are represented by single interaction sites, to capture many of the essential physical features of the system. By dramatically decreasing the number of degrees of freedom present, this model can engender improved basic understanding of field-driven polymer translocation over long length and time scales.

The remainder of this article is organized as follows: In section II, we define the interaction potential used in this study, the parametrization used to represent the phosphate backbone of DNA in salt–free aqueous solution, and the simulation methodology. In section III we present results obtained on the translocation process and on the influence of Bjerrum length and external electric field on the conformational characteristics of the polymer. Finally, conclusions drawn from this study as well as future research directions are discussed in section IV.

## 2 Methodology

### 2.1 Interaction potential

We study the translocation of a single negatively charged polymer chain designed to represent the phosphate backbone of a DNA fragment in an aqueous solution containing sodium counterions. Each molecule in this three-species mixture is represented as a set of spherical interaction sites. A coarse-grained approach in which groups of atoms are represented by single interaction sites is used in order to reduce the computational cost and to enable the study of the translocation process on long timescales.

We choose a simple interaction potential of the form: $`U(𝐫^N)=U_{\mathrm{str}}+U_{\mathrm{vdw}}+U_{\mathrm{coul}}+U_{\mathrm{ext}}+U_{\mathrm{wall}}`$. The first term decribes intramolecular valence interactions, with $`U_{\mathrm{str}}`$ representing bond stretching interactions. We chose to model the phosphate backbone as a flexible chain without bond angle bending or dihedral torsion interactions in our model. It must be pointed out that the chain is not fully flexible due to intramolecular van der Waals and Coulombic interactions between second and more distant neighbors. The last four terms describe nonbonded interactions, with $`U_{\mathrm{vdw}}`$, $`U_{\mathrm{coul}}`$, $`U_{\mathrm{ext}}`$ and $`U_{\mathrm{wall}}`$ representing van der Waals, Coulombic, external electric field, and molecular site – substrate interactions respectively. The pore is defined as a cylindrical channel (with axis along the $`z`$ direction) in a solid substrate, which is treated as a smooth van der Waals surface, and the external electric field $`E`$ driving polymer translocation is in the $`+z`$ direction. These contributions to the total potential are defined as:

$`U_{\mathrm{str}}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\mathrm{bonds}}{ij}}{}}{\displaystyle \frac{1}{2}}k_r(r_{ij}r_{\mathrm{eq}})^2`$ (1)
$`U_{\mathrm{vdw}}`$ $`=`$ $`{\displaystyle \underset{i<j}{}^{^{}}}4ϵ_{ij}\left[\left({\displaystyle \frac{\sigma _{ij}}{r_{ij}}}\right)^{12}\left({\displaystyle \frac{\sigma _{ij}}{r_{ij}}}\right)^6\right]`$ (2)
$`U_{\mathrm{coul}}`$ $`=`$ $`{\displaystyle \underset{i<j}{}^{^{}}}{\displaystyle \frac{q_iq_j}{r_{ij}}}`$ (3)
$`U_{\mathrm{ext}}`$ $`=`$ $`{\displaystyle \underset{i}{}}Eq_iz_i`$ (4)
$`U_{\mathrm{wall}}`$ $`=`$ $`{\displaystyle \underset{i}{}}4ϵ_{iw}\left[\left({\displaystyle \frac{\sigma _{iw}}{r_{iw}}}\right)^{12}\left({\displaystyle \frac{\sigma _{iw}}{r_{iw}}}\right)^6\right]`$ (5)
The primes on the sums in Eqs. (2) and (3) indicate that nearest–neighbor intramolecular nonbonded interactions are excluded from the sums. The distance $`r_{ij}`$ between sites $`i`$ and $`j`$ is defined as $`r_{ij}=|𝐫_{ij}|=|𝐫_j𝐫_i|`$ where $`𝐫_i`$ is the position of site $`i`$ and $`z_i`$ its vertical component.

The remaining step in building the molecular model is fixing the interaction parameters. These parameters need to be provided for every distinct combination of site types. In the next subsection, the site types and the interaction parameters used in the case of the coarse-grained model will be specified.

### 2.2 Parametrization of the coarse-grained model

An appropriate parametrization of the coarse-grained molecular model is obtained via an approximate mapping from the corresponding physical system. The DNA backbone is negatively charged and the concentration of positively–charged sodium ions is such that the system is neutral overall.

Water molecules are treated as point van der Waals (Lennard-Jones) interaction sites, with the Lennard-Jones (LJ) parameters determined from a ‘corresponding states’ mapping of the triple point of water onto that of the LJ system (Johnson et al., 1992) yielding $`\sigma _w=0.29`$ nm and $`ϵ_w=0.79`$ kcal/mol. These values are taken as the unit of length $`\sigma _0`$ and the unit of energy $`ϵ_0`$. The mass of a single water molecule is the unit of mass, $`m_0`$ = 18 a.u. Using these values, a unit of time $`t_0=\sqrt{m_0\sigma _{0}^{}{}_{}{}^{2}/ϵ_0}7\times 10^{10}s`$, a unit of charge $`q_0=\sqrt{ϵ_0\sigma _0}e/12`$, and a unit of electric field $`E_0=ϵ_0/(q_0\sigma _0)1400`$ V/$`\mu `$m are defined, where $`e`$ is the magnitude of the electron charge.

The sodium ions are treated as point interaction sites with the same LJ parameters as water and with a reduced charge $`q_{Na^+}=\alpha q_0`$. We have investigated systems with $`\alpha =1.2`$, $`2.4`$ and $`6`$ corresponding to a charge $`q_{Na^+}=e/10`$, $`e/5`$ and $`e/2`$. These values have been chosen in order to take into account dielectric screening effects in an approximate way. For simplicity, the mass of sodium ions is set equal to the mass of water molecules.

The phosphate backbone consists of alternating phosphate and deoxyribose groups. Each group is represented by a van der Waals point interaction site (Figure 1) with van der Waals parameters $`\sigma _{\mathrm{ph}}=2\sigma _0`$ and $`ϵ_{\mathrm{ph}}=ϵ_0`$. The phosphate sites carry reduced electric charge $`q_{\mathrm{ph}}=q_{\mathrm{Na}^+}`$ while the deoxyribose sites do not carry any electric charge. The bond stretching parameters are set to the following values: $`r_{\mathrm{eq}}=\sigma _0`$ and $`k_r=24`$ kcal/mol. The phosphate chain is quite flexible although not fully flexible due to the van der Waals and Coulombic interactions between second and more distant neighbors. The masses of both phosphate and deoxyribose sites are set to 5 $`m_0`$. Geometrical combination rules for $`\sigma `$ and $`ϵ`$ were applied to parametrize the LJ interactions between unlike species.

The pore is defined as a cylindrical channel in a solid substrate, which is treated as a smooth van der Waals surface, the atomic sites interacting with the nearest point on the surface of the pore or substrate. For simplicity, the LJ parameters of the wall sites are the same as those of water. The diameter of the pore is $`d_{\mathrm{pore}}=5\sigma _0`$ with a length $`l_{\mathrm{pore}}=2d_{\mathrm{pore}}`$. The diameter of the pore is $`2.5`$ times larger than that of the polymer in order to accomodate both the polymer chain and at most one solvation sphere. In Branton’s experiments, the alpha–hemolysin pore diameter has been estimated to be between $`18`$ to $`26\AA `$ (Kasianowicz et al., 1996; Meller et al., 2001) which is $`2`$ to $`2.5`$ times the diameter of a single-stranded DNA fragment ($`D_{ds}=20\AA `$, $`D_{ss}=10\AA `$ (Tinland et al., 1997)). The length of the pore was estimated to be $`l_{\mathrm{pore}}=52\AA `$, in good agreement with the value chosen in our simulation. Reduced variables defined in this subsection are used in the remainder of this article.

### 2.3 Simulation methodology

We carry out NVT (constant particle number, volume, and temperature) molecular dynamics simulations of a single charged polymer in solution, under an external field, with periodic boundary conditions. The unit cell contains a reservoir separated by an impenetrable substrate of thickness $`l_{\mathrm{pore}}`$ and surface area equal to the lateral dimension of the unit cell, and containing a cylindrical hole (nanopore) of diameter $`d_{\mathrm{pore}}`$ (Figure 2).

The chain consists of 40 sites, alternating between phosphate and deoxyribose groups, and the initial state is built such that the charged polymer is in the reservoir outside the pore cavity. The system contains, in addition to the polymer, 6000 water molecules and 20 Na<sup>+</sup> counterions with initial positions chosen randomly in the unit cell. The simulation was carried out at a reduced density $`\rho ^{}=0.85`$ and a reduced temperature $`T^{}=k_BT/ϵ_0=0.75`$, corresponding to the dense fluid phase of a LJ system. The applied external electric field of magnitude $`E^{}=E/E_0=0.5`$, corresponding to a field $`E=`$709 V/$`\mu `$m, is turned on after an initial equilibration run of duration $`t^{}=500`$. In the experiments conducted with the alpha - hemolysin nanopore, a voltage of 70-300 mV was applied across the system. The electrodes are located centimeters away from the pore, but the bulk of the voltage drop occurs across the biological membrane, leading to an estimated electric field of 15-60 V/$`\mu `$m, at least an order of magnitude smaller than the field used in the present work. Such a large field is used in the present study in order to ensure that field-driven translocation occurs on a time scale accessible to simulations.

A single-timestep velocity-Verlet MD integration scheme was used to integrate the equations of motion with a timestep $`\delta t=2.5\times 10^3`$, and the weak-coupling algorithm of Berendsen et al. (Berendsen et al., 1984) was used to maintain constant temperature. All LJ interactions were truncated at 2.5 $`\sigma _0`$ and no long-range corrections were applied. Long-range Coulomb interactions are evaluated to high accuracy using the particle-mesh Ewald (PME) method (Darden et al., 1993; Essmann et al., 1995). The PME technique, which is derived from the conventional Ewald method (Ewald, 1921; Toukmaji et al., 1996), makes use of the FFT to efficiently compute the long range part of the electrostatic interaction. In this study we used a highly optimized version of the PME with a relative accuracy of $`10^4`$.

We have carried out 3 simulations of systems with reduced charge magnitude $`q^{}=1.2`$, 2.4 and 6.0 for durations of $`4.2\times 10^6`$, $`3.0\times 10^6`$ and $`5.2\times 10^6`$ timesteps, respectively.

## 3 Results

The translocation process consists of three distinct stages: In the first stage, the polymer drifts through the reservoir, driven by the external electric field. In the second stage, the polymer is pushed against the lower substrate surface and behaves like a 2–dimensional chain wandering in search of the pore entrance. In the third stage, the polymer finds the entrance to the pore and translocation takes place.

Due to the relatively small system size, the time spent by the polymer in the second stage is small compared to that in the two other stages. We first study the electrophoresis of the charged polymer in the reservoir, then the kinetics of the translocation process and the possibility of identifying individual ions passing through the pore from the electrostatic signal produced.

### 3.1 Electrophoresis in the reservoir

Electrophoresis is a complex process dependent on the electrophoretic friction coefficient (Manning, 1981). This coefficient arises both from the macroion/solvent interactions, comprising the monomer/solvent interaction (the intrinsic friction) and from hydrodynamic interactions between monomer pairs, and from macroion/ion interactions. Since the macroion is negatively charged, it is surrounded by an oppositely charged liquid atmosphere (Fuoss et al., 1959). The force acting on the ions in the atmosphere is transmitted to the solvent, and therefore the atmospheric liquid behaves as a charged volume. Under an applied electric field, this charged liquid is subject to a volume force in the opposite direction from the drift velocity of the macroion. The macroion thus moves against a local hydrodynamic flow that slows its motion (relative to a hypothetical drift without the charged liquid). In addition, the electric field deforms the charged volume, increasing the atmospheric charge density at the end of the polymer opposite to the end which drift and decreasing ahead of the macroion. This effect induces an internal field (often referred in the literature as the relaxation field, or more accurately, the asymmetry field) opposing the external electric field, which also slows down the macroion motion.

Counterion condensation is in large part responsible for the structure and behavior of the liquid atmosphere and, as recognized in the late 60’s independently by Oosawa and Manning (Oosawa, 1971; Manning, 1978), plays a crucial role in the strong attractive interaction that acts between highly charged macroions such as DNA, leading to significant conformational changes. Counterion condensation results in a competition between Coulombic energy and entropy in minimizing the free energy of an aqueous solution containing mobile ions in the vicinity of an isolated macroion. For rodlike objects, whether the entropy or the Coulombic attraction dominates depends on the charge density of the macroion measured in terms of the Bjerrum length $`\lambda _B=q^2/k_BT`$, defined as the distance at which the Coulombic interaction between two charges is equal to the thermal energy $`k_BT`$. Counterion condensation occurs when the distance $`b`$ between charges on the macroion is small enough for the dimensionless parameter $`\mathrm{\Gamma }=\lambda _B/b`$ (the Manning parameter) to exceed unity. With the reduced units introduced in the previous section, a reduced Bjerrum length $`\lambda _B^{}=q_{}^{}{}_{}{}^{2}/T^{}`$ is defined. In the present model, the distance between 2 charges is $`b=2\sigma _0`$, leading to a Manning parameter $`\mathrm{\Gamma }=q_{}^{}{}_{}{}^{2}/2T^{}`$. Values for $`\lambda _B^{}`$ and $`\mathrm{\Gamma }`$ are listed in Table 1 for the systems with reduced charge $`q^{}=1.2`$, $`2.4`$, $`6.0`$. At the reduced temperature chosen in this study, Coulombic interactions dominate over the thermal energy for $`q^{}=6.0`$ and $`q^{}=2.4`$ systems, which thus should exhibit strong counterion condensation. For the $`q^{}=1.2`$ system, $`\mathrm{\Gamma }=0.96`$, which is close enough to unity to also exhibit counterion condensation. However, in this case we expect to have a much less strongly bound counterion - macroion complex than for $`q^{}=6.0`$ and $`q^{}=2.4`$.

It is difficult to define an unambiguous criterion for condensation. Indeed, as pointed by Oosawa (Oosawa, 1971), for a coiled polymer chain each charged group makes a sharp and deep potential hole at its position, each linear part of the chain makes a sharp and deep potential valley along its length and the coiled chain as a whole makes a potential trough in its apparent volume. Counterions located in these 3 regions are considered to be condensed, with the counterions at charged group holes being localized and the counterions in the other two regions being mobile. Therefore, an estimate of the degree of condensation depends on the criterion used to discriminate between free and condensed counterions. In the present study, a counterion is said to be condensed if it is within a distance $`r_c^{}`$ from an ion of the polymer.

Figure 3 shows representative configurations for the three systems both with and without applied external electric field $`E^{}`$. We first study the system behavior without an external applied field. The three systems exhibit counterion condensation but the nature of the resulting macroion–counterion complex varies with the magnitude of the charge. For $`q^{}=6.0`$, a large majority of the counterions are closely associated with the polymer chain, being mostly bound or in the narrow potential valley along the length of the chain, while for $`q^{}=1.2`$ the counterions are mostly located in the potential valley of the polymer chain. An intermediate situation occurs for $`q^{}=2.4`$. Figure 4 shows the radial pair correlation function $`g_{ic}(r^{})`$ between ions on the chain and counterions, for $`q^{}=1.2`$, 2.4 and 6.0 and $`E^{}=0`$. The common features exhibited are a main peak at $`r^{}1.5`$, corresponding to localized condensed counterions, and a shoulder around $`r^{}2.0`$, corresponding to mobile counterions in the narrow potential valley around the chain. The height of the main peak for $`q^{}=6.0`$ is roughly 7 times larger than for $`q^{}=1.2`$ and 4 times larger than for $`q^{}=2.4`$, confirming a stronger counterion localization for larger charge magnitude. For $`q^{}=1.2`$, Figure 5 shows the relative importance of longer range correlations with respect to short range correlations with the presence of secondary peaks around $`r2.3\sigma _0`$ and $`r3.5\sigma _0`$, indicating that a significant counterion fraction is located in the potential trough of the chain. Table 1 displays the average number of condensed counterions for the 3 systems for different values of $`r_c^{}`$ corresponding roughly to the bottom of the main peak in $`g_{ic}`$ ($`r_c^{}=1.8`$), the shoulder ($`r_c^{}=2.2`$) and a value taking into account the long range correlations ($`r_c^{}=6.0`$).

Under the applied external field, the counterions remain located very close to the macroion forming bound dipoles for $`q^{}=6.0`$. In addition, a conformational change occurs with the formation of a collapsed polymer structure. By contrast, for $`q^{}=1.2`$, the counterions are not condensed anymore (Figure 3). The force produced by the applied external field overcomes the attractive Coulombic force and strips the counterions from the chain. As a first approximation, the reduced critical charge at which the Coulombic force exactly balances the force produced by the external field is $`q_c^{}=E^{}\overline{r}^{\mathrm{\hspace{0.33em}2}}`$ where $`\overline{r}^{}`$ is the average reduced distance between a polymer ion and a counterion. Assuming a condensed state with all counterions bound to the polymer, the estimate $`\overline{r}^{}1.5`$$`2.0`$ leads to $`q_c^{}=1.125`$$`2.0`$. Thus, no counterions will be stripped from the chain for $`q^{}=6.0`$, while a large counterion fraction will be stripped away for $`q^{}=1.2`$. In the latter case, $`q_c^{}`$ is a lower limit since the counterions are not closely bound to the chain (see Figure 3). A Monte Carlo simulation study of the counterion condensation on a spherical macroion has shown the same phenomenon under a strong external field (Tanaka et al., 2002). With this estimate, the counterion fraction located in the potential trough of the chain will be stripped from the polymer for $`q^{}=2.4`$.

Figure 5 shows the ion–counterion pair correlation function when an external electric field is applied on the systems. For $`q^{}=1.2`$, the magnitude of the main peak for $`q^{}=1.2`$ is reduced by an order of magnitude, confirming that a large counterion fraction is stripped from the chain. The electric field has a much less drastic effect on the two other systems: A small fraction of counterions are stripped from the chain for $`q^{}=2.4`$ while for $`q^{}=6.0`$, the field enhances slightly the counterion binding to the polymer. Average values for the number of condensed counterions under external electric field are given in Table 1.

The main qualitative features exhibited during the time evolution are independent of the chosen value for $`r_c^{}`$. Figure 6 shows the time evolution of the number of condensed counterions for the three systems for $`r_c^{}=1.8`$ and $`r_c^{}=6.0`$. The $`r_c^{}=6.0`$ system does not exhibit any significant changes in the number of condensed counterions over the duration of the simulation. The counterions are mostly condensed on the chain, and the electric field slightly enhances this effect. By contrast, for the $`q^{}=1.2`$ system, a strong decrease in the number of condensed counterions occurs in an external field, the new equilibrium configuration being reached after a time $`t^{}2500`$. The $`q^{}=2.4`$ system exhibits an initial decrease of the number of counterions induced by the external electric field similar to that for $`q^{}=1.2`$ with a new equilibrium state also reached after a duration $`t^{}2500`$. This effect is weaker, however (see Table 1 and Figure 5) due to the fact that a larger fraction of counterions are condensed on the chain ($`q_c^{}<q^{}=2.4`$). The behavior during the translocation process is discussed in the next subsection.

We turn now to the study of the polymer conformations as a function of magnitude of the site charge $`q^{}`$ and electric field $`E^{}`$. A global measure of the polymer conformation is the mean square radius of gyration $`R_g^2`$, defined as:

$$R_g^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}<(𝐫_i𝐑)^2>$$
(6)
where $`𝐫_i`$ is the position of the $`i^{\mathrm{th}}`$ polymer site, $`N`$ is the number of sites in the chain, and $`𝐑`$ is the position of the center of mass of the chain. Average values of $`R_g^2`$ over a duration $`t^{}=2500`$ for the 3 systems with and without applied electric field are listed in Table 1. Figure 7 shows the time evolution of $`R_g^2`$ for $`q^{}=1.2`$, 2.4 and 6.0. For $`q^{}=1.2`$, the chain becomes more elongated under the action of the external field because the mobile counterions are stripped away by the field, leaving the bare negative ions on the polymer and leading to an electrostatic stiffening of the polymer. For $`q^{}=2.4`$, the first effect of the electric field, before translocation through the pore, is to strip away the fraction of counterions which are in the potential trough of the polymer, resulting in a loss of electrostatic screening of the negative charges on the polymer and an increase in $`R_g^2`$. Then the polymer chain is driven through the reservoir by the electric field and pushed against the lower substrate (Figure 3). In this regime, the polymer behaves like a 2-dimensional chain and its radius of gyration may be significantly reduced (Foo, 1997, 1998). Conformational changes during the translocation process are discussed in the next subsection.

For $`q^{}=6.0`$, the radius of gyration is of the same order of magnitude as for $`q^{}=1.2`$ and $`q^{}=2.4`$ when no external field is applied. Under an external electric field, no significant change in the magnitude of the radius of gyration is observed over a time interval $`t^{}5000`$, then a strong decrease in $`R_g^2`$ occurs, indicative of the tightly collapsed structure adopted by the chain (Figure 7). Additional studies are needed to determine whether the observed collapse is induced by the external field or is simply caused by strong counterion condensation. Biological processes and theoretical works tend to favor the latter hypothesis. The chain collapse is analogous to the packing of DNA into a cell (Bloomfield, 1996). This packing requires overcoming an enormous Coulombic barrier in a highly dilute aqueous solution containing a small concentration of polvalent cations. Simulations and theory have demonstrated short range attraction between two macroions modeled as charged cylinders (Gronbech–Jensen et al., 1997; Lyubartsev et al., 1995) or stiff polymers (Stevens, 1999). This attraction has been ascribed to correlated fluctuations of the counterions induced by counterion condensation (Ha et al., 1997; Manning et al., 1994). Self-attraction have been shown to occur also in flexible polyelectrolytes (Stevens et al., 1995; Schiessel et al., 1998; Brilliantov et al., 1998). Theoretical work using simple scaling arguments have shown that counterion condensation modifies the second virial coefficient of a polyelectrolyte due to the fact that at low enough temperature (or at high enough counterion valence) the counterions approach close enough to the macroion to form dipoles, leading to charge - dipole and dipole - dipole interactions (Schiessel et al., 1998). Similar ideas have been developed by Brilliantov et al. (Brilliantov et al., 1998), who also predict also a first order phase transition between a stretched polyelectrolyte and a strongly collapsed polyelectrolyte as a function of magnitude of electric charges due to counterion condensation.

### 3.2 Translocation process

The thermodynamics of electrophoresis in the presence of a narrow pore is complex and involves traversal of entropic and enthalpic barriers (Muthukumar, 1999, 2001; Ambjörnson et al., 2002; Slonkina et al., 2002; Lubensky et al., 1999; Sung et al., 1996; Boehm, 1999; Sebastian et al., 2000; Kumar et al., 2000; Lee et al., 2001). When the chain enters the pore, its conformational entropy is reduced, due to chain elongation. On the other hand, due to the small size of the pore aperture a significant fraction of the counterions are stripped from the chain and gain configurational entropy. The counterion unbinding also results in an increase in the Coulombic energy. In addition, the electrostatic energy of interaction of the macroion with the external field decreases as the negatively charged chain moves through the pore, and the electrostatic interaction energy of the free counterions also decreases. The situation is reversed when the polymer exits the pore. However, the gain in conformational entropy of the polymer upon exiting the pore is probably greater than the entropy loss associated with pore entry, since before pore entry the polymer is pushed against the lower substrate surface ($`z=z_0`$) and behaves like a conformationally restricted two-dimensional chain while it behaves like a three-dimensional chain upon exiting the pore. Detailed free energy calculations are needed to study the importance of these contributions in the translocation process. We focus here on the kinetics of the translocation process under a large electric field.

Translocation was only observed in the $`q^{}=2.4`$ simulation. As we have seen in the previous subsection, for $`q^{}=6.0`$, strong counterion condensation results in repulsive screened Coulomb interactions between the ions on the chain and self–attraction due to charge–dipole, dipole–dipole interactions and/or charge fluctuation along the chain. In the $`q^{}=6.0`$ system, the chain adopts a collapsed structure with an effective diameter larger than the pore diameter, effectively preventing translocation. There is, a priori, no reason for the $`q^{}=1.2`$ system not to translocate since the chain adopts an extended conformation. The simulation has to be performed over a long enough period of time to permit the polymer to find the pore entrance via diffusion along the substrate surface.

Due to the large external field, translocation of the polymer in the $`q^{}=2.4`$ system occurs after a short drift time spent in the reservoir. The polymer finds the entrance of the pore after a time interval $`\mathrm{\Delta }t^{}=4000`$ ($`\mathrm{\Delta }t=2.8`$ ns) after application of the field. The upper part of Figure 8 shows several stages in the translocation process. For clarity, only the polymer chain and the counterions are shown. The translocation starts when the first bead (head) enters the pore at $`z=z_0`$ and ends when the last bead (tail) exits the pore at $`z=z_0`$. The total translocation process occurs over a period of time $`\mathrm{\Delta }t^{}950`$, corresponding to 0.665 ns. Timescales reported for translocation of ions and water molecules are in the range $`10^610^9s`$, in reasonable agreement with our results. However this time is much shorter than the characteristic time reported for ssDNA translocation ($`10^6`$ s) (Kasianowicz et al., 1996; Meller et al., 2001; Henrickson et al., 2000) due to the larger electric field used in the present work.

We have decomposed the translocation process into three stages (Slonkina et al., 2002) in which the conformation of the polymer and its interactions with the pore are significantly different: the first stage corresponds to the translocation of the head from the entrance to the exit of the pore, the second stage corresponds to the exit of the head from the pore and the translocation of the tail from the reservoir to the entrance of the pore, and the third stage corresponds to the translocation of the tail from the entrance to the exit of the pore (see Figure 8). The duration of each stage of translocation is: $`\mathrm{\Delta }t_1^{}=t_{z=z_0}^{\mathrm{head}}t_{z=z_0}^{\mathrm{head}}=145.74`$ ($`\mathrm{\Delta }t_1=0.102`$ ns); $`\mathrm{\Delta }t_2^{}=t_{z=z_0}^{\mathrm{tail}}t_{z=z_0}^{\mathrm{head}}=424.75`$ ($`\mathrm{\Delta }t_2=0.297`$ ns) and $`\mathrm{\Delta }t_3^{}=t_{z=z_0}^{\mathrm{tail}}t_{z=z_0}^{\mathrm{tail}}=379`$ ($`\mathrm{\Delta }t_3=0.266`$ ns). The duration of the second and third stages are comparable, but are nearly three times as long as the first stage. These differences in time are probably due to variation in the magnitude of the total electric field $`E_{\mathrm{total}}^{}=E^{}+E_{\mathrm{int}}^{}`$ experienced by the polymer during the translocation, where $`E_{\mathrm{int}}^{}`$ is the internal field induced by the charge distribution.

The bottom part of Figure 8 shows the reduced electrostatic potential $`V^{}=V/V_0`$ ($`V_0=ϵ_0/q_0`$) due to ions and counterions in the $`(xz)`$ plane passing through the middle of the pore. We clearly see an inversion in the sign of the potential difference across the pore during the second stage of polymer translocation due to a change in the relative charge density on the two substrate surfaces. This change leads to an inversion in the direction of the internal field $`E_{\mathrm{int}}^{}`$, which in turn slows down the translocation process. The upper and lower substrate surfaces behave like a capacitor, with an excess of positive charges on the upper surface – the free counterions trying to translocate and pushed against the surface by the external electric field – and an excess of negative charges on the lower surface – the polymer chain wandering at the pore entrance. This situation is reversed when enough negatively charged monomers have translocated, analogous to the discharge of a capacitor.

Figure 9 shows quantitatively the time evolution of the reduced internal electric field $`E_{\mathrm{int}}^{}`$, defined as $`E_{\mathrm{int}}^{}=\mathrm{\Delta }V_{\mathrm{pore}}^{}/l_{\mathrm{pore}}^{}`$, across the pore during the translocation event, with $`\mathrm{\Delta }V_{\mathrm{pore}}^{}`$ being the average potential difference across the pore defined as $`\mathrm{\Delta }V_{\mathrm{pore}}^{}=1/4_{i=1}^4(V^{}(x_i,y_i,z_0)V^{}(x_i,y_i,z_0))`$ with ($`x_i=0`$; $`y_i=\pm r_0`$) for $`i=1,2`$ and ($`x_i=\pm r_0`$; $`y_i=0`$) for $`i=3,4`$ (see Figure 2). It confirms that the internal field decreases in magnitude and changes sign when the middle of the chain crosses the middle of the pore (i.e., in the second stage) resulting in a smaller magnitude of total electric field $`E_{\mathrm{tot}}^{}`$ and a slower translocation process.

We turn now to the identification of single translocation events involving either a counterion or charged sites on the polymer. In order to mimic the signal recorded at the gate of a transistor located on the inner cylindrical surface of the pore, we have computed the electrostatic potential $`V^{}`$ at 4 equidistant positions ($`x=0`$; $`y=\pm r_0`$) and ($`x=\pm r_0`$; $`y=0`$) on the surface of the mid–plane ($`z=0`$) of the pore. We have simultaneously recorded events corresponding to the passage of charged sites (either counterions or charged polymer beads) through the $`z=0`$ plane to investigate the correlations between signal and site positions. Figures 10 and 11 show the time evolution of $`V^{}`$ averaged over the four sensors as well as ion passage events during the translocation process. A positive spike records the passage of a charge in the same direction as the external electric field while a negative spike records the passage of a charge in the opposite direction. No counterions are pulled into the pore with the polymer so we do not have any counterions translocating against the electric field direction during polymer translocation. A strong correlation is found between counterion positions and the signal (sharp peaks in $`V^{}`$). It is also fairly easy to detect when the polymer translocation takes place (and its approximate duration), but the discrimination of individual charges on the polymer chain is more problematic due mainly to the small separation between neighboring charged sites on the chain, fluctuations in chain position, and the simultaneous translocation of counterions. Further analysis or filtering of the signal is required to achieve single-ion discrimination even with this simple model.

We notice that a significant larger number of counterions find the entrance of the pore during the third stage of polymer translocation. During the translocation process, the number of condensed counterions on the chain increases significantly (Figure 6). Due to the applied electric field, an excess counterion density is located near the upper substrate surface and a significant fraction of them – the ones near the aperture – condense back onto the polymer as soon as the head of the chain exits the nanopore (see Figure 8). The condensed counterions near the tail of the polymer are relatively mobile and have a significant probability of being stripped from the chain due to the external field and friction with the water molecules. The translocation of these counterions is facilitated by the polymer, which acts as a guide to the pore entrance.

In addition, the head of the polymer exiting the nanopore during the translocation process experiences a stronger electric field than inside the pore (see above). The polymer tail is still pushed by the field against the lower substrate surface creating an anchoring point on the substrate (see Figure 8) while the head of the chain is stretched by the electric field at the pore exit. This effect results in an increase in $`R_g^2`$ (Figure 7). After the head has exited the pore, counterion condensation back on the head, screening the Coulombic interactions and reducing the stretching force, as well as a weaker anchoring at the entrance of the pore due to the presence of fewer monomers, are responsible of the decrease in polymer extension, the chain behaving like a rubber band. As the translocation proceeds, the tail enters the pore and the release of the anchoring point induces a significant contraction of the chain. After the translocation, the chain behaves as a three-dimensional chain, resulting in an increase in $`R_g^2`$.

## 4 Conclusions

Using a coarse-grained bead–spring model, we have studied the electrophoretic translocation of a charged polymer through a nanopore for various charge magnitudes under the application of an external electric field. Three regimes can be identified, corresponding to drift of the polymer in the reservoir, diffusion of the polymer along the substrate surface in search of the pore entrance and polymer translocation through the pore. The three systems studied, with reduced charges, $`q^{}=1.2`$, 2.4 and 6.0, respectively, exhibit counterion condensation when no electric field is applied, but the $`q^{}=1.2`$ system presents a much less strongly bound macroion–counterion complex than the others two. Counterion condensation leads to significant changes in polymer confomation, with a rather extended chain conformation for $`q^{}=1.2`$ and a collapsed chain conformation for $`q^{}=6.0`$, due to self-attraction induced by charge–dipole/dipole–dipole interactions and/or counterion fluctuations along the chain. An external electric field can strip the counterions from the chain, depending on the relative magnitude of the field and the charge of the ions/counterions, inducing significant changes in polymer conformation. Due to the strong external electric field used, translocation occurs quickly for the $`q^{}=2.4`$ system and simulation over longer time scales should lead to the translocation of the $`q^{}=1.2`$ system. By contrast, translocation in the $`q^{}=6`$ system is inhibited due to the collapsed conformation adopted by the chain. The kinetics of the polymer translocation are not constant over the whole process and the chain translocates more slowly when its midpoint crosses the mid-plane of the pore due to an inversion in the direction of the internal electrostatic field produced by the time-varying distribution of charges. The upper and lower substrate surfaces behave qualitatively like a capacitor which undergoes a discharge. By recording the electrostatic potential we were able to identify the passage of individual counterions unambiguously. However, it is difficult to identify translocation events involving single charged sites on the polymer due mainly to the small distance between neighboring ions, fluctuations in polymer position, and simultaneous counterion translocation.

This work is mostly exploratory and many questions remain unanswered. Future work will focus on both the translocation phenomenon and on bulk electrophoresis. Using the simple model presented in this work, we will study the electrophoresis process in detail, including the effect of external electric field on counterion condensation as a function of the ion charge magnitude and macroion concentration. Special emphasis will be given to the modification of macroion conformation induced by the external field. In the framework of the simple model introduced in this study, we will study the translocation process statistically by steering the polymer through the pore and performing free energy computations to assess the relative importance of entropic and enthalpic contributions to the free energy of pore traversal. Since one of our aims is to study the feasibility of a biosensor for fast DNA sequencing we have initiated large-scale atomistic simulations of ssDNA in aqueous solution to study the electrostatic field produced by different nucleotides in the vicinity of the macromolecule, to explore the possibility of discriminating different bases based on their electrostatic signature.


format: text

text:
# Rotation curves in spiral galaxies due to intergalactic accretion flows onto the disc

WITHDRAWN:

The above mentioned paper is withdrawn. I have found that the approximation used in eq. (32) is valid for very high values of $`R`$, but not in the range 20-40 kpc, where the rotation curves were examined. The asymptotic regime of eq. (30) into eq. (31) is achieved in this range of galactocentric distances, but the asymptotic regime of eq. (26) into eq. (32) is not achieved in the same range of distances. Therefore, subsection 5.1 of the paper was not correct. New calculations of the amplitude of the U-warp produced by the proposed mechanism give much higher values, which are incompatible with the observed cases of U-warped galaxies. This fact convinced me that the proposed mechanism does not work to explain the flat rotation curves, unless some new mechanism is found to explain why the disks are stickier than derived from the calculation of the gravitational forces, for which we have not any evidence yet.

Acknowledgements: Thanks are given to F. J. Sánchez-Salcedo (UNAM, México) who has read the paper and found the inconsistency.


text:
# A route to computational chaos revisited: noninvertibility and the breakup of an invariant circle

## 1 Introduction

In an insightful 1989 paper entitled “Computational Chaos: a prelude to computational instability” E. N. Lorenz  reported on a one-parameter computational study of the dynamics of the noninvertible map

$`L_{(a,\tau )}:R^2R^2=\{\begin{array}{c}x_{n+1}=(1+a\tau )x_n\tau x_ny_n\hfill \\ y_{n+1}=(1\tau )y_n+\tau x_n^2,\hfill \end{array}`$ (3)
which arises from a simple forward Euler integration scheme ($`\tau `$ being the time step of the integration) of the two coupled nonlinear ODEs

$`dx/dt`$ $`=`$ $`axxy`$
$`dy/dt`$ $`=`$ $`y+x^2.`$ (4)
These ODEs are obtained by starting with the familiar Lorenz system : $`\dot{x}=\sigma (xy),\dot{y}=xz+\rho xy,\dot{z}=xy\beta z`$, letting $`\sigma \mathrm{}`$, and rescaling the variables and the remaining parameters. Lorenz fixed the value of $`a`$ at $`0.36`$, and varied the time step $`\tau `$. Although the corresponding differential equations (eq. (4)) exhibit an attracting equilibrium point, computer simulations indicated that the discrete approximation (eq. (3)) progressed from exhibiting an attracting fixed point, to an attracting invariant circle (IC), to a chaotic attractor (termed computational chaos) as the time step $`\tau `$ was gradually increased. Sequences of bifurcations similar to those described by Lorenz have been observed for other noninvertible maps of the plane , suggesting that this might be a universal “noninvertible route to chaos.” Some of the bifurcations along the route — local changes of stability, homoclinic and heteroclinic tangencies, crises caused by an attractor interacting with its basin boundary — are also observed in families of invertible maps. New bifurcations, however, unique to noninvertible families, as well as “invertible bifurcations with noninvertible complications,” are also observed along the route.

In this paper we perform a more detailed numerical study of Lorenz’s family, but focus primarily on a narrow range of parameters which includes the breakup of the invariant circle and the identification of $`\tau _{CD}`$ (the greatest lower bound on the values of the parameter $`\tau `$ for which the corresponding maps have “chaotic dynamics”), and $`\tau _{chaos}`$, the greatest lower bound on the values of the parameter $`\tau `$ for which the corresponding maps exhibit a “chaotic attractor.” Lorenz was interested in $`\tau _{chaos}`$ (labelled $`\tau _b`$ in ) as an indication of how poorly the Euler map $`L_{(0.36,\tau )}`$ approximated the original differential equation. He observed that, in the transition from smooth IC to chaotic attractor, the process began with the circle developing features with increasingly high curvature. There appeared to be a critical $`\tau `$ value at which the “IC developed cusps.” After this value, the “chaotic attractor” suggested by computer simulations appeared to have loops and a Cantor-like structure. It certainly was no longer a topological circle. He claimed (correctly) that if this scenario did in fact occur, then this critical $`\tau `$ value had to be an upper bound on $`\tau _{chaos}`$.

Our numerical investigations suggest that a smooth IC does not persist all the way to a “cusp” parameter. Instead, the IC is destroyed in a heteroclinic tangency initiating the crossing of a branch of the stable manifold and a branch of the unstable manifold of a periodic (here period-$`37`$) saddle point. (As for invertible maps, chaotic dynamics are apparently present during this crossing, so this tangency is an upper bound on $`\tau _{CD}`$.) It is the unstable manifold that subsequently develops cusps at a critical parameter value we label $`\tau _{cusp}`$, and loops beyond that value. A second manifold crossing results in the apparent appearance of the chaotic attractor, with loops inherited from the unstable manifold. We conclude that $`\tau _{cusp}`$ is not a single isolated bifurcation value separating smooth IC attractors from chaotic attractors.

Rather, $`\tau _{cusp}`$ is strictly above $`\tau _{CD}`$ and strictly below $`\tau _{chaos}`$. All three parameter values are part of the transition from a smooth IC attractor to a chaotic attractor. Because a cusp on an IC with an irrational rotation number would necessarily force a dense set of cusp points on the IC, it would make the existence of the IC itself unlikely. We therefore expect the existence of an entire transition interval of parameters, as opposed to a single bifurcation parameter, to be the generic scenario for the transition between smooth IC and chaotic attractor in the presence of noninvertibility.

Mechanisms for the breakup of ICs are relatively well established for invertible maps of the plane (see for example ). In the invertible scenario which appears to us to most closely match the noninvertible scenario of this paper, an IC is born in a Hopf (also called Neimark-Sacker) bifurcation, grows in size, coexists with a periodic orbit (after a saddle-node birth of the periodic orbit off the IC), is destroyed in a first global manifold crossing (also referred to in the literature as a crisis  where the IC collides with its basin boundary), and is reconstituted after a second global manifold crossing. The periodic orbit which persists through the destruction and reconstitution of the IC typically switches from outside (inside) the IC before it is destroyed to inside (outside) the IC after it is reconstituted. Chaotic invariant sets necessarily exist only during the manifold crossings. In contrast, the breakup of the IC in the noninvertible scenario is part of a transition to a chaotic attractor. In particular, the “chaotic attractor with loops,” which appears in this noninvertible case, is not a feature which appears in invertible maps. In fact, cusps and loops on iterates of smoothly embedded curves are not possible for smooth invertible maps, but are common features on (invariant) curves of noninvertible maps . Further, bifurcations such as the manifold crossings (crises) can have additional complications due to the presence of noninvertibility. For example, the unstable “manifolds” involved in the global crossings may have self intersections and cusps. Or stable and unstable “manifolds” may cross transversely at one homoclinic point, but fail to preserve transversality at other points along the homoclinic orbit.

In any case, the transition mechanism for the maps in eq. (3) is different from the invertible transitions, and the noninvertible nature of the map plays an important role: the IC — for parameter values for which it exists — interacts with (actually crosses) the locus on which the Jacobian of the linearized map becomes singular. This locus (and its images and preimages) is crucial in organizing the dynamics of noninvertible maps; it is termed “critical curve” in the pertinent literature  and constitutes the generalization (in two dimensions) of the critical point in unimodal maps of the interval . In particular, as explained in , the angle of intersection between a saddle unstable manifold and the critical curve underpins the transition of the local image of the manifold from being smooth and injective to nonsmooth and injective (with a cusp) to smooth but not injective (with loops).

As for invertible Hopf bifurcations (even in the noninvertible setting the dynamics near a Hopf bifurcation are necessarily locally invertible), an understanding of Lorenz’s $`a=0.36`$ one-parameter family is possible only by viewing it in the two-parameter Arnold horn context . We do this with $`L_{(a,\tau )}`$, using the second parameter ($`a`$) already provided in eq. (3). This leads to an examination of the internal Arnold resonance horn structure, to be compared to and contrasted with the invertible case. The interiors of the resonance horns have features which vary significantly from those in the horns of invertible maps (see preliminary results in ). Although our understanding of this internal structure is far from complete, examination of the partial picture still provides insight into the bifurcations observed in the Lorenz one-parameter cut. Developing computational tools to further investigate the internal structure is part of our ongoing research.

In this paper the observations of the transition in  are briefly summarized in Sec. 2, and then revisited in detail in Sec. 3. In Sec. 4, we place our one-parameter cut in context by examining the “resonance horns” in the larger two-parameter space. We discuss related noninvertible issues, including computational challenges, in Sec. 5, and state final conclusions in Sec. 6.

## 2 Observations

Figure 1 briefly summarizes the relevant observations in . For $`a`$ fixed at $`0.36`$, the differential equation (eq. (4)) has an attracting equilibrium point at $`(x,y)=(0.6,0.36)`$. Its basin of attraction is the right half plane. There is a symmetric attracting equilibrium for the left half plane at $`(0.6,0.36)`$. We will deal only with the right half plane attractors in this paper. The corresponding maps $`L_{(0.36,\tau )}`$ all have a fixed point at $`(x,y)=(0.6,0.36)`$. In a one-parameter cut with respect to $`\tau `$, the fixed point is stable for small $`\tau `$, but loses its stability via Hopf bifurcation at $`\tau =1.38889`$.

The resulting attractor is initially a smooth IC (fig. 1a, $`\tau =1.55`$) but then develops progressively sharper features (fig. 1b, $`\tau =1.775`$, where the circle is still smooth even though it appears to have cusps), then self-intersects, suggesting chaos (fig. 1c, $`\tau =1.785`$), and eventually exhibits clear signs of chaotic behavior (fig. 1d, $`\tau =1.91`$). Lorenz supports his claim of chaos with the computation of a positive Lyapunov exponent for $`\tau =1.91`$ . The last two figures are unusual in that the attractor apparently shows self-intersections, a clear sign of noninvertibility. All four figures were created by following a single orbit, after dropping an initial transient, from an initial condition near the fixed point. We use the term attractor loosely to denote the object generated by the computer via such a simulation. More formal definitions appear in the next section.

Lorenz proposed the following scenario as the simplest possible transition of the attractor from a smooth IC to a chaotic attractor. As a parameter is increased ($`\tau `$ in this example) and the IC crosses successively farther over the critical curve (called “$`J_0`$” below), causing it to develop successively sharper features. Simultaneously, its rotation number changes (increasing with $`\tau `$ in this case). By restricting to $`\tau `$ values for which the IC is quasiperiodic (and avoiding periodic lockings, defined below in Sec. 3.1), one obtains a sequence of $`\tau `$ values, each corresponding to a map with a smooth quasiperiodic IC, and limiting to $`\tau _{chaos}`$, the lower limit on chaotic behavior. At $`\tau _{chaos}`$, the attractor, which may or may not still be a topological circle, would develop cusps. Beyond $`\tau _{chaos}`$, the attractor would be chaotic, or at least have parameter values accumulating to $`\tau _{chaos}`$ from above for which the corresponding maps would exhibit chaos. We will now examine the implications of such a scenario.

## 3 Observations Revisited

By performing a more detailed numerical study, we show below that, while Lorenz’s scenario is the correct basic description of the transition, the IC actually breaks apart before the “cusp parameter” value is reached. The transition from smooth IC to chaotic attractor thus occurs over a range of parameter values, rather than at a single critical parameter value. We emphasize that Lorenz did not claim that his scenario did happen, only that there could be no simpler scenario. Our study shows that the proposed “simple” scenario does not happen in this particular family, nor should it be expected in the general transition from smooth invariant circle to chaotic attractor in the presence of noninvertibility.

Before beginning the description of our numerical investigation, we recall some terminology and preliminary results.

### 3.1 Noninvertible preliminaries

Let $`f:\mathrm{}^2\mathrm{}^2`$ be a smooth ($`C^{\mathrm{}}`$) map. Let $`x`$ be a period-$`q`$ point for $`f`$. As for invertible maps, we define the stable manifold of $`x`$ to be $`W^s(x)=\{y\mathrm{}^2:f^{kq}(y)x\mathrm{as}k\mathrm{}\}`$. Since the inverse map is not necessarily uniquely defined, we define the unstable manifold of $`x`$ as $`W^u(x)=\{y\mathrm{}^2:\mathrm{there}\mathrm{exists}\mathrm{a}\mathrm{biinfinite}\mathrm{orbit}\{y_j\}\mathrm{with}y_{i+1}=f(y_i),y_0=y\mathrm{and}y_{kq}x\mathrm{as}k\mathrm{}\}`$. The use of the term “manifold” is an abuse of terminology for noninvertible maps since, while the local stable and unstable manifolds are guaranteed to be true manifolds, the global manifolds are not . (Some authors address this problem by using “set” instead of “manifold” .) By iteration of the local smooth unstable manifold, the global unstable manifold can be computed as for invertible maps. If the manifold (or any curve) ever crosses the critical set ($`J_0`$, defined below) tangent to the zero eigenvector of the map at the crossing point, its image can have a cusp . More degenerate cases are also possible, but not considered here. The global unstable manifold can also have self intersections. Thus, it is globally a smoothly immersed submanifold, with smoothness violated only at forward images of these special crossing points. This is discussed further below. The global stable manifold, due to multiple preimages, can be disconnected, fail to be smooth, and even increase in dimension . The multiple preimages lead to interesting basins of attraction (see Section 5 and ).

A set $`S\mathrm{}^2`$ is invariant if $`f(S)=S`$. When $`S`$ is invariant and a topological circle we call it an invariant circle (IC). An invariant circle on which there exists a periodic orbit is called a frequency locked circle or circle in resonance. (The terminology is borrowed from return maps of periodically forced systems.) Typically we see a single attracting and a single repelling periodic orbit on a frequency locked circle, but multiple orbits are certainly possible. A periodic orbit is also called a periodic locking or a locked solution, although we usually reserve the use of periodic locking to refer to a periodic orbit on an invariant circle.

We have already indicated above that we are using the term attractor informally to indicate whatever is suggested via computer simulations. More formally, we define a compact set $`K`$ to be an attractor block if $`f(K)K^{}`$ (the interior of $`K`$). Then we define $`\mathrm{\Lambda }`$, the attracting set associated with attractor block $`K`$ by $`\mathrm{\Lambda }=_{n=1}^{\mathrm{}}f^n(K)`$ . In the cases we consider, including for $`\tau `$ in the range depicted in fig. 1, we take the attractor block $`K`$ to be an annulus containing the “attractor” in each picture: inner radius a small circle around the repelling fixed point $`O`$, and outer radius a (topological) circle which is big enough to contain the attractor, but still in the right half plane. The attracting set will always be taken to be associated with this annulus.

Any invariant circles or periodic points in the attractor block will also be contained in the attracting set $`\mathrm{\Lambda }`$. So will the unstable manifolds of the periodic points. The attractors displayed by computer simulation are generally (an approximation of) some subset of $`\mathrm{\Lambda }`$. If $`\mathrm{\Lambda }`$ exhibits at least some approximate recurrence property, such as a dense orbit (for example, a quasiperiodic IC), or a numerical orbit which appears dense (for example, a locked circle with a high period locking), then the attractor might be a good approximation of all of $`\mathrm{\Lambda }`$. In other cases, such as a locked circle with a low period locking, simulation would only display the sinks. Computation of the unstable manifolds of the saddles (crucial to our numerical investigation) often gives a better picture of $`\mathrm{\Lambda }`$ than does a simulation. Note that a locked circle is still an invariant circle, even though the whole circle is not the attractor. The circle is smooth by the unstable manifold theorem  through the saddles, but as with invertible maps, its smoothness where two branches of the unstable manifold meet at the sinks depends on the eigenvalues . Unlike invertible maps, as indicated above, the smoothness of an IC can also be violated by the development of cusps.

We now define the critical curve $`J_0`$ (termed “curve C” in  and often $`LC_1`$ for Ligne Critique in the recent noninvertible literature ) as the locus in phase space where the linearization of the map becomes singular. We call its image $`J_1`$ (“curve D” and $`LC`$, respectively). The two sets $`J_0`$ and $`J_1`$, as well as additional images and preimages, are known to be a key to understanding the dynamics of a noninvertible map. Consider, for example, Fig. 2a showing an enlargement of the IC, $`\mathrm{\Gamma }`$ for $`\tau =1.775`$, with three important additions: $`J_0`$, $`J_1`$, and $`\mathrm{\Gamma }_2^1`$, one additional first-rank preimage of the invariant circle $`\mathrm{\Gamma }`$. The map $`L_{(0.36,1.775)}`$ (abbreviated as $`L`$) has either one or three inverses depending on the phase-plane point in consideration (the term “$`Z_1Z_3`$ map” has been proposed for such maps, ). The invariant circle $`\mathrm{\Gamma }`$ has three first-rank preimages: itself, the curve $`\mathrm{\Gamma }_2^1`$ shown in the figure, and a third preimage $`\mathrm{\Gamma }_3^1`$, further away in phase space (not shown). The geometry of the map can be visualized by first folding the left side of the phase plane along $`J_0`$ and onto the right side. $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }_2^1`$ should now coincide. Next rotate roughly $`90^{}`$ so that $`J_0`$ maps to $`J_1`$. The image of $`\mathrm{\Gamma }`$ lands exactly on itself (before the folding), sending $`A`$ to $`L(A)`$ and $`B`$ to $`L(B)`$ in the process. Note that $`L`$ maps points that are in the intersection of the two regions bounded by $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }_2^1`$ to the region bounded by the pieces of $`\mathrm{\Gamma }`$ and $`J_1`$ which connect at $`L(A)`$ and $`L(B)`$. That is, some points inside $`\mathrm{\Gamma }`$ are mapped outside $`\mathrm{\Gamma }`$. Similarly, points inside $`\mathrm{\Gamma }_2^1`$, but not $`\mathrm{\Gamma }`$ are mapped from outside $`\mathrm{\Gamma }`$ to inside $`\mathrm{\Gamma }`$. Contrast this with the property that interiors of invariant circles for (orientation preserving) invertible maps are invariant.

Note also that $`\mathrm{\Gamma }`$ is tangent to $`J_1`$ at $`L(A)`$ and $`L(B)`$. The “folding” along $`J_0`$, typically creates a quadratic tangency along $`J_1`$ at images of curves crossing $`J_0`$, unless the curve crossing $`J_0`$ does so at an angle tangent to the zero eigenvector at that point. This phenomenon has been extensively discussed in . Our explanation here is based on the understanding reached in that paper. As a parameter (in this case $`\tau `$) is varied, the map $`L_{(0.36,\tau )}`$ (again abbreviated as $`L`$), the entire attracting set $`\mathrm{\Lambda }`$, the curves $`J_0`$ and $`J_1`$ as well as the intersection points of $`\mathrm{\Lambda }`$ with $`J_0`$ (e.g. $`A`$) and with $`J_1`$ (e.g. $`L(A)`$) vary. In our case, between $`\tau `$ values for figs. 2a and b, the invariant curve develops a “fjord’ which pushes across $`J_0`$, creating two new intersection points of $`\mathrm{\Lambda }`$ with $`J_0`$, and two new tangencies with $`J_1`$ at the two image points. We call the lower intersection point of the fjord $`C`$ (see fig. 2b). The points $`C`$ and $`L(C)`$ vary with $`\tau `$ as well. When the tangent of the “curve” in $`\mathrm{\Lambda }`$ at the intersection point $`C`$ becomes coincident with the null vector of the Jacobian of the map at $`C`$, the image of $`\mathrm{\Lambda }`$ ($`\mathrm{\Lambda }`$ itself) acquires a cusp touching $`J_1`$ at $`L(C)`$ . As the parameter value $`\tau `$ is further varied, the cusp becomes again a quadratic tangency, but $`\mathrm{\Lambda }`$ acquires a loop (figs. 2b and c). Iteration of this “loop” gives rise to an infinity of such loops (called “antennae” in ) on $`\mathrm{\Lambda }`$, suggesting chaotic behavior.

Further, Lorenz gives a heuristic justification of a simple criterion for an attractor (not formally defined, but assumed by us to be a compact invariant set with some sort of recurrence) to be a chaotic attractor (an attractor exhibiting sensitivity to initial conditions): the existence of two distinct points on the attracator which map to the same image point (also on the attractor) . In this context, the tracing out of an object with a single orbit suggests recurrence, and existence of a loop implies the existence of two points which map to the same image point. Thus, if this criterion is correct, attractors with loops would necessarily be chaotic attractors.

### 3.2 Numerical investigation

We set out to examine the transition from smooth IC to chaotic attractor in finer detail in figures 3 to 11. Before we embark on a detailed description, we should state that the real picture is yet more complex, involving similar phenomena on an even finer scale. For example, the same phenomena which cause the destruction of the “large” IC could also destroy the “small” IC’s which we refer to below. Nevertheless, this sequence of figures includes the essential features of the transition: coexistence of a periodic orbit with a nearby IC, destruction of the IC through a global bifurcation, appearance of loops on an unstable manifold, and the reappearance of an attractor, this time chaotic with loops. It might be useful to look ahead to the two-parameter bifurcation diagram in fig. 12 to put the one-parameter cut described in the current section in context.

We start from the smooth IC at $`\tau =1.775`$ (fig. 2a). At $`\tau 1.776243`$ a period-$`37`$ saddle and node pair appears in a saddle-node (turning point) bifurcation. The bifurcation occurs off of the IC, and results in the coexistence of two attractors: the stable invariant circle, $`\mathrm{\Gamma }`$, and an attracting period-$`37`$ sink. Although the sinks have real eigenvalues at the saddle-node bifurcation, the eigenvalues very quickly become complex and they have been labelled in fig. 3b for $`\tau =1.7765`$ by $`F_i`$ for “focus.” The separatrix between these two coexisting attractors consists of the stable manifolds of the period-$`37`$ saddles. Global unstable manifolds have been systematically constructed as (37th) forward iterates of local unstable manifolds as for invertible maps. In this case, we have plotted only the two branches emanating from saddle $`S_5`$. The right branch of the global unstable manifold directly approaches $`F_5`$. This is clear in the blowup of fig. 3b. The left branch very quickly approaches the IC $`\mathrm{\Gamma }`$. Successive 37th iterates travel slowly in a clockwise direction around the IC, eventually showing a good approximation to the whole IC. The approach to the IC is somewhat clearer in fig. 3b, where the branch appears to “join” the IC just to the left of $`S_5`$. The unstable manifold then continues almost “on” the IC: past $`F_{28}`$, $`S_{28}`$, and $`F_{14}`$, evolving around the IC, reentering from the right, just below $`F_5`$, before completing the (approximation to the) IC just to the left of $`S_5`$.

Notice also that 18 of the period-$`37`$ saddle-node pairs are “inside” (those marked by numbers inside $`\mathrm{\Gamma }`$) and the rest “outside” the invariant circle; this is another feature of noninvertible maps ). We discuss this again in Sec. 5.

As $`\tau `$ increases further, the consequences of global bifurcations are observed; these bifurcations involve interactions between the stable and unstable manifolds of the saddle period-$`37`$ points. The first such global bifurcation is evident in the sequence of figures 3 through 5: the left branch of the unstable manifold of $`S_{19}`$ misses $`S_5`$ to the left in fig. 3 and approaches the IC (the curve coming into the inset from the right and below $`F_5`$ is practically coincident with this unstable manifold branch); it makes the “saddle connection” with $`S_5`$ in fig. 4; it misses $`S_5`$ to the right in fig. 5 and approaches $`F_5`$.

The global crossing actually “starts” with a heteroclinic tangency (between a branch of a stable manifold of a period-$`37`$ saddle and and branch of an unstable manifold of another saddle on the same orbit) which destroys the IC at $`\tau 1.776878`$. This tangency is also called a crisis , since as $`\tau `$ increases toward the tangency parameter value, the IC (approached by the unstable manifold) and its basin boundary (the stable manifold) approach each other, colliding at the tangency parameter value. The attracting set during the crossing appears — as for invertible maps — to be a complicated topological object which at least includes the closure of the unstable manifolds of the period-$`37`$ saddles. Numerically, most orbits appear to approach the period-$`37`$ sink (fig. 4). Thirty-seventh iterates of each of the $`37`$ branches of the unstable manifold plotted in fig. 4a evolve clockwise toward the “next” saddle, suggesting that we are at least close to a heteroclinic crossing of stable and unstable manifolds. Although the stable manifolds of the saddles cannot be realistically computed because of the exponential explosion (“arborescence”) of the number of period-$`37`$ preimages, the crossing of the stable manifold is evident in the inset (fig. 4b). The left branch of the unstable manifold of $`S_{19}`$ heads toward saddle $`S_5`$, and although it is not clear from the figure, it begins to alternate back and forth, with folds on the right approaching $`F_5`$, and folds on the left approaching $`S_{28}`$. The portion approaching $`S_{28}`$ folds again, with one side heading toward $`F_{28}`$ (visible) and the other side heading toward $`S_{14}`$ (only the start is visible). Further iterates (not shown) would result in this one branch of the unstable manifold of $`S_{19}`$ eventually evolving all the way around the attracting set suggested in fig. 4a, repeatedly visiting every saddle and every node. As with invertible maps, this allows for long chaotic transients, if one could appropriately choose initial conditions on the unstable manifold. Most orbits, however, eventually are attracted to the period-$`37`$ sink.

We note that a transverse crossing of stable and unstable manifolds in the noninvertible setting does not imply transverse crossings along the rest of the orbit. For this reason, homoclinic or heteroclinic crossings do not necessarily imply the existence of a “horseshoe” or the chaotic transients referred to above. However, Sander  shows that the failure of a transverse crossing to provide a horseshoe is a codimension-two phenomenon, so we assume that our transverse crossings in our one-parameter study avoid such a problematic point, and behave as in the invertible case. Thus we expect the parameter value at the tangency initiating this heteroclinic crossing to be an upper bound on $`\tau _{CD}`$, marking the first appearance of chaotic dynamics. The invariant set with the chaotic dynamics does not appear to be an attractor, however, so the chaos is numerically evident only as transient behavior.

After the crossing, which ends in a heteroclinic tangency at $`\tau 1.776881`$, both branches of the unstable manifolds of the period-$`37`$ saddles approach the period-$`37`$ sink. The only attractor that remains is the stable period-$`37`$ sink, shown along with the period-$`37`$ saddles in fig. 5 for $`\tau =1.1778`$. The attracting set, however, is a “new” invariant circle, with a period-$`37`$ locking, formed by connecting the saddles and nodes with the unstable manifolds of the saddles. This invariant circle is seen in the blowup of fig. 5b. It is smooth except at the sinks, where two different branches connnect and the smoothness depends on the eigenvalues . Here the IC is only $`C^0`$ since the sinks are foci, with complex eigenvalues. There are apparently no cusps on the unstable manifolds at this parameter value.

As $`\tau `$ is increased, a secondary Hopf bifurcation of the period-$`37`$ foci at $`\tau =1.779444`$ renders them unstable. This can be seen in fig. 6 for $`\tau =1.78`$, where the unstable foci, are surrounded by stable period-$`37`$ invariant circles. It is important (and barely visible in the figure) that the unstable manifold intersects the small invariant circles (and itself!) repeatedly. This is a reminder of the noninvertibility of the map, and our acknowledged abuse of the term “manifold.” See  for a discussion of this type of unstable manifold self-intersection.

A number of high-period lockings on the period-$`37`$ invariant circles appear for very narrow ranges of the parameter. For example, at $`\tau =1.7835`$, a period-$`259`$ saddle-node pair (formed after a period-7 locking of the period-$`37`$ invariant circles) is visible (fig. 7); assuming the period-$`7`$ locking is on the small ICs, these small ICs are now frequency locked circles, composed of the chain of the unstable manifolds of the period-$`259`$ saddles (not shown). The period-$`37`$ unstable manifolds asymptotically approach the small ICs. The right branch clearly intersects itself and the IC it approaches (fig. 7b). Most individual orbits eventually approach the attracting period-$`259`$ sinks.

A second heteroclinic manifold crossing involving the “other” side of the global unstable manifold of the period-$`37`$ saddle “starts” at $`\tau 1.78373`$. The bifurcation is apparent in figures 7 through 9: the right branch of the unstable manifold of $`S_5`$ swings around $`F_5`$, and misses $`S_{19}`$ to the left in fig. 7, approaching the small IC around $`F_5`$; it makes the “saddle connection” with $`S_{19}`$ in fig. 8; it misses $`S_{19}`$ to the right in fig. 9, with successive iterates following essentially “on” the attractor with loops.

As for the earlier manifold crossing, understanding the saddle connection at $`\tau =1.784`$ requires some explanation. At this parameter value (fig. 8), a period-$`8`$ resonance on the small invariant circles results in the period-$`296`$ points apparently being the only attractors. The left branch of the unstable manifold of $`S_5`$ very quickly approaches the small IC surrounding $`F_{28}`$ and continues counterclockwise around it; it is not a homoclinic loop, although it appears close to being one. The right branch heads for $`S_{19}`$ after looping clockwise about $`3/4`$ of the way around $`F_5`$ (including the two high curvature “corners”); the global crossing is apparent since part of the manifold iterates to the left of $`S_{19}`$ around the IC (we iterated only half way around the IC), and part of the manifold iterates to the right of $`S_{19}`$ toward the right edge of the inset; this part to the right of $`S_{19}`$ continues the “long transient” behavior. It is important to notice that although both branches of the unstable manifold of $`S_{19}`$ are being followed by parts of the right branch of the unstable manifold of $`S_5`$, neither branch of the unstable manifold of $`S_{19}`$ has been directly plotted in this figure.

The slope of the unstable manifold at the point of intersection with $`J_0`$ continues to change as $`\tau `$ is increased. At $`\tau =\tau _{cusp}1.78428`$ the slope becomes equal to the slope of the eigenvector corresponding to the zero eigenvalue, and cusp points appear on the “manifold” (fig. 9c). Beyond this “cusp parameter” value, the unstable manifold contains loops. Since $`\tau _{cusp}`$ occurs while the unstable manifold is still intersecting the stable manifold, the picture for $`\tau `$ just beyond $`\tau _{cusp}`$, but not beyond the heteroclinic tangency which marks the end of the stable/unstable manifold crossing, is more complicated than might be expected in other families. In our case, the newly developed loops are damped out as they eventually approach the small period-$`37`$ ICs (possibly with locked solutions on the small ICs). Of course, the loops could follow a long transient before being attracted to the ICs. We do not have a figure corresponding to such a parameter value.

With further increase of $`\tau `$, we pass the “other end” of the manifold crossing (beyond the final heteroclinic tangency) and an apparently chaotic attractor appears. Thus, this final heteroclinic tangency parameter value appears to be an upper bound on $`\tau _{chaos}`$, marking the first appearance of a chaotic attractor. At $`\tau =1.785`$ the resulting attractor is seen in fig. 9 to coexist with a stable period-$`407`$ (born at a period-$`11`$ locking on the period-$`37`$ small invariant circles): the left branch of the unstable manifold of $`S_5`$ approaches the nearby small IC with the period-$`407`$ locking; the right branch approaches the large amplitude “messy” attractor that self-intersects and exhibits the loops it has inherited from the global unstable manifold. Such an attractor has been called a “weakly chaotic ring” (). Note that Lorenz plotted this weakly chaotic ring attractor by iterating a single orbit starting near the fixed point (fig. 1c); we plotted the same attractor by following a branch of the unstable manifold of the period-$`37`$ saddle orbit (fig. 9). Tracking the periodic saddles and their manifolds also allowed us to see the coexisting small IC with the period-$`407`$ attractor which was not apparent in fig. 1c.

The small amplitude invariant circles now begin to shrink, disappearing at $`\tau 1.78626168`$ in a second secondary Hopf bifurcation; the period-$`37`$ foci from return from sources to sinks; these continue to coexist with the large attractor that has become very “messy” with many more points of self-intersection (fig. 10b for $`\tau =1.7866`$). The stable period-$`37`$ foci approach the period-$`37`$ saddles in anticipation of the following period-$`37`$ saddle-node bifurcation, which happens at $`\tau 1.7866262`$.

Finally, at $`\tau =1.788`$ the saddle-node bifurcation of the period-$`37`$ solutions has occurred, and the pair has disappeared. Note that (as with invertible maps) the pairings have switched. For example, $`S_5`$ disappears with $`F_{28}`$, while $`F_5`$ disappears with $`S_{19}`$. Only the large chaotic attractor is left, seen in fig. 11.

Recap. The transition from smooth IC attractor ($`\tau =1.75`$, figs. 1b, 2a) to noninvertible “messy” chaotic attractor ($`\tau =1.788`$, fig. 11) is not a result of a smooth change of a smooth, large amplitude invariant circle which develops cusps at a critical parameter value. Instead, the transition occurs via a sequence of bifurcations, notably involving the interaction of the IC with nearby periodic orbits and the unstable manifolds of the saddles. Specifically, the IC as an attractor is destroyed in a global manifold crossing that leaves an attracting periodic orbit as the only attractor. The global unstable manifolds of the saddle period-$`37`$ points are the invariant objects that first develop cusps (at $`\tau _{cusp}1.78428`$, fig. 9c) and later loops (beyond $`\tau _{cusp}`$). When the cusps and loops first appear, they “damp out” as they approach the attracting periodic orbit (or small IC). When the large amplitude attractor reappears at after the second global manifold crossing, it is approached by the branch of the unstable manifold which now has loops. These loops are inherited by the large amplitude attractor which is now apparently chaotic.

## 4 The larger picture.

While a rich dynamical picture occurs along our one-parameter family with respect to $`\tau `$, it is clear that it only constitutes a part of a larger, actually two-parameter scenario. This is consistent with existing two-parameter descriptions of Hopf bifurcations for invertible maps . The primary features in the invertible case are resonance horns — inside which the corresponding maps have a periodic orbit of a certain period — emanating from the primary Hopf bifurcation curve. Since the fixed point at a Hopf bifurcation has eigenvalues on the unit circle, the IC close to its birth lies far away in phase space from $`J_0`$, where the map has a zero eigenvalue. Thus, we expect the two-parameter bifurcation diagram for the noninvertible Hopf bifurcation to strongly resemble the analogous bifurcation diagram for the invertible case — at least near the Hopf bifurcation curve.

We have embedded our one-parameter family ($`a=0.36`$ in eq. (3)) in the two-parameter family $`L_{(a,\tau )}`$. Fig. 12a shows part of the “big picture” in the $`\tau `$-$`a`$ plane: the period-one Hopf bifurcation curve together with period-$`37`$, period-$`5`$, and period-$`6`$ resonance horns emanating from it. The bifurcation points we encountered along the one-parameter cut “continue” to bifurcation curves in the two-parameter setting. The curves displayed in this figure are only a tiny sampling of the dozens of bifurcation curves we have computed, which are, in turn, only a tiny sampling of the bifurcation curves that can be computed for this family. (In particular, there are saddle-node curves in addition to those emanating from the primary Hopf bifurcation curve. See  for descriptions of such intraresonance horn features for invertible maps.)

As we move higher up in the horns — above the “IC tangent to $`J_0`$” curve in fig. 12a — the circle grows and starts to interact with $`J_0`$ and thus becomes “truly” noninvertible. It is above this curve where the corresponding maps can exhibit noninvertible behavior, such as mapping some of the interior of an IC to its exterior and vice versa (recall fig. 2a). More “interesting” behavior, such as coexisting attractors or the breakup of the IC, is suggested when the resonance horns start to overlap. In our one-parameter cut this overlap happens at least by the first period-$`37`$ sadde-node bifurcation at $`\tau 1.776243`$, although the three horns we have displayed do not overlap until higher $`\tau `$ values.

The enlargement in fig. 12b shows three of the multitude of bifurcation curves which subdivide the period-$`37`$ horn into “inequivalent” regions: two “equal eigenvalue” curves (just inside the horn), marking the transition between real and complex eigenvalues for a period-$`37`$ orbit, and a secondary Hopf curve (HB<sub>37</sub>), marking the birth of the “period-$`37`$ small ICs.” Somewhat unsatisfying in fig. 12b is the fact that none of the displayed bifurcation curves is unique to noninvertible maps. We remark that we did compute an “eigenvalue zero” curve (for the existence of a period-$`37`$ point with an eigenvalue of zero), but we did not display it in fig. 12b because it had no direct bearing on the bifurcation sequence we described (but see such a curve in a period-$`6`$ horn in fig. 13). However, this curve certainly has an indirect bearing on the bifurcations, because it is related to the crossing of $`J_0`$ by the attracting set (which includes any ICs, periodic orbits, and their unstable manifolds). If the attracting set never interacted with $`J_0`$, the bifurcations could be produced with invertible maps.

Other bifurcation curves, which our current computational tools did not handle easily, would separate the other $`\times `$’s in fig. 12b from each other. For example, the $`\times `$ labelled 4 is in a thin heteroclinic region; heteroclinic tangency curves separate it from 3 below and 5 above. Similarly, the $`\times `$’s labelled 8 and 9c are in a thin heteroclinic region, separating them from 7 below and 9 above. A “cusp curve” (unique to noninvertible resonance horns) runs through 9c. The $`\times `$ labelled 7 is in a thin “period-$`259`$ secondary resonance horn,” which emanates from the secondary Hopf curve HB<sub>37</sub> and corresponds to a locking on the small period-$`37`$ ICs.

A more complete knowledge of the internal structure of the “noninvertible” resonance horns in the two-parameter setting is clearly the key to extending our current understanding of the transition from a smooth invariant circle to a chaotic attractor.

## 5 Discussion

### 5.1 Typical one-parameter cuts

A natural question which arises following the detailed numerical investigation of Sec. 3 is how “typical” are the features of that specific one-parameter family. We take the liberty in the next few paragraphs of using the intuition we have gained via our numerical explorations, including the partial two-parameter bifurcation diagram of the previous section, to provide at least a partial answer to this question. Even though somewhat speculative, we feel that these observations are useful for placing the one-parameter cut in a broader context.

A key observation accompanying our understanding of the IC breakup is that an IC with a cusp and an irrational rotation number would have the cusp iterate forward to a dense set of image points, which would also be cusps. Even if such an object is topologically possible, it appears to us to be highly unlikely. Because of this, we expect the invariant circle to break up before it reaches a parameter value where it develops a cusp. Thus, the IC must interact with some other object in the attracting set. In our experience, the other object in this “pre-cusp” breakup is a nearby periodic orbit and its accompanying stable and unstable manifolds. The related manifold crossings appear to be the “natural” way to first destroy an existing IC and to later reconstitute “itself” as a chaotic attractor.

This observation is consistent with the sequence of transitions we described in Sec. 3. Extending this one-parameter transition interval to the two-parameter setting, we imagine a transition band, roughly parallel to the ‘IC tangent to $`J_0`$’ curve, and passing between parameter values labelled ‘1b’ and ‘1d’ in fig.12a. Below the lower boundary of the transition band a smooth invariant circle exists (possibly locked, possibly coexisting with periodic orbit(s)). Above the upper boundary of the transition band, the corresponding maps are chaotic. Loosely, the transition band is a “thick version” of the “cusp curve” suggested by Lorenz.

Further supporting this transition band scenario is the prominent “opening” of the resonance horns in fig. 12a. This would be even more prominent if more resonance horns were plotted: most of the resonance horns appear to open up and overlap near the top of the upper boundary of the transition band. This is suggestive of, or at least consistent with, chaos.

The existence of the transition band also suggests why the period-$`37`$ orbit might have been the specific orbit that was involved in our one-parameter scenario: that was the horn which was passed through as the transition to chaos band was crossed. Vertical parameter cuts for values of $`a`$ other than $`a=0.36`$ could pass through a different periodic horn during this narrow transitional parameter interval. We note that Lorenz specifically chose $`a=0.36`$ because it appeared to “avoid periodic windows.” In the two-parameter context, this can be interpreted as avoiding low-period resonance horns. Thus, we might expect a “typical” one-parameter cut to often interact with a periodic orbit of a period less than $`37`$.

Alternative viable bifurcation sequences are suggested by the blowup of the period-$`37`$ horn in fig. 12b. For example, moving the cut to the left of $`a=0.3599`$ could avoid the secondary Hopf bifurcation curve HB<sub>37</sub>, and thus avoid the small invariant circles that appeared in the $`a=0.36`$ cut. Other differences in varying $`a`$ would be more obvious if we were able to compute the region(s) of the period-$`37`$ horn corresponding to the global bifurcations, and curves corresponding to the development of cusps. (We expect to be able to do this in the future.) We know the $`\times `$ labelled ‘$`9c`$’ in fig. 12b is both on a cusp curve (in the parameter space) and in a heteroclinic region. Moving the vertical cut either left or right would probably prevent these two from coinciding. That is, the fact that the cusp parameter value occured during the manifold crossing in the $`a=0.36`$ cut is probably not typical of most one-parameter cuts. If the cusp parameter value occured first (with respect to $`\tau `$), then the phase portrait for the first appearance of loops would be simpler: the loops would simply approach the small attracting ICs. (If we encountered the cusp parameter outside the secondary Hopf curve, it would be still simpler: the loops would damp out at the attracting period-$`37`$ nodes. See fig. 15 for comparison.) On the other hand, if the cusp parameter value occured after the heteroclinic region had been crossed, we claim (with no proof) that there must be another periodic orbit interacting with the IC as well. Perhaps its unstable manifold will be the one to develop a cusp and later undergo a global manifold crossing to reconstitute the chaotic attractor with loops.

### 5.2 Related noninvertible features

The relatively high period of the locking encountered along the one-parameter cut above makes the elucidation of the phase and parameter space structures cumbersome. To illustrate some of these phenomena more succinctly we perform computations in the neighborhood of a simpler (period-$`6`$) resonance horn.

Fig. 13 is a blow-up of this period-$`6`$ resonance horn, with two additional types of bifurcation curves plotted inside: (a) two global bifurcation curves (approximating a thin heteroclinic manifold crossing region associated with the stable and unstable manifolds of the period-$`6`$ saddle points) and (b) the curve on which the periodic points acquire an eigenvalue zero (i.e. one periodic point lies on $`J_0`$ and one on $`J_1`$) . This latter curve provides an indication of the interaction of the dynamics with noninvertibility. An additional curve (‘IC tangent to $`J_0`$’), extending beyond the resonance horn, is plotted. This is the locus, in parameter space, where the invariant circle starts interacting with noninvertibility, by first touching the critical curve $`J_0`$. Note that the ‘eigenvalue zero’ curve is tangent to the ‘IC tangent to $`J_0`$’ curve.

In fig. 14, we show a phase portrait for a parameter value inside the period-$`6`$ horn and above the ‘IC tangent to $`J_0`$’ curve ( $`(a,\tau )=(0.465,1.472)`$ ). The invariant circle crosses $`J_0`$ transversely at two points $`A`$ and $`B`$; it is easy to see that it will then become quadratically tangent to $`J_1`$ at the image points $`L(A)`$ and $`L(B)`$ (see ). Also illustrated in the figure is the relative location of the attracting invariant circle and one more of its three first-rank preimages. Clearly, one of the period-$`6`$ saddle-node pairs is located inside the attracting invariant circle while the remaining five are located outside. (Note that iterates of points to the left of $`J_0`$ switch sides of the IC; thus, there must be an even number of periodic points to the left of $`J_0`$.) We have also illustrated the portion of the interior of the invariant circle that is mapped “outside” after one iteration of the map, as well as the portion of phase space outside the invariant circle that gets mapped inside.

Our next observation is for another parameter value inside the period-$`6`$ horn: $`(a,\tau )=(0.4838,1.409)`$. In fig. 15 we have an example of an unstable manifold of a period-$`6`$ saddle with a cusp which damps out at the period-$`6`$ attracting orbit. Recall that in our one-parameter cut, the cusp on the period-$`37`$ unstable manifold developed in the middle of a heteroclinic crossing. Iterates of the cusp could therefore follow a long transient before eventually being damped out as it approached a period-$`37`$ sink. Fig. 15 is a simpler version of the damping out of cusps on an unstable manifold in the approach to a sink.

Finally, we have made an attempt to illustrate the complexity of the basins of attraction of the coexisting attractors: the stable period-$`6`$ points and the stable invariant circle. While for the purposes of this paper we will not dwell on the details of this complicated structure (for a discussion of the basic building blocks of such a picture see again  as well as ), we emphasize a few key elements:

* the basins of attraction are not simply connected;
* the stable “manifolds” consist of a large number (possibly, even an infinity) of branches (resulting from the multiple backward in time trajectories). Only the computed “immediate” local manifold of the period-$`6`$ saddles is shown in fig. 16; the images and preimages of them form the boundaries of the grey disconnected “islands” that belong to the basin of attraction of the period-$`6`$ foci.
* notice that the stable set intersects $`J_0`$. The images and preimages of the intersection points have been called “points of alternance”  on parts of the stable manifolds: these are points (on $`J_0`$ and its preimages) where two preimages of a branch of the stable set come together; the “direction of relative movement” of nearby points switches at such points of alternance.

### 5.3 Computational issues

As a final comment on the current state of our work on noninvertible systems, we note that computational tools for the accurate approximation and visualization of objects such as stable and unstable manifolds are crucial in both understanding what is possible and in analyzing specific examples. The possibility of multiple inverses can cause an exponential explosion in the number of distinct segments of an invariant curve or of stable manifolds of periodic points. The same is true for the critical set: all the images and preimages of the critical curve $`J_0`$. In three dimensions the computational issues are even further complicated: the $`J_0`$ curve becomes a surface, and its iterates and preiterates must be constructed. Standard one-dimensional continuation tools have to be extended (to simplicial continuation, or through software like PISCES ) to effectively compute, triangulate and visualize these surfaces and their intersections. Constructing scientific computing, database management and visualization tools for these objects in two and three dimensions is an important task, and to our knowledge, an open area of research. Development of these tools would facilitate the analysis of many physical systems that are inherently noninvertible: discrete-time adaptation, control and decision problems, ranging from Model Reference Adaptive Control systems to signal processing applications -in addition to more traditional noninvertible models arising in population dynamics and economics. These tools will also be applicable in the case of multiple forward trajectory systems (such as those arising in implicit integrators, or in model predictive control), as well as, more generally, in relations .

## 6 Conclusions

The breakup of invariant circles in invertible maps of the plane possess a relatively well-studied dynamic phenomenology; homoclinic tangles of the invariant manifolds of nearby periodic saddle points, which arise from lockings on or off the invariant circle, play an important role in the two-parameter picture. Invariant circles for noninvertible maps of the plane possess, in some sense, an even “richer” phenomenology. Initially, close to the Hopf bifurcation giving rise to them, they are locally invertible, and the Arnold horn scenario holds. Large amplitude noninvertible invariant circles, however, may interact with the critical curves of the map. This can (and does) lead to new transitions and bifurcations. We demonstrated that the one-parameter transition discovered in  actually constitutes a part of the interior structure of the Arnold horns in the noninvertible two-parameter setting; we placed it in the appropriate context by pointing out the role of nearby periodic points and their stable and unstable manifolds. The picture we have presented here is far from complete, however we hope that our work is a step toward a more complete understanding of the transition from smooth invariant circle to chaotic attractor in the presence of noninvertibility.

Acknowledgments. Besides being motivated by the original work of Professor E. Lorenz, we acknowledge extensive discussions with him when this paper was first being developed in the early and mid 1990’s. He declined, however, being a co-author of the paper, something that we believe he deserved, and we note this here. We also acknowledge discussions with Prof. R. Rico-Martinez of the Instituto Tecnologico de Celaya (Mexico) and Prof. C. Mira of the Complexe Scientifique de Rangueil, Toulouse, France. We also acknowledge the support of the National Science Foundation (IGK and BBP: Grant No. DMS-9973926), and the Swiss Office of Energy (CEF). Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the funding agencies.