text:
# A Disk Galaxy of Old Stars at z ∼ 2.51footnote 11footnote 1Based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation.

## 1. Introduction

Our understanding of the mechanisms and timescales of the formation of the very first luminous galaxies in the universe is still very uncertain. Many detailed studies of local luminous early-type galaxies favor formation at very early epochs and within a very short time span (e.g., Peebles 2002 and references therein). In particular, the $`\alpha `$ elements are enhanced in luminous ellipticals relative to Fe (e.g., Worthy, Faber, & González 1992), this enhancement is strongly correlated with the central velocity dispersion $`\sigma `$ (Trager et al., 2000), and the relation is similar for field and cluster ellipticals (Bernardi et al., 1998; Colless et al., 1999; Thomas, Maraston, & Bender, 2002). Galaxies with the oldest stellar populations are both the most luminous and the most $`\alpha `$ enhanced (Thomas, Maraston, & Bender, 2002). If these claims are valid, the most luminous galaxies formed within a quite short period very early in the history of the universe, less luminous early-type galaxies formed more gradually and somewhat later, and differences in star formation histories between ellipticals in the field and those in clusters are relatively minor. This last statement is also supported by fundamental-plane analyses (e.g., van Dokkum & Ellis 2003; Treu et al. 2001; but see also Treu et al. 2002). Such conclusions appear to be difficult to reconcile with standard interpretations of the results of the semi-analytic cold-dark-matter (CDM) models (e.g., Thomas & Kauffmann 1999; Kauffmann & Haehnelt 2000), which predict that, for the currently favored low-density cosmology ($`H_070`$, $`\mathrm{\Omega }_m0.3`$, $`\mathrm{\Omega }_{vac}0.7`$; assumed throughout this paper), most of the present-day bright ellipticals will have undergone final assembly between $`z=1`$ and $`z=2`$. However, there are at least two caveats to this statement: (1) One of the least certain parts of the semi-analytic CDM models is the treatment of star formation and resulting feedback effects, and it is possible that different assumptions in this area could mitigate at least some of the discrepancy; and (2) assembly at relatively low redshifts is not necessarily inconsistent with formation of most of the stars at much higher redshifts. Nevertheless, both the color-luminosity relationship for spheroids (Bower, Lucey, & Ellis, 1992; Ellis et al., 1997; van Dokkum et al., 2000) and the correlation of $`\alpha `$-element enhancement with luminosity show that the formation of the stellar population and the assembly of the galaxy cannot be completely independent processes. As Peacock (1999) has succinctly put it, “It seems as if the stars in ellipticals were formed at a time when the depth of the potential well that they would eventually inhabit was already determined.”

In this context, it is of considerable interest to seek out and to investigate the galaxies with the oldest stellar populations at any given redshift. While such galaxies may not be typical, they can establish one end of the range of realistic formation scenarios, and they likely have less complicated star-formation histories than most. Furthermore, they may, in fact, be more representative of the formation of the bulk of the stars in early-type galaxies generally than might at first seem likely, since a miniscule contamination of young stars in a galaxy can drastically alter its rest-frame ultraviolet colors and eliminate such galaxies from color-selected samples.

There has been much recent activity in observations and discussion of so-called extremely red objects (EROs; e.g., Liu et al. 2000; Daddi et al. 2000, 2003; Stockton 2001; Cimatti et al. 2003; Yan & Thompson 2003; Miyazaki et al. 2003; Gilbank et al. 2003; earlier references can be found in these papers). These objects, which are often defined to have $`RK6`$ or $`IK>4`$, include two disparate classes. Some are highly reddened high-redshift starbursts (e.g., Dey et al. 1999). The remainder are almost exclusively galaxies with $`z1.2`$ whose stellar populations are already ancient and have little admixture from more recent star formation (Stockton, Kellogg, & Ridgway, 1995; Dunlop et al., 1996; Spinrad et al., 1997; Liu et al., 2000; Soifer et al., 1999; Stockton, 2001). For convenience, we shall henceforth simply refer to this latter class as “old galaxies” (OGs); but note once again that this term refers to the age of the stellar population, with no necessary implication regarding the age of assembly of the galaxy itself, aside from the sorts of constraints mentioned above. In any case, these OGs seem to represent the earliest major episodes of star formation in the universe.

The actual ages of the dominant stellar population in the small number of identified OGs that have been studied in detail have generated some controversy because of disagreements about spectral synthesis models and their interpretation (see, e.g., Dunlop 1999 and references therein). However, a minimum age of 2 Gyr for the oldest objects at $`z1.5`$ is accepted by virtually everyone, and there are good reasons for believing that some, at least, have ages between 3 and 4 Gyr (Dunlop, 1999; Stockton, 2001; Nolan, Dunlop, & Jimenez, 2001; Nolan et al., 2002). If these latter age estimates are correct, then the precursors of these galaxies will still have essentially fully formed stellar populations at higher redshifts. For example, for our assumed cosmological parameters, galaxies at $`z=2.5`$ will be $`1.6`$ Gyr younger than at $`z=1.5`$. This means that $`z=1.5`$ galaxies for which the last major episode of star formation occurred between 3 and 4 Gyr earlier might be expected to be fairly well settled and have small average internal extinction at $`z=2.5`$, and probably at even higher redshifts.

Many surveys of EROs have made the assumption, explicit or not, that a dominant spheroidal component (in mass, if not in light) is a necessary condition for the presence of an old stellar population at high redshifts (e.g., Moriando, Cimatti, & Daddi, 2000; Smith et al., 2002b). In this paper, we describe the first object among those found in our search for OGs in radio source fields at $`z2.5`$ that we have been able to investigate in some detail, and we show that it is a luminous disk galaxy, with little or no bulge component, comprising very old stars.

## 2. Selecting Galaxies with Old Stellar Populations in Radio-Source Fields at $`z2.5`$

We have recently concluded a search for OGs in the fields of quasars at $`z1.5`$ (Stockton, 2001, Stockton et al., in preparation). Radio source fields are generally likely to include regions of higher-than-average density in the early universe (e.g., Best et al. 2003), in which processes of galaxy formation and evolution are expected to have proceeded more rapidly than in the field. More importantly, looking for companions to radio sources at a specific redshift allows us to choose redshifts for which the photometric diagnostics from standard broadband filters give the cleanest separation between old galaxy populations and other possible contaminants, such as the dusty starbursts mentioned above. The relevant key redshifts for our present purposes are $`z1.5`$ (4000 Å break just shortward of the $`J`$ band) and $`z2.5`$ (break between $`J`$ and $`H`$ bands).

Assume that, for the brighter OGs from our sample at $`z1.5`$, essentially all of their stars formed within 500 Myr of the Big Bang (corresponding to a stellar population age of 3.7 Gyr at $`z=1.5`$). For passive evolution of the stellar population, one finds from solar metallicity models (see Bruzual & Charlot 1993, 2003; we have used the 2000 versions of the models for analysis here unless otherwise specified) that the precursors of such galaxies at $`z2.5`$ would have $`J22.4`$, $`H20.5`$, and $`K^{}19.5`$. This assumption is deliberately an extreme one: if our assumed age (corresponding to $`z_{formation}9`$) is too high, then these magnitudes are all a bit pessimistic. Our approach has been first to obtain imaging to $`K^{}20`$ ($`10\sigma `$) and $`J23`$ ($`5\sigma `$), looking for objects in the field with $`JK^{}3`$. If any likely candidates are present, we obtain deeper imaging at $`J`$ and $`K^{}`$, and include imaging in the $`H`$ band. Most of the infrared observations (including all of those relevant to this paper) have been obtained with the infrared camera CISCO (Motohara et al., 2002) on the 8.2 m Subaru telescope on Mauna Kea. The detector is a Rockwell $`1024\times 1024`$ Hg-Cd-Te Hawaii-2 array, with a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$105 and thus a field of just under 2′. Additional $`R`$-band imaging was obtained with ESI on the Keck II telescope.

Our sample comprises 87 radio-source fields from the Texas catalog (Douglas et al., 1996) selected with the help of the NASA/IPAC Extragalactic Database (NED) to have redshifts in the range $`2.3<z<2.7`$ and $`A_B<1`$, based on the Galactic extinction estimates of Schlegel, Finkbeiner, & Davis (1998). The coverage we have been able to achieve so far on these fields has been largely determined by telescope scheduling and weather. We currently have adequate observations for 10 fields, of which at least 3 have objects that have well characterized spectral-energy distributions (SEDs) that closely match those expected for old stellar populations at the redshifts of the radio sources, with little intrinsic reddening. General results for this survey will be given elsewhere; here we concentrate on one galaxy, found in the field of the radio galaxy 4C 23.56, for which we have both strong constraints on the SED and adaptive-optics (AO) imaging. A finding chart for the object is shown in Fig. 1. Although it was found independently in our survey, this galaxy was previously noted as a very red galaxy by Knopp & Chambers (1996) in their comprehensive study of the field of 4C 23.56. The galaxy is number 68 in their Table 4, and we will refer to it as 4C 23.56 KC68, or, more briefly, simply as KC68.

The photometry for 4C 23.56 KC68 is shown in Fig. 2. The magnitudes, prior to Galactic reddening correction, are $`R>26.7`$ ($`3\sigma `$), $`I>26.5`$ ($`3\sigma `$), $`J=23.38\pm 0.12`$, $`H=20.91\pm 0.06`$, $`K^{}=19.90\pm 0.03`$. We have used the photometric redshift code Hyperz222http://webast.ast.obs-mip.fr/hyperz/ (Bolzonella, Miralles, & Pelló, 2000) to explore possible fits to the photometric data, and some of these model SEDs are also shown in Fig. 2. The fitting procedure allowed redshifts $`0z6`$, internal extinction corrections $`0A_V6`$ (assuming a Calzetti et al. 2000 reddening law), and a full range of solar-metallicity stellar population models, including both single-burst models as well as models incorporating various rates of continuing star formation. Within this parameter space, several general conclusions emerged:

1. Only single-burst models produced good fits. This is presumably because significant continuing star formation would weaken the sharp inflection shortward of the $`H`$-band point and/or violate the $`R`$ and $`I`$-band upper limits.
2. The redshift is strongly constrained to the range $`2.1z2.7`$. Figure 3 shows the $`\chi ^2`$ statistic for the fits as a function of redshift.
3. When redshift was left as a free parameter, the fits tended to home in on $`2.35`$, with a 99% confidence interval $`2.1z2.7`$. If one assumes the galaxy to be a member of a cluster associated with 4C 23.56 and therefore fixes the redshift at 2.483, a fairly respectable solution exists with an age of 2.6 Gyr and a modest internal extinction $`A_V=0.16`$ mag.
4. Although the lowest reduced $`\chi ^2`$ was obtained for an old (3.0 Gyr) population with $`A_V0.14`$, to some extent it is possible to trade younger populations for increased reddening. However, this exchange cannot be carried too far: for ages $`1`$ Gyr the rapid rise in the flux at rest-frame 2000 Å, coupled with the bluer colors longward of 4000 Å, prevents the underlying continuum from being massaged to give a reasonable fit to the photometry by any amount of Calzetti-like extinction.

Table 1 gives parameters for the three models shown in Fig. 2, which are representative of the range of plausible fits to the photometry.

Note that, because of the well-known age-metallicity degeneracy, we can obtain younger ages by allowing super-solar metallicities. Such a trade-off means, for example, that even though some of the ages in Table 1 exceed the age of the universe for the given redshift and our assumed cosmological parameters (i.e., 2.75 Gyr at $`z=2.34`$; 2.59 Gyr at $`z=2.48`$), these excess ages can be considered proxies for high metallicities. At a redshift of 2.48, an $`L^{}`$ elliptical with a 2.1-Gyr-old solar-metallicity stellar population would have $`K^{}20.5`$ in the observed frame, so the observed $`K^{}=19.8`$ (after correction for Galactic extinction) for KC68 implies a luminosity of about $`2L^{}`$. This mass associated with this luminosity is not strongly dependent on metallicity, the differential effect between solar and twice solar metallicities being a reduction in the inferred mass by about 10%.

## 3. Adaptive-Optics Imaging of 4C 23.56 KC68

The galaxy 4C 23.56 KC68 happens to lie $`20\mathrm{}`$ from a 14th mag star, and we have been able to obtain adaptive-optics (AO) imaging, using the Subaru 36-element, curvature-sensing AO system and the Infrared Camera and Spectrograph (IRCS; Kobayashi et al. 2000). The $`1024\times 1024`$ Aladdin3 array had a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$0225, giving a field of 23$`\stackrel{}{\mathrm{.}}`$0. The total exposure was 12,600 s, all in the $`K^{}`$ filter. Figure 4 shows a section of the AO-imaged field. The image cores have FWHM = 0$`\stackrel{}{\mathrm{.}}`$10, and there is no evidence for significant anisoplanicity effects across the relevant part of the field. Insets in Fig. 4 show larger-scale images of KC68 at two contrast levels. The galaxy clearly has a high degree of symmetry, showing that $`2L^{}`$ galaxies, comprising relaxed old stellar populations, can be in place by the time the universe is only 2.6 Gyr old.

We can explore this morphology in more detail by making model fits to the image. We have used C. Y. Peng’s GALFIT333http://zwicky.as.arizona.edu/$``$cyp/work/galfit/galfit.html/ routine (Peng et al., 2002) to fit two-dimensional models corresponding to de Vaucouleurs ($`r^{1/4}`$-law), exponential, and general Sérsic profiles(Sérsic, 1968; Caon, Capaccioli, & D’Onofrio, 1993), all convolved with the PSF. Figure 5 shows these models and the residuals after subtraction from the observed image. Although a pure exponential model is an excellent fit, the formal best fit is for a Sérsic model with an index $`n=1.49`$, an effective radius $`r_e=0\stackrel{}{\mathrm{.}}22`$ ($`=1.8`$ kpc), and an axial ratio $`b/a=0.33`$. These numbers tend to confirm one’s visual impression that the galaxy is disk-like in morphology. While one cannot exclude the possibility of a small $`r^{1/4}`$-law bulge component, all of our attempts to model the galaxy as a two-component ($`r^{1/4}`$-law + exponential) system resulted in the $`r^{1/4}`$-law-component being given an extremely large effective radius (regardless of the initial parameter settings), essentially smearing it into the sky background.

For a more quantitative view, Fig. 6 shows the surface-brightness profile of KC68, along with seeing-convolved model fits. It is clear from the left panel that an $`r^{1/4}`$-law model cannot adequately fit the profile, and that an exponential model gives a much better fit. The right panel shows the best Sérsic-profile fit.

## 4. Summary and Discussion

The identification of this galaxy and the other two (mentioned in §2) at $`z2.5`$ that we have found so far with strikingly dominant old stellar populations demonstrates that at least some massive galaxies have completed virtually all of their star formation at sufficiently early times that they are no longer heavily dust-enshrouded at this redshift. Of these, at least 4C 23.56 KC68 has also achieved a relaxed regular morphology. The surprise is that this relaxed configuration seems to be very close to a pure disk, rather than to a spheroid, or even to a spheroid + disk.

There is abundant evidence that the large ($`40`$%) fraction of S0 galaxies seen in rich clusters at the present epoch drops dramatically at higher redshifts (Dressler et al., 1997; Fasano et al., 2000), with the slack being taken up by spirals, and with the proportion of elliptical galaxies remaining roughly constant. The mechanisms for the apparent conversion of spirals into S0 galaxies have been widely discussed ever since Butcher & Oemler (1978) first noticed the increasing fraction of blue galaxies in clusters at high redshifts (see, e.g., Bicker, Fritze-v. Alvensleben, & Fricke 2002 for a recent example).

However, there is also evidence to support the classical view that at least a significant fraction of the S0 population consists of galaxies that formed coevally with the oldest ellipticals (Tinsley & Gunn, 1976; Ellis et al., 1997). Iye et al. (2000, 2003) have described an S0-like galaxy at $`z1.5`$. In addition, recent large surveys find that a substantial fraction of objects identified as EROs have disks (Smith et al. 2002a; Yan & Thompson 2003; Gilbank et al. 2003; see also Smith et al. 2002b), and the only galaxy in our own sample of galaxies at $`z1.5`$ with old stellar populations (Stockton, 2001) for which we have good AO imaging also is dominated by an exponential profile (Stockton et al., in preparation). At a higher redshifts, Labbé et al. (2003) have identified a sample of large, disk-like galaxies, one of which has fairly red colors (Franx et al., 2003) and is at $`z=2.94`$.

All of these cases appear still to have significant (and usually dominant) bulge components. Indeed, one of the principal points of Smith et al. (2002b) is that the galaxy they describe may be a spheroid in the process of acquiring a disk, validating a commonly proposed scenario for the evolution of early-type spirals. In other words, the reason that this particular galaxy is classified as an ERO lies in the dominance of the bulge population and the weakness and low star-formation rate of the disk component. Many other galaxies classified as EROs also show discrete star-forming regions, but in these cases, too, the bulge population seems to dominate the colors used for the classification.

The galaxy we are discussing here poses a much more radical challenge to our understanding of the formation of old stellar populations at high redshifts: KC68 appears to be a $`2L^{}`$ galaxy of old stars having a disk-like morphology, with no significant bulge component. As indicated below, this galaxy seems so difficult to fit into the standard cold-dark-matter (CDM) picture of early galaxy formation that it is worth briefly revisiting the steps by which we have reached this conclusion.

Could the redshift be seriously in error?—Certainly for low extinction, the photometric redshift seems quite unambiguous (see Fig. 3). For the Calzetti et al. (2000) reddening law, this conclusion remains for the highest extinctions we have considered. We cannot exclude the possibility that one could contrive a reddening law that could result in a reasonable fit at a different redshift, for sufficiently high extinction. However, the high degree of symmetry shown by the galaxy does not fit the fairly chaotic morphology usually shown by galaxies having large amounts of dust. We believe that a redshift outside of the range $`2.1z2.7`$ is extremely unlikely.

Are the stars really old?—We have seen in §2 that reddening alone cannot reproduce the observed SED, including the sharp inflection shortward of the $`H`$ band and the tight constraints on the flux at rest-frame 2000 Å and 2400 Å given by the $`R`$ and $`I`$-band upper limits. (This conclusion strictly assumes a Calzetti et al. 2000 reddening law; but the statement almost certainly holds for any plausible reddening law). It is necessary to have a strong break in the underlying spectrum, even if reddening is substantial. Younger ages produce progressively worse fits, and ages $`<1`$ Gyr do not seem at all likely, provided only that the stellar initial mass function (IMF) is not pathologically deficient in high-mass stars. Of course, one can make the stellar population as young as one wishes by an ad hoc truncation of the IMF above about $`1.8M_{}`$ (e.g., Broadhurst & Bouwens 2000). For standard IMFs (e.g., Salpeter or Miller-Scalo), increasing the metallicity above solar can potentially give a slightly younger age. Experiments with the recent Bruzual & Charlot (2003) models show that the gross SEDs of solar-metallicity models can be reproduced by models at 2.5 times solar metallicity with ages reduced by about 30%. However, any real stellar population will have a range of metallicities, and chemically consistent evolutionary synthesis models (e.g., Fritze-v. Alvensleben 1999; Jimenez et al. 2000) show that populations with mass-weighted metallicities that are much higher than solar can still have integrated fluxes dominated by stars with sub-solar metallicities in the rest-frame near-UV, where line blanketing is important. This effect means that ages estimated from solar-metallicity models are unlikely to be serious overestimates (Dunlop, 1999). In the case of KC68, there seems to be no way to avoid requiring almost all of the stars to be close to 1 Gyr old or older. Furthermore, such solutions are about a factor of 3 less probable than the best fit, which is for a very old (2–3 Gyr) population with little reddening. Thus, it is likely that the most recent significant star formation in 4C 23.56 KC68 occurred when the universe was $`<1`$ Gyr old, although we cannot completely exclude somewhat later times. (Note that if we were to assume the fairly extreme solution of a 1-Gyr-old population with internal reddening of $`A_V1.1`$ mag, the star formation would have to be essentially complete by $`z3.8`$, when the age of the universe would have been 1.6 Gyr). The important point remains, in any case, that there cannot have been significant recent star formation.

Is the morphology really disk-like?—As we have stated in §3, we cannot exclude the possibility of a small $`r^{1/4}`$-law component, but any such component seems to be insignificant in terms of the total luminosity of the system. The small projected axial ratio and the good fit to a pure exponential profile indicate that the system is strongly dominated by a disk.

In sum, the most straightforward and plausible interpretation remains that 4C 23.56 KC68 is a well-settled disk of old stars at a redshift close to that of 4C 23.56 itself. If this interpretation is indeed correct, what are the consequences? Firstly, something very close to the classical picture of monolithic collapse must have occurred, with the important proviso that vigorous star formation did not begin until after the gas had settled into a disk. Secondly, star formation must have proceeded at an extremely rapid rate in order to form the required $`3\times 10^{11}M_{}`$ of stars over a time span that could not possibly be much more than $`10^9`$ years and that was probably more like a few $`10^8`$ years. Thirdly, the star formation and processes associated with it must have been very efficient at processing and expelling the gas in order to have strongly quenched any subsequent significant star formation.

In hierarchical formation models, including the standard CDM picture, galaxy morphologies are largely a reflection of their merger and accretion histories up to the epoch at which they are observed (e.g., Steinmetz & Navarro 2003). Semi-analytic models of galaxy formation (e.g., Kauffmann & Charlot 1998) suggest that most massive spheroids underwent final assembly between redshifts of 1 and 2, and that most large spiral disks were gradually acquired by smooth accretion of cold gas at redshifts $`1`$ and lower, after the last major merger (Katz & Gunn, 1991; Navarro & White, 1994). There appears to be little place in this scenario for a galaxy like 4C 23.56 KC68, for at least two reasons: (1) the implied prodigious star formation rate, which has to have been sustained for a few $`10^8`$ years without the help of a major merging event (as indicated by the lack of a detectable bulge), and (2) the building up of a massive galaxy apparently by smooth accretion of cold gas, without any evidence for significant hierarchical merging.

We of course cannot say from this one example whether similar galaxies are common at high redshifts. If they are, they would in one sense be the perfect progenitors for spheroidals: they comprise natural reservoirs of stars formed over very short periods of time in deep potential wells at early epochs. In short, their stars have all of the properties needed to explain the color-luminosity and alpha-element-luminosity relations for spheroidals, provided only that the morphological transformation from disks to spheroidals results from mergers of a small number of massive disk galaxies (Kauffmann & Charlot, 1998). One obvious course for the future would be to explore the morphologies of galaxies selected in the same way as KC68 was with high resolution imaging, either with NICMOS on HST, or with laser-guide-star AO on large ground-based telescopes.

We thank Jenny Patience for assisting with some of the Subaru observations, and we thank the staff of the Subaru and Keck telescopes for their excellent support for our observing programs. We have benefited from discussions with Josh Barnes, John Kormendy, and Brent Tully. This research has been partially supported by NSF grant AST03-07335. It made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors recognize the very significant cultural role that the summit of Mauna Kea has within the indigenous Hawaiian community and are grateful to have had the opportunity to conduct observations from it.


format: text

text:
# 1 Introduction

## 1 Introduction

Peptide folding is currently attracting considerable attention. Recent advances in this area include the de novo design of two monomeric three-stranded antiparallel $`\beta `$-sheet peptides, Betanova and Beta3s. Peptides that have the ability to fold on their own and are well characterized experimentally are valuable not least as a testbed for theoretical models for protein folding. $`\beta `$-sheet peptides are particularly interesting in this respect, as $`\beta `$-sheet formation is more challenging to model than $`\alpha `$-helix formation. Therefore, it is no surprise that both Betanova and Beta3s have become the subject of computational studies. Simulations of peptide sequences that are somewhat similar to these and occur in natural proteins, so-called WW domains, have been reported, too. For a recent review of computational studies of peptide folding, see Granakaran et al.

Here we present a study of the C-terminal $`\beta `$-hairpin from the protein G B1 domain and a triple mutant of Betanova called LLM. The original Betanova, which is less stable than the peptide LLM, is considered too. These different sequences are studied using an all-atom model with a simplified interaction potential. An earlier version of this model was tested on the same $`\beta `$-hairpin and an $`\alpha `$-helix, the designed so-called F$`_\text{s}`$. The model was able to fold these two sequences and the folded population showed, in both cases, a temperature dependence comparable with experimental data. It should be pointed out that different sequences are studied using exactly the same parameters; the interaction potential is, like that of Kussell et al. but unlike many other simplified potentials for protein folding, entirely sequence-based. This is of importance even if only one sequence is studied, because it ensures that the formation and breaking of non-native bonds is not a neglected part of the dynamics.

## 2 Materials and Methods

The model we study is a revised version of an earlier model. It contains all atoms of the polypeptide chain, including hydrogens, but no explicit water molecules. All bond lengths, bond angles and peptide torsion angles ($`180^{}`$) are held fixed, so each amino acid has the Ramachandran torsion angles $`\varphi `$, $`\psi `$ and a number of side-chain torsion angles as its degrees of freedom (for Pro, $`\varphi `$ is held fixed at $`65^{}`$). All bond lengths and bond angles are the same as in the original model.

The potential function

$$E=E_{\text{ev}}+E_{\text{loc}}+E_{\text{hp}}+E_{\text{hb}}$$
(1)
is composed of four terms. The remaining part of this section describes these different terms, with emphasis on what is new compared with the earlier model. Energy parameters are quoted in dimensionless units. To set the energy scale of the model, we use the midpoint temperature for the $`\beta `$-hairpin as determined by Muñoz et al., $`T_\text{m}=297`$ K, which corresponds to $`kT0.440`$ in the model.

The first term in Eq. 1, $`E_{\text{ev}}`$, represents excluded-volume effects and has the form

$$E_{\text{ev}}=\kappa _{\text{ev}}\underset{i<j}{}\left[\frac{\lambda _{ij}(\sigma _i+\sigma _j)}{r_{ij}}\right]^{12},$$
(2)
where $`\kappa _{\text{ev}}=0.10`$ and $`\sigma _i=1.77`$, 1.75, 1.55, 1.42 and 1.00 Å for S, C, N, O and H atoms, respectively. The role of the parameter $`\lambda _{ij}`$ is to reduce the repulsion between non-local pairs; $`\lambda _{ij}=1`$ for all pairs connected by three covalent bonds and $`\lambda _{ij}=0.75`$ otherwise. The reason for using $`\lambda _{ij}<1`$ for non-local pairs is partly computational efficiency, and partly the restricted flexibility of chains with only torsional degrees of freedom. To speed up the calculations, the sum in Eq. 2 is evaluated using a cutoff of $`r_{ij}^\text{c}=4.3\lambda _{ij}`$ Å.

The second interaction term, $`E_{\text{loc}}`$, is new compared with the earlier model. By introducing this term and modifying $`\sigma _i`$ for C and N, we slightly adjusted the shape of the Ramachandran $`\varphi `$, $`\psi `$ distribution. $`E_{\text{loc}}`$ is a local electrostatic energy given by

$$E_{\text{loc}}=\kappa _{\text{loc}}\underset{I}{}\rho _I\left(\frac{q_iq_j}{r_{ij}^{(I)}/\mathrm{\AA }}\right),$$
(3)
where the outer sum runs over all non-Pro amino acids along the chain, and the inner sum represents the interaction between the partial charges of the backbone NH and CO groups within one amino acid (the sum has four terms: NC, NO, HC and HO). The partial charges are $`q_i=\pm 0.20`$ for H and N and $`q_i=\pm 0.42`$ for C and O. We put $`\kappa _{\text{loc}}=125`$, which corresponds to a dielectric constant of $`ϵ_r2.0`$ if $`\rho _I=1`$. The factor $`\rho _I`$ reduces the interaction strength for the two end amino acids and Gly, which can be viewed as a simple form of context dependence; $`\rho _I=0.25`$ for end amino acids, $`\rho _I=0.5`$ for Gly, and $`\rho _I=1`$ otherwise. A similar factor is used for $`E_{\text{hb}}`$ (see below).

The third term in Eq. 1, $`E_{\text{hp}}`$, is an effective attraction between hydrophobic side chains that are not nearest or next-nearest neighbors along the chain. It has the pairwise additive form

$$E_{\text{hp}}=M_{IJ}C_{IJ},$$
(4)
where $`C_{IJ}`$ is a measure of the degree of contact between side chains $`I`$ and $`J`$, and $`M_{IJ}`$ sets the energy that a pair in contact gets. The contact measure $`C_{IJ}`$ is a number between 0 and 1, defined as before. The interaction matrix $`M_{IJ}`$ is given in Table I and differs from that used in our earlier study, which was based on the Miyazawa-Jernigan contact energies. With an all-atom representation, this cannot be expected to be a good choice for more general sequences, since the Miyazawa-Jernigan contact energies were derived using a different, reduced chain representation. The new matrix $`M_{IJ}`$ has a simplified structure in that the hydrophobic amino acids are grouped into three classes (see Table I). The $`M_{IJ}`$ values are taken to be large for the aromatic class (Phe, Trp, Tyr), which in part is an attempt to compensate for the fact that it is relatively difficult for these large side chains with few degrees of freedom to make proper contacts.

The last term of the potential, the hydrogen-bond energy $`E_{\text{hb}}`$, is given by

$$E_{\text{hb}}=ϵ_{\text{hb}}^{(1)}\underset{\mathrm{bb}\mathrm{bb}}{}\rho _{ij}u(r_{ij})v(\alpha _{ij},\beta _{ij})+ϵ_{\text{hb}}^{(2)}\underset{\mathrm{sc}\mathrm{bb}}{}\rho _{ij}u(r_{ij})v(\alpha _{ij},\beta _{ij}),$$
(5)
where the two terms represent backbone-backbone interactions and interactions between the backbone and charged side chains, respectively. The second term in Eq. 5 does not include any side chain-side chain interactions, as it did in our earlier study. Apart from that, the only difference compared with the earlier model is the factor $`\rho _{ij}`$, which like $`\rho _I`$ in Eq. 3 can be seen as a simple form of context dependence. We put $`\rho _{ij}=0.25`$ if any of the two amino acids involved is an end amino acid, $`\rho _{ij}=0.5`$ if any of them is a Gly, and $`\rho _{ij}=1`$ otherwise. The constants $`ϵ_{\text{hb}}^{(1)}=3.1`$ and $`ϵ_{\text{hb}}^{(2)}=2.0`$ as well as the functions $`u`$ and $`v`$ are exactly the same as before.

To study the thermodynamic behavior of this model, we use simulated tempering, in which the temperature is a dynamical variable. Details on our implementation of this method can be found elsewhere. For a review of simulated tempering and other generalized-ensemble techniques for protein folding, see Hansmann and Okamoto. Eight different temperatures are studied, ranging from 284 K to 371 K. For the backbone degrees of freedom, we use three different elementary moves: first, the pivot move in which a single torsion angle is turned; second, a semi-local method that works with up to eight adjacent torsion angles, which are turned in a coordinated way; and third, a symmetry-based update of three randomly chosen backbone torsion angles. For the side-chain degrees of freedom, we use simple Metropolis updates of individual angles.

For each peptide, eight independent Monte Carlo runs were performed, starting from random conformations. Each run required a few days on a standard desktop computer, and contained several folding/unfolding events. The similarity between the results from the different runs strongly suggest that the simulations did map out all relevant free-energy minima of the model. All statistical errors quoted are 1$`\sigma `$ errors obtained from the variance between the runs. The fits of data discussed in the next section were carried out by using a Levenberg-Marquardt procedure.

For a given protein structure, there generally exist alternative structures with similar secondary-structure content but different overall topologies. This holds true even for a small $`\beta `$-hairpin, for which a flip of the side chains gives rise to a topologically distinct structure. To make models discriminate between different topologies is a delicate task. To assess whether or not a model is able to do that, it is necessary to make a suitable choice of observables. In our calculations, we monitor two variables that can be used for this purpose: first, the root-mean-square deviation (rmsd) from the folded structure, $`\mathrm{\Delta }`$, calculated over all non-H atoms (a backbone rmsd is much less informative); and second, the number of native backbone-backbone hydrogen bonds, $`N_{\text{hb}}^{\text{nat}}`$. Figure 1 illustrates which hydrogen bonds we take to be present in the native states of the peptides studied.

In our calculations, a hydrogen bond is considered formed if the energy is less than $`ϵ_{\text{hb}}^{(1)}/3`$ (see Eq. 5).

Using the original model, we studied the $`\alpha `$-helical F$`_\text{s}`$ peptide and a $`\beta `$-hairpin. Here, we study the same $`\beta `$-hairpin and two three-stranded $`\beta `$-sheet peptides, LLM and Betanova. Before turning to these results, it should be pointed out that the F$`_\text{s}`$ sequence still makes an $`\alpha `$-helix in the revised model, as can be seen from the free energy $`F(\mathrm{\Delta },E)`$ in Fig. 2. $`F(\mathrm{\Delta },E)`$ has a pronounced, dominating minimum at $`\mathrm{\Delta }2`$–3.5 Å, which corresponds to $`\alpha `$-helix. In addition, there are weakly populated minima corresponding to $`\beta `$-sheet structures at $`\mathrm{\Delta }9`$ Å and $`\mathrm{\Delta }12`$ Å.

## 3 Results and Discussion

### 3.1 $`\beta `$-Hairpin

Using the model described in the previous section, we first study the 16-amino acid C-terminal $`\beta `$-hairpin from the protein G B1 domain. An important quantity monitored in our earlier study of this peptide was the hydrophobicity energy $`E_{\text{hp}}`$. This variable should be strongly correlated with Trp fluorescence, which Muñoz et al.  used to characterize the melting behavior of this peptide. The temperature dependence of $`E_{\text{hp}}`$ was found to be in reasonable agreement with the data of Muñoz et al. Several other groups have performed atomic simulations of the same $`\beta `$-hairpin, with or without explicit water. In contrast to ours, most models seem to require further calibration in order not to show a temperature dependence much weaker than that of experimental data.

Figure 3a shows the temperature dependence of $`E_{\text{hp}}`$ in the revised model.

The line is a fit of the data to a simple (first-order) two-state expression. The parameters of the fit are the midpoint temperature $`T_\text{m}`$, the energy difference $`\mathrm{\Delta }E`$, and two baselines. We use the parameter $`T_\text{m}`$ to set the energy scale of the model; this parameter is taken as $`T_\text{m}=297`$ K as determined by Muñoz et al. For the energy difference, we then obtain $`\mathrm{\Delta }E=13.1`$ kcal/mol. These values of the two-state parameters $`T_\text{m}`$ and $`\mathrm{\Delta }E`$ correspond to a native population of 74% at $`T=284`$ K, which agrees well with the result of Muñoz et al., 72% at $`T=284`$ K. The NMR analysis of Blanco et al. gave, by contrast, a lower native population, 42% at $`T=278`$ K. A possible explanation of this discrepancy would be that this peptide does not show a clear two-state behavior; the apparent native population may then very well depend on which quantity is studied. At first glance, this explanation may seem unlikely, given that the temperature dependence of the Trp fluorescence data to a good approximation showed two-state character. Let us therefore stress that, despite that the two-state fit in Fig. 3a looks quite good, this sequence does not show ideal two-state behavior in our model. This can be seen, for example, from the energy distribution, which lacks a clear bimodal shape. This was shown in our earlier study, and holds true in the revised model as well. Similar results have also been obtained in simulations of a designed, fast-folding three-helix-bundle protein.

In Fig. 3b we show the temperature dependence of the number of native hydrogen bonds, $`N_{\text{hb}}^{\text{nat}}`$, which we expect to be more strongly correlated than $`E_{\text{hp}}`$ with the NMR measurements of Blanco et al. For $`N_{\text{hb}}^{\text{nat}}`$, a two-state fit is not meaningful; for that, further data at lower temperatures would be needed. On the other hand, the quantity $`N_{\text{hb}}^{\text{nat}}`$ can be used as a direct measure of nativeness. Based on inspection of many examples, we use as a criterion for nativeness that at most two of the native hydrogen bonds should be missing, which can be used for the two other peptides too (see below). For the $`\beta `$-hairpin with seven native hydrogen bonds (see Fig. 1a), this criterion ($`N_{\text{hb}}^{\text{nat}}5`$) gives a native population of 39% at $`T=284`$ K. This value is close to the estimate of Blanco et al., 42% at $`T=278`$ K. Due to uncertainties about the precise definitions of nativeness and of when a hydrogen bond is formed, this agreement could be somewhat accidental. There is no doubt, however, that the native population obtained using $`N_{\text{hb}}^{\text{nat}}`$ is significantly lower than that obtained above using $`E_{\text{hp}}`$. Figure 4 shows the probability distributions of $`N_{\text{hb}}^{\text{nat}}`$ at $`T=284`$ K and $`T=306`$ K. The number of native hydrogen bonds is seen to rapidly decrease with increasing $`T`$, as it should.

The two-state parameter $`\mathrm{\Delta }E`$ extracted from the $`E_{\text{hp}}`$ data is somewhat smaller here, $`\mathrm{\Delta }E=13.1`$ kcal/mol, than it was in our earlier study, $`\mathrm{\Delta }E=16.1`$ kcal/mol. The reason for this is not so much that the model has changed, but rather that the fits were done in different ways. In our previous study, $`T_\text{m}`$ was held fixed at the specific heat maximum. Here, following the analysis of Muñoz et al. more closely, we take $`T_\text{m}`$ to be a parameter of the fit. The fitted value of $`T_\text{m}`$ turns out to lie slightly below (1–2%) the specific heat maximum. Our new analysis improves the agreement with the result of Muñoz et al., which was $`\mathrm{\Delta }E=11.6`$ kcal/mol.

Although the precise shape of the structures with lowest energy is sensitive to the details of the model, it is also interesting to make an rmsd-based comparison with experimental data. For this purpose, we use the NMR structure for the full protein G B1 domain (PDB code 1GB1, first model), as the NMR restraints for the isolated $`\beta `$-hairpin were insufficient to determine a unique structure. Figure 5a shows the free energy $`F(\mathrm{\Delta },E)`$ calculated as a function of rmsd, $`\mathrm{\Delta }`$, and energy, $`E`$, at $`T=284`$ K.

Three distinct, highly populated minima can be seen. The two minima with lowest $`E`$ are found at $`\mathrm{\Delta }2.0`$Å and $`\mathrm{\Delta }3.1`$Å, respectively. Both these correspond to $`\beta `$-hairpin structures with a high $`N_{\text{hb}}^{\text{nat}}`$. That $`N_{\text{hb}}^{\text{nat}}`$ is high implies, in particular, that the topology of the $`\beta `$-hairpin is the native one. The main difference between these two minima lies in the shape of the turn. The third minimum, at $`\mathrm{\Delta }4.0`$Å, is somewhat higher in $`E`$ than the first two. This minimum is also dominated by $`\beta `$-hairpin structures with the native topology and many hydrogen bonds, but the two strands tend to be out of register with each other, so $`N_{\text{hb}}^{\text{nat}}`$ is low. Largely, it is the existence of this third minimum that makes the apparent native population depend on which of the observables $`E_{\text{hp}}`$ and $`N_{\text{hb}}^{\text{nat}}`$ we use. Finally, there are also two weakly populated free-energy minima corresponding to $`\beta `$-sheet structures with the non-native topology ($`\mathrm{\Delta }5.3`$ Å) and $`\alpha `$-helix ($`\mathrm{\Delta }8`$–10 Å), respectively.

### 3.2 Three-Stranded $`\beta `$-Sheets

The de novo design of the 20-amino acid three-stranded antiparallel $`\beta `$-sheet peptide Betanova was reported in 1998. Recently, mutants of this peptide with higher stability were created by López de la Paz et al. Among the most stable mutants found was the triple mutant LLM (Val5Leu, Asn12Leu, Thr17Met). The peptide LLM and the original Betanova were estimated to have native populations of 36% and 9%, respectively, at $`T=283`$ K, based on NMR data. Melting curves have, as far as we know, not been reported for these peptides.

Our simulations of LLM show first of all that this sequence does make a three-stranded antiparallel $`\beta `$-sheet in this model. This can be seen from Fig. 5b, which shows the free energy $`F(\mathrm{\Delta },E)`$ at $`T=284`$ K. The free energy has a broad minimum at $`\mathrm{\Delta }3`$–5 Å, corresponding to $`\beta `$-sheet structures with the native topology and a high $`N_{\text{hb}}^{\text{nat}}`$. The shape of the $`\beta `$-sheet varies within the minimum. At $`\mathrm{\Delta }3.4`$ Å, where the free energy is lowest, the $`\beta `$-sheet has a bent shape, which enables the chain to make strong hydrophobic contacts. At $`\mathrm{\Delta }4.5`$ Å, the $`\beta `$-sheet tends to be much flatter, which is hydrophobically disfavored but makes it possible for the chain to form more perfect hydrogen bonds. There is also a free-energy minimum at $`\mathrm{\Delta }6.5`$ Å, which corresponds to three-stranded antiparallel $`\beta `$-sheet structures with the non-native topology. However, the native topology is the thermodynamically favored one. Note that the native and non-native topologies exhibit non-overlapping sets of backbone-backbone hydrogen bonds, so $`N_{\text{hb}}^{\text{nat}}`$ is low at the $`\mathrm{\Delta }6.5`$ Å minimum.

The main reason why the model favors the native topology over the non-native one lies in the side-chain orientations for the hydrophobic pairs Trp3-Leu12 and Leu5-Tyr10. The Cα-Cβ vectors of these pairs point inwards in the non-native topology, which makes it difficult to achieve proper contacts between the side chains. This is much easier to accomplish in the native topology, where the Cα-Cβ vectors point outwards. Interestingly, the situation is similar for the $`\beta `$-hairpin above. The $`\beta `$-hairpin also has two pairs of hydrophobic side chains that are ‘bow-legged’ in the native topology and ‘knock-kneed’ in the non-native one.

Next we estimate the native population for LLM. As we want to compare with the NMR-based results of López de la Paz et al., we consider $`N_{\text{hb}}^{\text{nat}}`$ rather than $`E_{\text{hp}}`$. Figure 6a shows the $`N_{\text{hb}}^{\text{nat}}`$ distribution at $`T=284`$ K. In addition to the native and non-native peaks at high and low $`N_{\text{hb}}^{\text{nat}}`$, respectively, this distribution exhibits a third peak at $`N_{\text{hb}}^{\text{nat}}=4`$. The typical conformation at this peak contains only the first of the two native $`\beta `$-turns (see Fig. 1b). The second $`\beta `$-turn is less stable, as will be discussed below. Using the criterion that at most two native hydrogen bonds should be missing ($`N_{\text{hb}}^{\text{nat}}6`$), we obtain a native population of 38% at $`T=284`$ K for LLM, which agrees well with the result of López de la Paz et al., 36% at $`T=283`$ K. We also performed simulations of the original Betanova, and Fig. 6a shows the result for this sequence too. From this figure it is evident that Betanova is less stable than LLM. The probability that $`N_{\text{hb}}^{\text{nat}}6`$ is 14% for Betanova at $`T=284`$ K, which means that this criterion gives a native population close to the NMR-based result of López de la Paz et al.  not only for LLM but also for Betanova.

That the model predicts LLM to be more stable than Betanova is not surprising because LLM has a more pronounced hydrophobic core. The agreement with experimental data is, nevertheless, remarkably good, especially since these calculations do not involve any adjustable parameter; the energy scale of the model is fixed using melting data for the $`\beta `$-hairpin and is then left unchanged.

Figure 6b shows the frequencies of occurrence for the different native hydrogen bonds (see Fig. 1b) for LLM and Betanova. For Betanova, there is a clear difference between the hydrogen bonds involved in the first $`\beta `$-turn (1–4) and those involved in the second $`\beta `$-turn (5–8). The latter four occur infrequently, showing that the second $`\beta `$-turn is quite unstable, which is in line with the conclusions of López de la Paz et al. For LLM, the difference in stability between the two $`\beta `$-turns is less pronounced. However, hydrogen bond 7, which connects Met17 to Tyr10 (see Fig. 1b), is quite unstable. The reason for this is that the side chain of Met17 can make better contacts with other hydrophobic side chains if the strand is slightly bent. This bend makes it difficult for hydrogen bond 7 to form.

Finally, in Fig. 7 we show the temperature dependence of $`E_{\text{hp}}`$ and $`N_{\text{hb}}^{\text{nat}}`$ for LLM. As in the $`\beta `$-hairpin case, we find that a simple two-state fit provides a good description of the data for $`E_{\text{hp}}`$. The fitted values of the parameters $`T_\text{m}`$ and $`\mathrm{\Delta }E`$ are $`T_\text{m}=303`$ K and $`\mathrm{\Delta }E=13.0`$ kcal/mol, which means that the native population obtained from this fit is significantly higher than that obtained from the $`N_{\text{hb}}^{\text{nat}}`$ distribution (see Fig. 6a). So, the model predicts that the apparent native population depends on which observable is used for this sequence, too. We are not aware of any existing experimental data that support, or refute, this conclusion for LLM.

## 4 Conclusion

Using a novel all-atom model with a simplified sequence-based potential, we have investigated the equilibrium behaviors of three $`\beta `$-sheet peptides. We determined native populations for these peptides in two ways, from the distribution of the number of native hydrogen bonds ($`N_{\text{hb}}^{\text{nat}}`$) and from the temperature dependence of the hydrophobicity energy ($`E_{\text{hp}}`$). These two estimates were compared with experimental results based on NMR and Trp fluorescence, respectively. This comparison is summarized in Table II. The agreement with experimental data is good, which in particular means that the model to a good approximation is able to reproduce the relative stabilities of these three peptides, as obtained from the NMR measurements. In line with the experimental results on the $`\beta `$-hairpin, we find that the apparent native population depends on whether we use $`N_{\text{hb}}^{\text{nat}}`$ or $`E_{\text{hp}}`$. This reflects the fact that the melting transition is not a clear two-state transition in our model (for any of these three sequences). It is also worth noting that, despite that the two-state picture is an oversimplification, the temperature dependence of $`E_{\text{hp}}`$ is quite well described by a simple two-state expression (see Figs. 3a and 7a). Computational studies of the $`\beta `$-hairpin have also been performed by many other groups, but the temperature dependence obtained was typically too weak, as has been pointed out by Zhou et al.  Our model shows a temperature dependence which is in good agreement with experimental data.

Our study of these three different peptides was carried out using one and the same set of parameters. In addition, we showed that the F$`_\text{s}`$ peptide makes an $`\alpha `$-helix for this choice of parameters. While these results are very encouraging, it is important to stress that we do not expect the model to be directly applicable to other sequences. However, by confronting the model with new sequences, we hope it will be possible to refine the potential, and thereby further extend its applicability. The present study was a first step in this direction, in which the model was improved by studying LLM and Betanova. To make the model able to fold these sequences, many changes were made, several of which were minor. The two perhaps most important changes were the replacement of the old hydrophobicity matrix ($`M_{IJ}`$), and the introduction of a simple form of context dependence for the hydrogen bonds. Whether it will be possible to carry on this process to a point where the model correctly reproduces the thermodynamics of small proteins remains to be seen. One thing that probably will be necessary in order to achieve this goal is to include multibody effects in the hydrophobicity potential; the present pairwise additive potential is likely to become insufficient as the chains get larger. Computationally, there is room for extending the calculations to larger chains; the calculations presented here required about two weeks on a standard desktop computer for each peptide.

Acknowledgments: We thank Luis Serrano and Manuela López de la Paz for providing NMR data for LLM and Betanova. This work was in part supported by the Swedish Foundation for Strategic Research and the Swedish Research Council.


text:
# Effects of Finite Deformed Length in Carbon Nanotubes

## Abstract

The effect of finite deformed length is demonstrated by squashing an armchair (10,10) single-walled carbon nanotube with two finite tips. Only when the deformed length is long enough, an effectual metal-semiconductor-metal heterojunction can be formed in the metallic tube. The effect of finite deformed length is explained by the quantum tunnelling effect. Furthermore, some conceptual designs of nanoscale devices are proposed from the metal-semiconductor-metal heterojunction.

The discovery of carbon nanotubesSI  offers a natural quasi one dimensional material to explore and verify quantum physics theory, and the advantage in design of nanodevice. Recently, an interesting topic is how a mechanical deformation of a single-walled carbon nanotube (SWNT) affects its electronic property.TWT ; EDM ; JC ; AM ; MBN  Some previous studiesAM ; MBN  show that in many cases, a local deformation can not drive a metal-to-semiconductor transition (MST) in armchair SWNTs. At the same time, our recent work JQL  indicates that a MST can be achieved by a radial deformation along the holistic armchair SWNTs. The different conclusions drawn from the two models show that the finiteness of the deformed length plays an important role in the transport properties of SWNTs. So, for designing electronic devices based on nanotubes, the deformed length along the tube will be a key parameter to ensure their functions. Here, in this letter, we will demonstrate the effects of the deformed length by squahing an armchair (10,10) SWNT.

In order to simplify the theoretical analysis, we consider an armchair (10,10) SWNT symmetrically squashed by two identical tips (with a width of $`d_x=4.0`$ Å and a length of $`d_z`$) along the $`\pm y`$ direction till $`d_y=2.2`$ Å,JQL  as shown in Figs. 1a and 1b. The atomic structures, or the shape of tubes, for different tip lengths $`d_z`$, are optimized by the tight-binding (TB) molecular dynamics (MD) code.CHX ; JQL  Fig. 1a shows an example optimized structure for tip length $`d_z=6.0`$ Å. Noticeably, the cross section of the tube changes from a dumbbell (Fig. 1b, corresponding to fully deformed part of the tube) to an ellipse (Fig. 1c, less deformed part), and then to a circle (Fig. 1d, undeformed part). The smallest interatomic distance between the upper face and the lower face is $`d_{AA^{}}=2.1`$ Å, and the mirror symmetry is broken due to small $`d_y`$. In this case, the squashed part should work as semiconductors.JQL

After obtaining the optimized structure with different tip lengths, we employ the TB Green’s function methodJCC ; JW ; JW2  to study the electronic transport properties of the deformed tubes. Fig. 2a presents the conductance curves near the Fermi energy $`E_F`$ of tubes deformed by tips with different lengths $`d_z`$. We can see that no energy gap is opened when $`d_z`$ is $`3.0`$ nm (the pink short dash line). In other words, the deformed tube keeps its metallic behavior. Moreover, the conductance of the deformed tube almost keeps constant near $`E_F`$. Only when the tip length increases, the conductance near $`E_F`$ will begin to decrease, as shown in Fig. 2a. It is easy to find that the conductance near $`E_F`$ approaches zero when $`d_z=25.2`$ nm (the red dash line). As the limit, when the length of tips approaches infinite, the whole tube is deformed, and then the conductance near $`E_F`$ will be zero.

The above results demonstrate that the deformed length in SWNTs plays an important role in their electronic properties. In order to reveal the effects of finite deformed length in SWNTs, we present a visual picture by the local density of states (LDOS)JW  per unit cell layer along the tube axes. Each unit cell layer consists of 40 atoms which distribute over two atomic layers. When $`d_z=25.2`$ nm, the LDOS of 120 unit cell layers (4800 atoms) is plotted in Fig. 2b, as the deformed part of the tube contains about 100 unit cell layers. It is obvious that the LDOS of the deformed part is almost zero near $`E_F`$. Energy gap exists near $`E_F`$ within the deformed part. While, the LDOS crosses $`E_F`$ continuously in the nearby part. Consequently, the squashed tube forms a metal-semiconductor-metal (MSM) heterojunction. Meanwhile, when $`d_z=3.0`$ nm, the LDOS of 30 unit cell layers (1200 atoms) is plotted in Fig. 2c, since the deformed part only includes about 12 unit cell layers. The LDOS crosses $`E_F`$ continuously in both the deformed part and its nearby. Hence, the tube keeps its metallic behavior.

Physically, the effect of finite deformed length can be understood by the quantum tunnelling effect. It is well known that the tunnelling probability attenuate exponentially with increasing width of energy barrier. Obviously, electrons can easily tunnel through the energy barrier if it is narrow enough. This also occurs for nanoelectronic devices. If the size of the device is too small, the electron will tunnel through the device, and consequently, the device can not function well as expected. The size effect is an essential problem for nanoelectronic devices design, as the tunnelling effect is universal.

To give a vivid sight of the tunnelling effect in the deformed tube, we present the LDOS of each atom within one unit cell layer of the dumbbell part in Figs. 3a and 3b. It is well known that the graphene sheet and hence the nanotube have two equivalent sublattices, which may be labelled as A and B sublattices. It has been pointed out that, when a MST is achieved in a squashed armchair SWNT, the LDOS distributing on the two sublattices are distinguishable, resulting in a discontinuity energy spectrum near $`E_F`$.JQL  This phenomena can also be found in the dumbbell part of tube deformed by tips with length $`d_z=25.2`$ nm, as shown in Fig. 3a. For the case of $`d_z=3.0`$ nm, the LDOS distribution of the dumbbell part is presented in Fig. 3b. We can find only a slight distinction between the LDOS of the two sublattices. While the LDOS crosses $`E_F`$ continuously. The LDOS near $`E_F`$ in the dumbbell part should result from the tunnelling effect.

More interesting thing is the change of electronic structure of the tube due to the finite deformation. The LDOS within a unit cell layer of the ellipse part (as shown in Fig. 1c) of the tube deformed by tips with $`d_z=25.2`$ nm is plotted in Fig. 3c. Though there is no energy gap near $`E_F`$ with the absence of the bond formation between the flattened faces, it is still interesting to find that the LDOS over the two sublattices are also distinguishable, which results from the distinction of the two sublattices in dumbbell part. To reveal this effect, in Fig. 3d, we plot the LDOS within a circle unit cell layer, which is far from the deformed center with distance 40.6 nm. To our surprise, again, the LDOS of the two sublattices are distinguishable. And more strange thing is the LDOS on one certain atom oscillates with energy, as indicated by the red dash line in Fig. 3d. Further studies show that the LDOS of the two sublattices within a unit cell layer are always distinguishable when the distance between the layer and the deformed center increases, and the oscillation of LDOS with energy becomes more exquisite. For an ideal SWNT, the LDOS near $`E_F`$ distribute over each carbon atom homogeneously. With the finite deformation, the $`\pi `$ and $`\pi ^{}`$ wave functions in the tube are mixed due to the breaking of the mirror symmetry.JQL  At the same time, as the translation symmetry along the tube is also broken, the standing wave functions are the combination of the incoming and reflecting wave functions,JW2  and the LDOS of them can be expressed as $`2[1\pm \mathrm{cos}(2\mathrm{\Delta }kx)][1+R\pm 2\sqrt{R}\mathrm{cos}(2k_F^{^{}}x+\delta )]`$.HFS  Obviously, the distinction between the LDOS of the two sublattices results from the difference between $`1\pm \mathrm{cos}(2\mathrm{\Delta }kx)`$. For a certain $`x`$, the oscillation of LDOS with energy results from $`\mathrm{cos}(2\mathrm{\Delta }kx)`$. As the energy dispersion relations near $`E_F`$ is approximately linear with $`\mathrm{\Delta }k=2E/(\sqrt{3}a\gamma _0)`$, the period of the oscillation is $`T_E=\sqrt{3}a\gamma _0\pi /(2x)`$.

From the above discussion, we know that an armchair (10,10) SWNT can form a MSM heterojunction when the deformation length is large enough ($`d_z20`$nm). For the conducting electrons, the middle semiconductor layers act as a energy barrier. Thus, based on the MSM heterojunction, we can easily structure some SWNT-based nanoscale devices.

Double barrier – To structure nanoscale double barrier, we only need to deform an armchair SWNT at two different positions, as shown schematically in Fig. 4a. As well known, the typical behavior of double barrier structure is the resonant tunnelling. We find that the conductance exhibits such a behavior (see Fig. 4a). The small conductance peaks lying in the energy gap are due to the resonant tunnelling effect. And, it is easy to control the electronic transport property of the nanoscale double barrier through changing the widths of the two energy barriers and the middle energy trap (i.e., changing the deformed length and the distance between the two deformed positions). For the example case shown in Fig. 4a , the tip length $`d_z=25.2`$ nm, and the distance between the two deformed positions $`l=57.1`$ nm.

Superlattice – If an armchair SWNT is deformed periodically along its axis, as shown schematically in Fig. 4b, it will behave as a superlattice. The conductance of the nanoscale superlattice is also presented in Fig. 4b. The period of the superlattice is $`l=57.1`$ nm, and the tip length is $`d_z=25.2`$ nm. Each conductance peak is split into many small peaks, as shown in the inset. This is the feature of superlattice structure, named as resonant splitting effect. It is also found that the conductance vanishes near $`E_F`$ between the conductance bands. We can control the electronic transport through changing the range of the conductance gap, which can be achieved by adjusting the structure parameters.

Transistor – For a logical circuit, transistor is one of the most important components. Designing a nanoscale transistor will be one of key steps to realize the nanologic. Recently, Bachtold et al. achieved a nanoscale transistor experimentally.AB  They adopted a semiconducting SWNT. The SWNT needs to be connected with two gold electrodes. Here, based on the MSM heterojunction, we can design pure carbon nanoscale transistor.

Fig. 4c is the schematics of a conceptual design of nanoscale transistor ($`d_z=25.2`$ nm). If a gate voltage $`V_g`$ is applied to the tips, the MSM heterojunction will act as a nanoscale transistor. When $`V_g=0`$, as the deformed part behaves as semiconductor, the conductance is zero and the transistor is OFF. With the change of $`V_g`$, the energy gap of the deformed part shifts. When $`E_F`$ does not lie in the energy gap any longer, the transistor begins to be ON. That is the basic principle of our nanoscale transistor. The dependence of the conductance with $`V_g`$ is also presented in Fig. 4c. We can find the logical voltage of the transistor is about $`0.10`$ V.

In summary, we demonstrate the effect of finite deformed length through squashing an armchair (10,10) SWNT. Only a deformation with enough long length can lead to an effectual MSM heterojunction in the metallic tube. The effect of finite deformed length is a general problem and is explained by the quantum tunneling effect. The change of the electronic structure, or the LDOS distributions, due to the finite deformation, shows the correlation between the nanodevices and leads. It should be considered in design of nanologic. Furthermore, based on the MSM heterojunction, some conceptual designs of nanoscale devices are proposed.

The work is supported by the National High Technology Research and Development Program of China (Grant No. 2002AA311153), the Ministry of Education of China, and the National Natural Science Foundation of China.


text:
# Robustness of non-abelian holonomic quantum gates against parametric noise

## I Introduction

The use of uniquely quantum phenomena to process information has led to surprising results in quantum key distributions crypto , information transfer protocols teleportation  and computation algos . From the point of view of the actual implementation of these theoretical protocols a main challenge is posed by the fact that generically quantum states are very delicate objects quite difficult to control with the required accuracy. The interaction with the many degree of freedom of the environment causes the loss of information (decoherence) and moreover errors in processing the information lead to a wrong output state (control errors). The first problem has been extensively studied over the past few years and a few ways to overcome it have been proposed and experimentally realized. These strategies include error avoiding error\_avoiding , error correcting strategies error\_correcting  and decoupling techniques bang-bang .

A new approach called topological quantum computation has been argued to be able to effectively solve both of them and open new ways to inherently robust quantum computation topological . Information is encoded in topological degrees of freedom of a system which are not sensitive to the local environment-noise effects and then are robust against decoherence tns . Unfortunately to the date, no simple feasible physical system has been identified to this aim; in fact, the systems proposed are usually complicated many-particle ones living over a macroscopic non-trivial structure (e.g. torus or cylinder topology). On the other hand, we can develop topological information processing, where the operator used depends on topological controls and then are robust against a the unwanted fluctuations of the driving fields. In this case, a important intermediate step is the geometrical quantum computation and particularly promising is the fully-geometrical approach called Holonomic Quantum Computation (HQC) HQC .

At variance with topological information processing, for geometrical QC several implementation proposals have been made; indeed the holonomic structure shows up in a variety of quantum systems, both in its Abelian (Berry) abelian  and non-Abelian version HQC\_proposal ; duan ; paper1 ; long\_pap . For this reason interest around HQC has recently growth leading also to proposals in which the adiabatic request can be relaxed and for implementing geometrical non-adiabatic quantum computation non\_adiabatic . One of the apparent advantages of HQC is that the gating time for holonomic gates does not depends on the logical operator applied but only on the adiabatic request; then HQC may lead to a new approach to implement complex operators difficult to construct with the standard dynamical gates (as discussed in long\_pap ). In view of its geometric nature i.e., dependence on areas spanned by loops in the control parameter manifold, HQC has been suggested to be robust against some class of errors preskill ; ellinas . Nevertheless thorough studies aimed to address this important issue are still relatively few and certainly not exhaustive HQC\_noise .

In this paper we will deal with the noise due to imprecise control of the system parameters during the evolution. This error source will be referred to as parametric noise. We will show that the fidelity of the holonomic gates displays three regimes, depending on the ratio between the adiabatic time and the noise correlation time. This results can be understood in view of the geometrical dependence of the holonomic operator. We will study in detail the class of holonomic gate proposed in paper1  in which the physical system are semiconductor quantum dots, the logical qubits are excitonic quantum state controlled by ultrafast lasers.

In section II, after a brief review of the holonomic approach, we describe the system used and the logical gate studied; moreover, it is discussed how we model the noise in the control parameters. In Sec. III we give a description of our simulations and show the results with different kind of noise processes for two single qubit gates and for a two qubit gate. A comparison with dynamical gates subject to the same noise is given too. Section IV contains the conclusions.

## II Holonomic quantum gates with parametric noise

Let us consider a family $``$ of isodegenerate Hamiltonians $`H(\lambda )`$ depending on $`m`$ dynamically controllable parameters $`\lambda `$, in the HQC approach HQC  one encodes the information in a $`n`$fold degenerate eigenspace $``$ of an Hamiltonian $`H(\lambda _0)`$. Varying the $`\lambda `$’s and driving adiabatically $`H(\lambda )`$ along a loop in the $`\lambda `$ manifold we produce a non-trivial transformation of the initial state $`|\psi _0U|\psi _0`$. These transformations, known as holonomies, are the generalization of Berry’s phase and can be computed in terms of the Wilczek-Zee gauge connection W-Z : $`U(C)=𝐏exp(_CA)`$ where $`C`$ is the loop in the parameter space and $`A=_{\mu =1}^mA_\mu d\lambda _\mu `$ is the $`u(n)`$valued connection. If $`|D_i(\lambda )`$ ($`i=1,\mathrm{},n`$) are the instantaneous eigenstates of $`H(\lambda )`$, the connection is $`(A_\mu )_{\alpha \beta }=D_\alpha |/\mathrm{\Omega }^\mu |D_\beta `$ ($`\alpha `$, $`\beta =1,\mathrm{},n`$). The set of holonomies associate with a given connection is known to be a subgroup of the group of all possible $`n`$-dimensional unitary transformations; when the dimension of this holonomy group and coincides with the dimension of $`U(n)`$ one is able to perform universal quantum computation with holonomies HQC .

For concreteness in this paper we will focus on the class of holonomic quantum gates analyzed in paper1 ; long\_pap . Logical qubits are given by polarized excitonic states controlled by femtosecond laser pulses. The parameters we have used in performing our simulations are those relevant to this specific kind of physical systems.

First we concentrate to one-qubit gates. Despite this might look at first as major limitation, we observe that, as far as the holonomic structure is concerned, the two-qubit gates are very similar. So we expect that most of the results we are going to present here e.g., the existence of separate regimes, should, to a large extent, hold true for two-qubit gates too.

The time-dependent interaction Hamiltonian in the interaction picture is

$`H_{int}`$ $`=`$ $`\mathrm{}(\mathrm{\Omega }_+|E_L^++\mathrm{\Omega }_{}|E_L^{}+\mathrm{\Omega }_0|E_L^0)G|`$ (1)
$`+`$ $`h.c.`$
where $`|E_L^i`$ ($`i=+,,0`$) are the polarized excitonic states (two logical and one ancilla) and $`|G`$ is the ground state (absence of exciton). This Hamiltonian family admits two dark states i.e., $`H_{int}(\mathrm{\Omega })|D_i(\mathrm{\Omega })=0,(i=1,2).`$ This two-fold degenerate manifold represent contains our encoded logical qubit: $`|0_L:=|E_L^+,|1_L:=|E_L^{}.`$ In Refs. duan ; paper1  it has been shown that the Wilczek-Zee connection associated to the Hamiltonian family (1) allows to construct universal one-qubit gates. These are realized by giving an explicit prescription for driving the control parameter $`\mathrm{\Omega }`$’s along suitable adiabatic loops.

We suppose now to add to the control $`\mathrm{\Omega }`$-field a ’small’ noise which perturbs our trajectory on the control manifold. We test the robustness of the geometrical mixing  single-qubit gate proposed in Ref. paper1 . To obtain this gate, we made the following loop in the parameter space : $`\mathrm{\Omega }_{}(t)=\mathrm{\Omega }\mathrm{sin}\theta \mathrm{cos}\phi `$, $`\mathrm{\Omega }_+(t)=\mathrm{\Omega }\mathrm{sin}\theta \mathrm{sin}\phi `$ and $`\mathrm{\Omega }_0(t)=\mathrm{\Omega }\mathrm{cos}\theta `$ (where the $`\mathrm{\Omega }`$ parameters is fixed and the $`\theta (t)`$ and $`\phi (t)`$ are time dependent) and the holonomic operator obtained at the end of the loop is $`U=e^{i\varphi \sigma _y}`$, (where $`i\sigma _y=|E_L^+E_L^{}||E_L^{}E_L^+|`$). The geometrical parameter $`\varphi =\mathrm{sin}\theta d\theta d\phi `$ is the solid angle spanned by the parameter vector $`\stackrel{}{\mathrm{\Omega }}=(\mathrm{\Omega }_+,\mathrm{\Omega }_{},\mathrm{\Omega }_0)`$ on the parameter manifold (sphere). Changing the relation between $`\theta `$ and $`\phi `$ we change the loop and then the value of $`\varphi `$. We choose a loop in order to obtain $`\varphi =\pi /2`$ and $`U=exp(i\pi \sigma _y/2)=|E_L^+E_L^{}||E_L^{}E_L^+|`$.

The logical operator $`U`$ depends only on geometrical parameters (i.e. solid angle swept on the parameters manifold), every perturbation that changes the trajectory in the control manifold changes the operator leading to computational errors. This notwithstanding, the perturbations leaving (almost) unchanged this angle will not affect the holonomic operator. Then we expect that even strong fluctuations - provided their time scale is sufficiently fast - average out, leaving in this way the angle unchanged and thus not affecting the computation.

The control parameters are the intensities and the phases of the lasers (since we suppose to be in resonant condition, we have fixed frequencies) perturbed by external noise (one for the phases and one for the intensities). To clarify which of these mainly affects the gate operation, first we separate the two types of errors (intensity and phase fluctuations) and then we apply both of them.

We use a straightforward model for the noise : we extract a random number from a probability distribution, we add a constant noisy field for the time $`T_n`$ to the evolving control field, then we extract another random number and so on. To simplify the simulations we choose $`T_n`$ in such a way that $`T_{ad}`$ is a multiple of $`T_n`$ ($`T_{ad}=n_rT_n`$). The fundamental parameter is the noise time $`T_n`$ that is the lapse of time of each random extraction; i.e. it represents the time scale of each random fluctuation.

For the ’intensity’ noise we modify only the value of the Rabi frequencies $`\mathrm{\Omega }`$. We have three lasers turned on and we suppose they have independent fluctuations $`\delta \mathrm{\Omega }_i`$ (i.e. every $`T_n`$ we extract three random numbers). The evolution on the control manifold is described by:

$$\{\begin{array}{ccc}\mathrm{\Omega }_{}(t)\hfill & =& \mathrm{\Omega }\mathrm{sin}\theta \mathrm{cos}\phi +\delta \mathrm{\Omega }_{}(t)\hfill \\ \mathrm{\Omega }_+(t)\hfill & =& \mathrm{\Omega }\mathrm{sin}\theta \mathrm{sin}\phi +\delta \mathrm{\Omega }_+(t)\hfill \\ \mathrm{\Omega }_0(t)\hfill & =& \mathrm{\Omega }\mathrm{cos}\theta +\delta \mathrm{\Omega }_0(t)\hfill \end{array}$$

In this case, the Rabi frequencies remain real parameters, but, if we introduce a ’phase’ noise ($`\mathrm{\Omega }_je^{i\xi _j}\mathrm{\Omega }_j`$), they acquire an imaginary part. The random numbers are taken as before and the evolution in the control manifold is:

$$\{\begin{array}{ccc}\mathrm{\Omega }_{}(t)\hfill & =& e^{i\xi _{}}\mathrm{\Omega }\mathrm{sin}\theta \mathrm{cos}\phi \hfill \\ \mathrm{\Omega }_+(t)\hfill & =& e^{i\xi _+}\mathrm{\Omega }\mathrm{sin}\theta \mathrm{sin}\phi \hfill \\ \mathrm{\Omega }_0(t)\hfill & =& e^{i\xi _0}\mathrm{\Omega }\mathrm{cos}\theta \hfill \end{array}$$

The most general and complicated situation is when both ’intensity’ and ’phase’ noise are present.

## III Simulations

For all the simulation we choose the parameters used in Refs. paper1 ; long\_pap  which satisfy the adiabatic condition : $`\mathrm{\Omega }=0.02`$ fs-1, $`T_{ad}=7.5`$ ps and $`\mathrm{\Omega }T_{ad}=150`$. The probability distribution for the noise is a Gaussian with zero mean and $`<\sigma >=<\frac{\delta \mathrm{\Omega }}{\mathrm{\Omega }}>=0.1`$ where $`\delta \mathrm{\Omega }`$ is the fluctuation of the Rabi frequency. Though this value is far below the experimental control, we use it to understand the robustness of holonomic gates against strong perturbations.

Once the evolution of the state with the noise has been computed, we must compare it with the ideal one in which the noise is absent. In order to make such a comparison quantitative we introduce the fidelity

$$=\sqrt{\psi _{out}^{id}|\rho _{out}|\psi _{out}^{id}}$$
(2)
where $`|\psi _{out}^{id}`$ is the final state without noise and $`\rho _{out}`$ is the density matrix of associated to the noisy final state. In our case the evolution is unitary, then $`\rho _{out}=|\psi _{out}^{noise}\psi _{out}^{noise}|`$ and the fidelity reduce to be a scalar product between the noisy and the ideal state. To eliminate the dependence of (2) on the initial states we make a sampling of the initial state space and then average the results. Even if we have a four dimensional working space, the initial state space has dimension two; in fact, in the ideal gates (once satisfied the adiabatic condition) we always start and end in a superposition of the logical states $`|E_L^+|E_L^{}`$. This simplifies the sampling procedure because we can take the Bloch sphere as sampling space. We sample the Bloch sphere with $`18`$ states sample  and for each of them calculate the fidelity. Fixed the random numbers extracted during the evolution (i.e. fixed $`T_n`$), we have also to take into account the dependence from the series of the random number extracted. For every sampled state we make five different realization with different extracted random numbers and calculate the average fidelity for the particular initial state. Finally, we make the average fidelity for all the sampled states to obtain the fidelity of our gate for a fixed $`T_n`$.

In the Figures presented here on the $`y`$ axis is plotted the fidelity and on the $`x`$ axis is plotted the number of extractions (that is the the ratio $`T_{ad}/T_n`$ between the adiabatic and noise time).

In Figures 1 we report the simulated evolution for two regimes of noise. In fig. 1 (A) we plot the fidelity when we extract up to $`100`$ random numbers during the evolution ($`T_{ad}/100T_nT_{ad}`$). Up to $`30`$ random numbers the average fidelity is $`0.875`$, while it decreases up to a minimum of about $`0.5`$ (with total average of $`0.632`$) if we extract more random numbers.

A possible interpretation of this effect can be given looking at Figure 2 where we show the evolution on the parameter sphere. In figure 2 (A) $`n_r=2`$ ($`T_n=T_{ad}/2`$) and we change the noise field twice during the evolution. Despite the intense noise, the shape of the loop is still clearly visible, it is simply shifted with respect to the ideal one. The value of the resulting solid angle swept is near to the ideal one. In figure 2 (B) we extract $`n_r=70`$ random extractions during the adiabatic evolution ($`T_n=T_{ad}/70`$). The fluctuations are too intense and too few to cancel out then, as we expected, the solid angles, swept, respectively during the ideal and the noisy loops, are different. This can explain the result in figure 1 (A).

As stated before, even if we have strong fluctuations, we expect that if we extract many random numbers the noise in average does not affect the solid angle and leave the holonomic operator unaffected. This is what seems to be confirmed by the simulations illustrated in Figure 1 (B) where we extract from $`50`$ to $`5000`$ random numbers. The fidelity increases and it is even better of those in figure 1 (A); the average fidelity is $`0.918`$ but it increases to $`0.956`$ if take into account the values from $`1000`$ to $`5000`$ random extractions. The relative loop in the parameter space is shown in figure (3) where we extract $`5000`$ random numbers during the evolution: the fluctuations have $`\delta \mathrm{\Omega }=0.1\mathrm{\Omega }`$ but are so quick that cancel out.

From these simulations it is evident that we have three regimes in our model which can be explained as the sign of the geometrical dependence of the holonomic operator:

* Slowly varying random fluctuations ($`T_{ad}/T_n1`$): the loop basically maintains its shape and it is simply shifted. This situation does not affect the gate too much.
* Intermediate regime ($`50T_{ad}/T_n100`$): the intense fluctuations badly modify modify the loop shape and alter the gate operator.
* Fast varying random fluctuations ($`T_{ad}/T_n1`$): the fluctuations effectively cancel out and do not change the operator.

These geometrical features are independent from the ratio $`\delta \mathrm{\Omega }/\mathrm{\Omega }`$ and persist even for small values of the fluctuation $`\delta \mathrm{\Omega }`$. If we decrease the intensity of the noise our adiabatic gate improves as shown in Figures (4) where, with the same parameters (i.e. Rabi frequency and adiabatic time), we choose $`\delta \mathrm{\Omega }=0.01\mathrm{\Omega }`$. We note that the fidelity increases but the features are the same of the previous simulations (i.e. three regimes).

In Figure (5) the noise is applied only to the phase of the control field with the same adiabatic parameters of the previous simulations. Since the way of producing the operator is always the same (i.e. loop in the parameter space with a noise), we expect that the effects of noise are the same of those for ’intensity’ noise. This can be clearly seen in figures (5) where we find the same features of the previous plots. The fidelity is much better respect to ’intensity’ noise and we can say the the ’phase’ noise does not affect our gate.

We finally discuss the case of both intensity and phase noise (Figure 6). As we expect, the main part of the error is given by the intensity noise . The three regimes, discussed previously, are evident also in this case.

In Figures (7) we show the population of the non-logical states ($`|G`$ and $`|E_L^0`$) at the end of the gate application as function of $`T_{ad}/T_n`$. In an ideal adiabatic gate (for $`T_{ad}\mathrm{}`$) these states are not populated. In Figure (7) (A) these populations increase with the $`T_{ad}/T_n`$ ratio due to the strong and fast fluctuations of the Rabi frequencies which perturb the Hamiltonian. For very fast varying fluctuations (Figure 7 (B)) the undesired populations decrease because of cancellation effects. From this analysis we can say that for slowly and fast varying random fields (i.e. the two extreme regimes) we remain in the logical computational space.

We wish to compare this holonomic gate with a standard dynamical gate, the latter one being characterized by the same unitary operator of the holonomic one and having different (typically shorter) gating time (the time need for the application of the gate). In our system this means to consider as logical states $`|G=|0`$ and $`|E_L^i=|1`$ and apply a $`\pi `$pulse laser sequence to produce a transition $`|G|E_L^i`$.

We made similar simulations for dynamical gate with the same parameter ($`\mathrm{\Omega }=0.02`$ fs-1) and with the same noise.

In Figures (8) we make a comparison between holonomic (solid line) and dynamical gates (dashed line) subject to the same intensity and phase noise. The comparison is not direct since the gating times are different. This means that the ratios $`T_{ad}/T_n`$ and $`T_{dyn}/T_n`$ are different for adiabatic and dynamical gates. To compare the effect of the gate subject to the same noisy field (i.e. with the same noise time $`T_n`$), we have to take into account that the dynamical gates, in our model, are about 100 times faster than the adiabatic ones (see Refs. paper1 ; long\_pap ). This means that if during the adiabatic evolution the noise changes $`n_r^{ad}=T_{ad}/T_n`$ times, for the dynamical gates it changes only $`n_r^{dyn}=T_{ad}/(100T_n)`$ times (if for the dynamical noise $`n_r^{dyn}<1`$, it changes only once). In other words, if for dynamical gate we extract $`n_r^{dyn}`$ random numbers during the evolution, for the holonomic gate we extract $`n_r^{ad}=100n_r^{dyn}`$ random numbers. This means that to compare the fidelities we have to look in the fast varying fluctuating noise ($`T_{ad}/T_n1`$) region of the previous figures.

For this reason in Figures 8 (A) and (B) we put two different scales for the adiabatic ($`T_{ad}/T_n`$) and dynamical ($`T_{dyn}/T_n`$) gates and the plots have been rescaled taking into account for the different gating time.

Numerical simulations in Fig. 8 show that the performance of holonomic gates and dynamical gates are comparable in the region where the first ones are reliable; that when $`T_{ad}/T_n1`$ (see discussion above).

Both dynamical and holonomic gates can be further improved. Since the dynamical gates are not subject to adiabatic constraints, we can choose different parameters in order to minimize the effect of the noise but this can affect the gating time.

For holonomic gates, given a noise with fixed correlation noise time $`T_n`$, we can try to change adiabatic time in order to modify the $`T_{ad}/T_n`$ ratio. This should allow us to fall in a ’good’ regime i.e. fast or slowly varying fluctuations. Decreasing adiabatic time to get in the small $`T_{ad}/T_n`$ region can produce new errors due to the lack of adiabaticity during the evolution and then it must be treat carefully. Increasing the adiabatic time to enter in the region of great $`T_{ad}/T_n`$ leads both to a better precision in the adiabatic gates and to the cancellation of noisy fluctuations but gives longer gating times.

To complete the set of universal quantum gates, we need another single qubit gate and a two qubit gate.

In the first case, we apply an intensity noise to the phase gate presented in Ref. paper1  with $`\mathrm{\Omega }_{}=0`$, $`\mathrm{\Omega }_+=\mathrm{\Omega }\mathrm{sin}(\theta /2)e^{i\phi }`$ and $`\mathrm{\Omega }_0=\mathrm{\Omega }\mathrm{cos}(\theta /2)`$. The holonomic operator is $`U=e^{i\varphi |E^+E^+|}`$ where $`\varphi =\frac{1}{2}\mathrm{sin}\theta d\theta d\phi `$. The results are shown in Figures (9) and the structures discussed above are evident.

Finally, in Figures (10) we present the simulations of the phase shift two qubit gate in Ref. paper1  with an intensity noise. In this case we use two exciton states ($`|E_L^iE_L^j`$ and $`|GG`$) and a two photon interaction Hamiltonin similar in structure to (1) but with Rabi frequency $`\mathrm{\Omega }_{eff}=2\mathrm{}\mathrm{\Omega }^2/\delta `$ (where $`\delta `$ is the laser detuning we need to avoid single photon transition and create two exciton states).

In Figure (10) (A) where the second regime (where the fidelity decreases) is present but less evident and seems to be compressed in the slowly varying random fluctuations zone ($`10T_{ad}/T_n30`$); for great $`T_{ad}/T_n`$ ratio the fidelity decisively increases as in the previous figures (Figure (10) (B)). Moreover, we note that the fidelity between the adiabatic final states with and without noise is high and close to $`1`$ even if we have chosen the adiabatic parameter ($`\mathrm{\Omega }_{eff}T`$) smaller than the one used for single qubit gates. These can be explained as consequence of the fact that in the imperfect adiabatic evolution unwanted states gets populated. On the other hand the effect of the fluctuating noise is to induce undesired transition as well. These two effects are superposed and with smaller adiabatic parameter the effect of the noise seems to be less important.

## IV Conclusions

We numerically studied the robustness of a non-abelian holonomic quantum gate against stochastic errors in control parameters. The robustness of logical gates show three regimes upon the variations of the noise correlation time $`T_n`$. A possible interpretation of these regimes can be given on the basis of the the geometric (i.e. solid angle swept in the parameter space) dependence of the holonomic operator. For fast random varying fluctuations we have a good robustness of the holonomic gate since, as argued in other papers preskill ; ellinas , the fluctuations during the loop tend to cancel out. For random varying fluctuations in the intermediate regime the holonomic gates are significantly corrupted because the fluctuations strongly deform the parameter loop. For slowly random varying fluctuations the performance improves again. This fact is not surprising as it may seem, indeed the loop in the parameter space turns out in this case to be simply is shifted rather than deformed; then similar solid angles are swept. Our analysis suggests that the main noise source is given by fluctuations in the intensity of the control parameters whereas the phase fluctuations do not seem to sizeably affect the gate studied. The effect of the noise decreases with the variance of the intensity of the fluctuations and for $`\delta \mathrm{\Omega }/\mathrm{\Omega }=0.01`$ we have average fidelity close to $`1`$. A first comparison shows that holonomic and dynamical gates have comparable performance in the $`T_{ad}/T_n1`$ region.

We performed similar simulations for two single qubit gates and for a two qubit gate in order to complete the set of universal quantum gates. For the single qubit gates we obtain similar results. For the two qubit gate the three regimes descibed above are present but less evident.

We believe that our analysis and conclusions should be extended rather easily to different sort of system proposed for HQC: for example, it definitely extends to the model proposed in ref. duan  since the involved holonomic structure is isomorphic to the one here analyzed. The general features of our results should hold also more general situations, since they apparently do not rely on the detailed features of Hamiltonian (1), rather on the general structure of holonomic computations. A related, though logically distinct, issue is the robustness of HQC against environmental decoherence florian . This is clearly an important subject to be addressed in future investigations.

P.Z. gratefully acknowledges financial support by Cambridge-MIT Institute Limited and by the European Union project TOPQIP (Contract IST-2001-39215)