text:
# Testing cosmological defect formation in the laboratory*

## 1 Introduction

For quite a few years, my research has been mainly on the interface between particle physics and cosmology. So it may be quite surprising that I am talking at a conference on Vortex Matter! But there are good reasons. New connections have been forged in the last few years with condensed matter physics — connections in which vortices play a central role, and that form the theme of the ESF Programme on Cosmology in the Laboratory. I want to explain how this came about. (I regret that Grisha Volovik, who co-chairs that Programme with me, and who hoped to be here too, was unable to come.)

Our present understanding of fundamental particle physics leads us to believe that very early in its history, the Universe underwent a series of phase transitions. The full symmetry of the underlying theory is apparent only at extremely high energies. As the Universe expands and cools, the symmetry is progressively broken. For example the electroweak symmetry $`\mathrm{SU}(2)\times \mathrm{U}(1)`$ is manifest only at energies above a few hundred GeV; below that the symmetry is broken by the Higgs mechanism to the $`\mathrm{U}(1)`$ of electromagnetic gauge invariance. At even higher energies, the symmetry may be larger still. There may be a GUT transition at an energy scale of about $`10^{15}`$ GeV, above which the strong, weak and electromagnetic interactions are all combined in a grand unified theory with a symmetry group such as $`\mathrm{SU}(5)`$ or $`\mathrm{SO}(10)`$. If, as is widely believed, the underlying theory is a superstring or M-theory, the implied supersymmetry must be broken at a transition scale of perhaps a few TeV.

But how can we test these ideas? The energy scales are far too high to be accessible to accelerator experiments. The only place such energies are found is the early Universe, in the first fraction of a second after the Big Bang. But of course the early Universe was opaque; we have no direct observational access to it. We have to look for surviving traces of these very early events. There are various possibilities, but one in particular that I want to talk about. A common feature of symmetry-breaking phase transitions is the formation of topological defects of one kind or another: monopoles, cosmic strings or domain walls. All these have analogues in condensed-matter systems: hedgehogs, vortices or flux tubes, and solitons.

Because of their topological stability, defects may have survived long enough to be observable . Indeed, monopoles and domain walls, if formed, could have survived all too well, and already be in conflict with observation, though there are possible ways out of this. Cosmic strings, however, are attractive from a cosmological perspective. They might help, for example, in explaining baryogenesis — creating the observed matter–antimatter asymmetry of our Universe — or magnetogenesis — seeding the magnetic field of galaxies.

To make definite predictions that can be tested against astronomical observation, we need to answer several questions:

* What kinds of defects can be formed in the hypothesized phase transitions?
* How many defects would be formed at the phase transition?
* How would they evolve as the Universe expands?
* How would they interact with matter and radiation to generate observable signatures?

We have what we believe are good answers to these questions, based on various assumptions and approximations. But can we rely on them, when our methods are untested?

Cosmology is an exciting subject, but it suffers from one major shortcoming: we have no means of conducting controlled experiments, so we cannot directly test our calculational methods. It is here that the analogy with condensed matter comes in. The mathematical descriptions of topological defects in particle-physics models and in condensed-matter systems are often very similar. The last of the four questions above is a purely cosmological one, but for the others analogies with condensed matter may be very instructive. Our methods of computing, for example, the number of defects formed at a cosmological phase transition can also be applied to a transition in a suitably chosen condensed-matter system, on which we can do real experiments, an idea first suggested by Zurek . Following through this idea has led to some very exciting developments in condensed-matter physics.

## 2 Cosmic strings

The simplest model in which strings or vortices are formed is described by a complex scalar field $`\varphi `$ with $`\mathrm{U}(1)`$ symmetry, $`\varphi \varphi \mathrm{}^{\mathrm{i}\alpha }`$. This may be a global symmetry or a local gauge symmetry. In the latter case, $`\varphi `$ interacts with a gauge potential $`A_\mu `$ transforming according to $`A_\mu A_\mu \frac{1}{e}_\mu \alpha `$. The scalar potential $`V`$ must be $`\mathrm{U}(1)`$-invariant. If we choose

$$V(\varphi )=\lambda (|\varphi |^2\eta ^2)^2,$$

(often called the Mexican hat potential), this is the Abelian Higgs model, the relativistic analogue of the Ginzburg–Landau model. In this case, $`\varphi =0`$ is a maximum of the potential, so the symmetry is broken. There is a degenerate ground state, labelled by the phase angle $`\alpha `$: $`\varphi =\eta \mathrm{}^{\mathrm{i}\alpha }`$.

This model (with or without the gauge field) exhibits a phase transition. Above a critical temperature $`T_\mathrm{c}\eta `$ (I use units in which $`c=\mathrm{}=k_\mathrm{B}=1`$) there are large fluctuations in $`\varphi `$ about a mean value of zero. As the system is cooled through the transition, $`\varphi `$ falls into the trough of the potential and acquires a nonzero average value. In so doing it has to choose a phase $`\alpha `$. But in a large system such as the Universe, there is no reason why this choice should be the same everywhere; $`\alpha `$ will vary randomly in space. Indeed, there can clearly be no correlation between the directions of $`\alpha `$ in regions beyond the causal horizon, which have had no previous causal contact.

Now it may happen that if we traverse some large loop in space, the value of $`\alpha `$ will change by $`2\pi `$ (or some multiple of $`2\pi `$). In that case, somewhere inside the loop $`\varphi `$ must go through zero. Indeed, it must vanish along a curve that threads its way through the loop; this is the cosmic string. It is the precise analogue of the Abrikosov vortex in the Ginzburg–Landau model.

Because $`\varphi `$ has to climb over the central hump of the potential, there is excess energy trapped on the string. In the local-symmetry case, this yields a string tension of order $`\eta ^2`$. When the symmetry is global, the energy per unit length of an isolated string is logarithmically divergent, but in a real system the divergence is cut off at the average string separation (or the system size). The strings are topologically stable. Strings will tend to shorten with time, under the effect of the string tension, and small closed loops of string may shrink and disappear, but they cannot break (except by a very exotic and unlikely process involving the formation of a pair of black holes). Because $`\alpha `$ is uncorrelated over large distances, the phase transition will generate a random tangle of cosmic string. The density of string will tend to decrease with time, but some string is likely to survive a long time, long enough to have observable consequences.

This $`\mathrm{U}(1)`$ model is merely a simple example, but strings also appear in many more realistic models, in particular in grand unified theories. Their dynamics and likely evolution are not usually very different, however — though there is one caveat: some strings can support persistent currents carried by fermions trapped in zero modes on the string, strongly influencing their behaviour.

## 3 The Zurek predictions

As the Universe cools through the relevant transition, we expect a random tangle of string to be formed, with some characteristic scale $`\xi _{\mathrm{str}}`$. In other words, in any randomly chosen volume $`\xi _{\mathrm{str}}^3`$, we expect to find on average a length $`\xi _{\mathrm{str}}`$ of string. The question is: what is it that determines this scale? Of course, it must be related in some way to the correlation length $`\xi `$ of the scalar field. But in the neighbourhood of the transition, $`\xi `$ is changing very rapidly. Indeed, at a second-order transition, the equilibrium correlation length $`\xi _{\mathrm{eq}}`$ goes to infinity at the transition temperature $`T_\mathrm{c}`$. So this raises the question: at what temperature, or what time, should we equate $`\xi _{\mathrm{str}}`$ to $`\xi `$?

Zurek  has provided an answer to this question based on a causality argument (see also ). When the system goes through a real transition, at a finite rate, it is clear that $`\xi `$ cannot become infinite. There is a maximum speed, $`c`$, with which correlations in the phase of the scalar field can propagate. In the relativistic case, this is the speed of light. In a non-relativistic system it is some characteristic speed of the system, for example the speed of second sound in superfluid helium-4. Correlations can never extend beyond a finite range, determined by a balance between the relaxation rate of the scalar field and the rate $`\dot{T}/T`$ of the transition.

Zurek’s argument may be paraphrased thus: As the system cools, $`\xi `$ more or less keeps up, so long as it can do so, with the equilibrium correlation length, $`\xi _{\mathrm{eq}}`$. But once $`\xi _{\mathrm{eq}}/t`$ becomes larger than $`c`$, it can no longer do so. From then on, $`\xi `$ does not change much until the point after the transition when it again becomes equal to the decreasing $`\xi _{\mathrm{eq}}`$. That is the time, now often called the Zurek time, $`t_\mathrm{Z}`$, at which we should identify $`\xi _{\mathrm{str}}`$ with $`\xi `$. An almost equivalent statement (up to a factor of order 1) is that $`t_\mathrm{Z}`$ is the time at which correlations starting from zero at the transition and propagating with speed $`c`$, can reach the distance $`\xi _\mathrm{Z}=\xi _{\mathrm{eq}}(t_\mathrm{Z})`$.

This argument leads to a definite prediction for the defect density, $`l`$, that is, the average length of string per unit volume. At $`t_\mathrm{Z}`$, we expect $`l(t_\mathrm{Z})=k/\xi _\mathrm{Z}^2`$, where $`k`$ is a constant roughly of order one. Numerical simulations  suggest that $`k`$ is actually somewhat less than one, perhaps of order $`0.1`$.

## 4 Tests in Liquid Crystals and Helium-4

Zurek  originally suggested testing these predictions by looking at a rapid quench in superfluid helium-4. Starting at high pressure in the normal phase just above the ‘lambda line’, a rapid redution of pressure takes the sample through the transition into the superfluid phase.

The first experiments, however , were in thin samples of nematic liquid crystals, where networks of defects were seen to be formed when the system was rapidly cooled from the normal to the nematic phase. The nematic transition is first-order, so Zurek’s argument as given above is not directly applicable. In that case, what one expects  is that $`\xi _{\mathrm{str}}`$ is approximately (a few times) the mean distance between nucleation centres of bubbles of the new phase. This indeed seems to be the case.

The helium-4 experiment, which provides a more direct test of Zurek’s predictions, has now been performed, twice, by a group in Lancaster, using an apparatus comprising a small chamber filled with helium whose sides were formed of bellows so that it could be rapidly expanded. The object of the exercise was to detect the presence of any vortices formed during the transition by monitoring the absorption of second sound, which is strongly attenuated by vortices. The first experiment  did apparently show attentuation of the signal after the quench, falling off exponentially with a timescale of the order of a hundred milliseconds, as the vortices gradually disappeared, and compatible in magnitude with Zurek’s prediction for the number of vortices. But it was inconclusive for various reasons. Firstly, the detection equipment was swamped and unable to measure the attenuation until about 50 ms after the quench, so the inferred vortex density had to be extrapolated back to the Zurek time. Secondly, there were other possible sources of vorticity. In particular, vortices might have been nucleated by hydrodynamic effects at the bellows, and a particular concern attached to the capillary tube used to fill the chamber which was closed at its outer end, so that fluid was injected into the chamber on each expansion.

For these reasons, the Lancaster group decided to build an improved apparatus, in which the chamber was made much smoother, and in particular the injection tube was closed at the point of entry into the chamber, to minimize the risk of extraneous vortex formation. However, results with this improved apparatus  rather disappointingly failed to reveal any vortex formation. I shall discuss possible reasons for this later.

## 5 Tests in Helium-3

Meanwhile, experiments were performed using helium-3 at two different laboratories, Helsinki and Grenoble (with collaboration from Lancaster). There are considerable advantages in using the lighter isotope. One is that the correlation length is much larger, say 40–100 nm as compared with less than 1 nm in $`{}_{}{}^{H}4`$. This means that a continuum description, of Ginzburg–Landau type, is a much better approximation. It also means that relatively speaking vortex formation requires much more energy, so that extraneous vortex formation is less likely. But perhaps the most important difference is that, because $`{}_{}{}^{H}3`$ is an excellent neutron absorber, one can perform a temperature- rather than pressure-driven quench. Both experiments made use of this technique. The absorption of a neutron in a sample in the superfluid $`{}_{}{}^{H}3`$-B phase heats up a small region above the critical temperature. It then cools rapidly, in a period of about 1 $`\mu `$s, back through the transition into the superfluid phase. One expects the formation of a random tangle of vortices during this process.

In other respects the experiments were very different. The Helsinki experiment  used a sample at a temperature not far below the transition temperature, in a rotating cryostat. If the rotation speed $`v`$ is less than some limit, no vortices are formed spontaneously at the walls, so the superfluid component is actually completely stationary, while the normal component is rotating with the container. The relative velocity between the two components means that any vortex formed following neutron absorption will be subject to the transverse Magnus force. If a vortex loop is big enough and appropriately oriented, the effect will be to expand it until it meets the walls of the container, after which the ends migrate to the top and bottom surfaces, and the vortex joins a central cluster parallel to the rotation axis.

Detection of these captured vortices is made possible by another of the key features of $`{}_{}{}^{H}3`$, its non-zero nuclear spin, which makes it possible to use nuclear magnetic resonance. In the NMR trace one can actually see each individual vortex joining this central cluster. It turns out that the Zurek scenario leads to a very simple prediction for the dependence of the number of vortices captured per neutron event on $`v`$. There is a critical velocity $`v_\mathrm{c}`$ below which no vortices are captured. When $`v>v_\mathrm{c}`$,

$$n=\gamma \left\{\left(\frac{v}{v_\mathrm{c}}\right)^31\right\},$$

where $`\gamma `$ is a constant. Remarkably the entire dependence on the pressure, the bulk temperature and the magnetic field is contained in the single parameter $`v_\mathrm{c}`$. The results clearly confirm this: plotted against $`v^3`$, they fit well to straight lines, with a common intercept at $`\gamma `$. Moreover the dependence of $`v_\mathrm{c}`$ on the bulk temperature also fits the prediction.

The Grenoble experiment  was complementary. Essentially it involved calorimetry. They used a sample at a much lower temperature, but in a non-rotating container, and sought to measure the energy release every time a neutron was absorbed. The absorption reaction is

$$n+{}_{}{}^{H}3p+{}_{}{}^{H}3+764\mathrm{keV}.$$

The total energy released in the form of quasiparticles was measured, and found to peak in the range 575–650 keV, depending on the pressure. Some energy, around 50 keV, is also released in the form of ultraviolet radiation. However, there is a very clear overall deficit, which is attributed to vortices. It is very difficult to think of any other mechanism of energy loss. Moreover, the magnitude and the dependence on pressure are entirely consistent with Zurek’s prediction.

## 6 Tests in Superconductors

From the point of view of analogy with cosmological phase transitions, tests in superconductors are particularly important, because they provide an example of symmetry breaking in a local gauge theory.

At least two such experiments have in fact been performed, by a group at Technion. In the first , they used a 1 cm2 thin film of the high-temperature superconductor, YBCO. The film was heated to above $`T_\mathrm{c}`$ by optical irradiation, and then allowed to cool back through the transition. To determine the number of flux tubes formed during this process, the experimenters measured the total flux using a SQUID close to the sample, but separated from it by a mylar sheet. To be precise, this process does not measure the number of fluxons, but only the net total flux, the number of fluxons minus the number of antifluxons. A straightforward application of Zurek’s technique would lead to an estimate for the net flux of about 100, whereas Carmi and Polturak’s results  were consistent with zero, with an error of about 10. Moreover, they argue that because of the different behaviour of a gauge theory, the actual prediction should be higher, around 104, in which case there is a very clear contradiction. However, this prediction does not seem to be on very firm ground.

The same group have done another, similar experiment involving a loop of semiconductor with a large number of Josephson junctions in series. In that case they do find definite evidence  of the trapping of magnetic flux during the cooling process. Here the loop becomes superconducting before being linked by the junctions, so one might expect the phases to be random. The experiments actually show a larger flux than predicted on that basis, suggesting that the phase differences are not uniformly distributed, though it is not entirely clear why.

Experiments have been suggested  using an annular Jospehson junction, in which cooling can lead to the trapping of radial flux lines. A group in Salerno has already done some experiments with an apparatus of this kind, but the experimental parameters were not optimized for this particular type of measurement: the predicted average number of trapped quanta is less than unity. Future experiments along these lines should yield interesting results.

## 7 Theoretical Interpretation

Zurek’s predictions of the numbers of defects formed are based on an admittedly rather crude argument, and it is certainly not surprising to find that it needs modification.

At first sight, the difference between the behaviour of helium-3 and helium-4 may seem surprising, though in fact theoretically they are very different types of system. There is a critical region below $`T_\mathrm{c}`$ but above the so-called Ginzburg temperature, within which there is a large population of transient thermal loops. This region is much wider in helium-4 than in helium-3. Work by Karra and Rivers  has shown that this has a marked effect on the Zurek predictions. On the other hand, there may be a more immediate reason for the failure of the helium-4 experiment to detect any vortices: they may simply not last long enough to have been seen. Vorticity does decay with time, with a lifetime that is not well known. The authors  estimated the lifetime on the basis of measurements of decay of turbulent vorticity, but this may well be inapplicable . It would be very desirable to measure the vorticity at a much earlier time after the transition. The situation should be clarified by future experiments .

Understanding the results in superconductors also presents a challenge. It is puzzling that the experiment with a thin film of YBCO saw nothing. One of the problems is that there is much more theoretical uncertainty about the predicted defect density in the case of a local gauge theory. There is a competing mechanism for the formation of defects, namely the thermal fluctuations in the magnetic field. Hindmarsh and Rajantie  have shown that that this mechanism yields a very different density and distribution of defects. It is unclear which mechanism will dominate. Neither mechanism seems to support the estimate that led to the conclusion of a disagreement by a factor of 1000 in the experiment with a thin superconducting film, but there remains at least a factor of 10 that is hard to explain.

The Josephson junction results seem to be positive, but if anything show too large a trapped flux. Why there should be this difference between the two experiments with superconductors is again something of a puzzle. More experiments along these lines too are very desirable.

## 8 Conclusions

The results so far are distinctly confusing. By far the best evidence in favour of the Zurek scenario comes from the experiments in helium-3, both of which yielded positive results that are hard to interpret in any other way. Results in nematic liquid crystals are also positive, but not a direct test of the Zurek scenario. In superconductors, it is puzzling that the experiment with a thin YBCO film saw no fluxons, while more than expected were seen in the Josephson-junction array.

It does seem that vortex formation in a rapid transition is now well established, but there is still a lot of work to do to understand the details.

## Acknowledgements

I am very grateful to the organziers of this meeting for inviting me to speak. The work reported here was supported by the European Science Foundation under the COSLAB Programme (Cosmology in the Laboratory).


format: text

text:
# 1 Introduction

## 1 Introduction

In recent years, several important insights have been gained into the physical principles of protein folding . Still, in terms of quantitative predictions, it is clear that it would be extremely useful to be able to perform more realistic folding simulations than what is currently possible. In fact, most models that have been used so far for statistical-mechanical simulations of folding rely on one or both of two quite drastic approximations, the lattice and Gō  approximations.

The reason that lattice models have been used to study basics of protein folding is partly computational, but also physical — on the lattice, it is known what potential to use in order for stable and fast-folding sequences to exist (a simple contact potential is sufficient). How to satisfy these criteria for off-lattice chains is, by contrast, largely unknown, and therefore many current off-lattice models  use Gō-type potentials  where non-native interactions are ignored. The use of the Gō approximation has some support from the finding that the native structure is a determinant for folding kinetics . However, it is an uncontrolled approximation, and it is, of course, useless when it comes to structure prediction, as it requires prior knowledge of the native structure.

In this paper, we discuss an off-lattice model that does not follow the Gō prescription. Using this model, we perform extensive folding simulations for a small helical protein. The force field of the model is simple and based on hydrogen bonds and effective hydrophobicity forces (no explicit water). There exist other non Gō-like models with more elaborate force fields that have been used for structure prediction with some success . However, it is unclear what the dynamical properties of these models are.

The original version of our model was presented in Ref.  and has three types of amino acids: hydrophobic, polar and glycine. This version was applied to a designed three-helix-bundle protein with 54 amino acids . For a suitable relative strength of the hydrogen bonds and hydrophobicity forces, it was found that this sequence does form a stable three-helix bundle, except for a twofold topological degeneracy, and that its folding transition is first-order-like and coincides with the collapse transition (the parameter $`\sigma `$ of Ref.  is zero).

Here, we extend this model from three to five amino acid types, by taking alanine to be intermediate in hydrophobicity between the previous two hydrophobic and polar classes, and by introducing a special geometric representation for proline, which is needed to be able to mimic the helix-breaking property of this amino acid. Otherwise, the model is the same as before. The modified model is tested on a real three-helix-bundle protein, the 10–55-amino acid fragment of the B domain of staphylococcal protein A. The structure of this protein has been determined by NMR , and an energy-based structure prediction method has been tested on the sequence . The folding properties have been studied too, both experimentally  and theoretically . In particular, this means that we can compare the behavior of previous Gō-like models to that of our more realistic model.

## 2 Materials and Methods

### 2.1 Geometry

Our model is an extension of that introduced in Ref. . It uses three different amino acid representations: one for glycine, one for proline and one for the rest. The non-glycine, non-proline representation is illustrated in Fig. 1a, and is identical to that of hydrophobic and polar amino acids in the original model. The three backbone atoms N, Cα and C are all included, whereas the side chain is represented by a single atom, a large Cβ. The remaining two atoms, H and O, are used to define hydrogen bonds. The representation of glycine is the same except that Cβ is missing.

The representation of proline is new compared to the original model. The side chain of proline is attached to the backbone not only at Cα, but also at N. A well-known consequence of this is that proline can act as a helix breaker. For the model to be able to capture this important property, we introduce a special representation for proline, which is illustrated in Fig. 1b. It differs from that in Fig. 1a in two ways: first, the Ramachandran angle $`\varphi `$ is held constant, at $`65^{}`$; and second, the H atom is replaced by a side-chain atom, Cδ. This more realistic representation of proline is needed when studying the protein A fragment which has one proline at each of the two turns.

All amino acids except proline have the Ramachandran torsion angles $`\varphi `$ and $`\psi `$ (see Fig. 1a) as their degrees of freedom, whereas $`\psi `$ is the only degree of freedom for proline. All bond lengths, bond angles and peptide torsion angles ($`180^{}`$) are held fixed. Numerical values of the bond lengths and bond angles can be found in Ref.  and Fig. 1b.

The helix-breaking property of proline manifests itself clearly in the shape of the $`\psi `$ distribution for amino acids that are followed by a proline in the sequence (with the proline on their C side). Helical values of $`\psi `$ are suppressed for such amino acids. This is illustrated in Fig. 2a, where the peak on the left corresponds to $`\alpha `$-helix. From Fig. 2b, it can be seen that the model shows a qualitatively similar behavior.

### 2.2 Force Field

Our energy function

$$E=E_{\text{loc}}+E_{\text{sa}}+E_{\text{hb}}+E_{\text{col}}$$
(1)
is composed of four terms. The first two terms $`E_{\text{loc}}`$ and $`E_{\text{sa}}`$ are local $`\varphi ,\psi `$ and self-avoidance potentials, respectively (see Ref. ). The third term is the hydrogen-bond energy $`E_{\text{hb}}`$, which is given by

$`E_{\text{hb}}`$ $`=`$ $`ϵ_{\text{hb}}{\displaystyle \underset{ij}{}}\left[5\left({\displaystyle \frac{\sigma _{\text{hb}}}{r_{ij}}}\right)^{12}6\left({\displaystyle \frac{\sigma _{\text{hb}}}{r_{ij}}}\right)^{10}\right]v(\alpha _{ij},\beta _{ij})`$ (2)
$`v(\alpha _{ij},\beta _{ij})`$ $`=`$ $`\{\begin{array}{cc}\mathrm{cos}^2\alpha _{ij}\mathrm{cos}^2\beta _{ij}\hfill & \alpha _{ij},\beta _{ij}>90^{}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (5)
where $`i`$ and $`j`$ represent H and O atoms, respectively, and where $`r_{ij}`$ denotes the HO distance, $`\alpha _{ij}`$ the NHO angle, and $`\beta _{ij}`$ the HOC angle.

The last term in Eq. (1), the hydrophobicity or collapse energy $`E_{\text{col}}`$, has the form

$$E_{\text{col}}=ϵ_{\text{col}}\underset{i<j}{}\mathrm{\Delta }(s_i,s_j)\left[\left(\frac{\sigma _{\text{col}}}{r_{ij}}\right)^{12}2\left(\frac{\sigma _{\text{col}}}{r_{ij}}\right)^6\right],$$
(6)
where the sum runs over all possible CβCβ pairs and $`s_i`$ denotes amino acid type. To define $`\mathrm{\Delta }(s_i,s_j)`$, we divide the amino acids into three classes: hydrophobic (H; Leu, Ile, Phe), alanine (A; Ala) and polar (P; Arg, Asn, Asp, Gln, Glu, His, Lys, Pro, Ser, Tyr). Cys, Met, Thr, Trp and Val do not occur in the sequence studied. There are then six kinds of CβCβ pairs, and the corresponding $`\mathrm{\Delta }(s_i,s_j)`$ values are taken to be

$$\mathrm{\Delta }(s_i,s_j)=\{\begin{array}{cc}\hfill 1& \text{for HH and HA pairs}\hfill \\ \hfill 0& \text{for HP, AA, AP and PP pairs}\hfill \end{array}$$
(7)
The main change in the force field compared to Ref.  is that alanine forms its own hydrophobicity class, besides the previous two hydrophobic and polar classes. Alanine is taken as intermediate in hydrophobicity, meaning that there is a hydrophobic interaction between HA pairs but not between AA pairs. In addition, the interaction strength $`ϵ_{\text{col}}`$ is increased slightly, from 2.2 to 2.3.The energy unit is dimensionless and such that $`kT_\text{c}=0.62`$, $`T_\text{c}`$ being the collapse temperature (see Sec. 3). Finally, in the self-avoidance potential, the Cδ atom of proline is assigned the same size as Cβ atoms. Otherwise, the entire force field, including parameter values, is exactly the same as in Ref. .

With these changes in geometry and force field, we end up with five different amino acid types in the new model. First, we have hydrophobic, alanine and polar which share the same geometric representation but differ in hydrophobicity, and then glycine and proline with their special geometries.

In this paper, we test this model on the 10–55-amino acid fragment of the B domain of staphylococcal protein A. Calculated structures are compared to the minimized average NMR structure  with PDB code 1bdd. Throughout the paper, this structure is referred to as the native structure.

As a first test of our model, two different fits to the native structure were made. The first fit is purely geometrical. Here, we simply minimized the root-mean-square deviation (rmsd) from the native structure, $`\delta `$ (calculated over all backbone atoms). This was done by using simulated annealing, and the best result was $`\delta =0.14`$ Å. In the second fit, we took into account the limitations imposed by the first three terms of the potential, by minimizing the function

$$\stackrel{~}{E}=E_{\text{loc}}+E_{\text{sa}}+E_{\text{hb}}+\kappa \underset{i}{}(𝐫_i𝐫_i^0)^2,$$
(8)
where $`\kappa =1`$ Å-2 and $`\{𝐫_i^0\}`$ denotes the structure obtained from the first fit. The minimum-$`\stackrel{~}{E}`$ structure had $`\delta =0.32`$ Å. These results show that our model, in spite of relatively few degrees of freedom, permits a quite accurate description of the real structure.

### 2.3 Numerical Methods

To simulate the thermodynamic behavior of this model, we use simulated tempering , which means that the temperature is a dynamical variable (for details, see Refs. ). The temperature update is a standard Metropolis step. Our conformation updates are of two different types: the simple non-local pivot move where a single torsion angle is turned, and the semi-local biased Gaussian step proposed in Ref. . The latter method works with the Ramachandran angles of four adjacent amino acids. These are turned with a bias toward local rearrangements of the chain. The degree of bias is governed by a parameter $`b`$. In our thermodynamic simulations, we take $`b=10`$ (rad/Å)2, which gives a strong bias toward deformations that are approximately local .

Figure 3 shows the evolution of the energy in a simulated-tempering run that took about two weeks on an 800 MHz processor. Data corresponding to all the different temperatures are shown (eight temperatures, ranging from $`kT=0.54`$ to $`kT=0.90`$). We see that there are many independent visits to low-energy states, which is necessary in order to get a reliable estimate of the relative populations of the folded and unfolded states. To test the usefulness of the semi-local update, we repeated the same calculation using pivot moves only. The difference in performance was not quantified, but it was clear that the sampling of low energies was less efficient in the run relying solely on pivot moves.

For our kinetic simulations, we do not use the pivot update but only the semi-local method. The parameter $`b`$ is taken to be 1 (rad/Å)2 in the kinetic runs, which turned out to give an average change in the end-to-end vector squared of about 0.5 Å2.

## 3 Results and Discussion

### 3.1 Thermodynamics

We begin our study of the model defined in Sec. 2 by locating the collapse transition. In Fig. 4, we show the radius of gyration (calculated over all backbone atoms) against temperature for both the protein A sequence and three random sequences with the same length and composition. The random sequences were generated keeping the two prolines of the protein A sequence fixed at their positions, one at each turn. The remaining 44 amino acids were randomly reshuffled.

Naively, one may expect these sequences to show similar collapse behaviors, since the composition is the same. However, the protein A sequence turns out to collapse much more efficiently than the random sequences (see Fig. 4). The native structure has a radius of gyration of 9.25 Å, which is significantly smaller than one finds for the random sequences in this temperature range. The specific heat (data not shown) has a pronounced peak in the region where the collapse occurs. Taking the maximum as the collapse temperature $`T_\text{c}`$, we obtain $`kT_\text{c}=0.62`$ for the protein A sequence.

The chain collapse is not as abrupt for the protein A sequence as for the designed sequence studied in Ref. . This is not surprising, as that sequence has a hydrophobicity pattern that fits its native structure perfectly. The protein A sequence does not have a fully perfect hydrophobicity pattern, but still the collapse behavior is highly cooperative, as can be seen from the comparison with the random sequences.

Next, we turn to the structure of the collapsed state. As a measure of similarity with the native structure, we use

$$Q=\mathrm{exp}(\delta ^2/100\mathrm{\AA }^2),$$
(9)
where $`\delta `$, as before, denotes rmsd. An alternative would be to base the similarity measure on the number of native contacts present, rather than rmsd. The problem with such a definition is that it does not provide an efficient discrimination between the two possible topologies of a three-helix bundle  — the third helix can be either in front of or behind the U formed by the first two helices. This problem is avoided by using rmsd.

In Fig. 5a, we show the free-energy profile $`F(Q)`$ in the collapsed phase at $`kT=0.54`$.

We see that there is a broad minimum at $`Q0.8`$–0.9, with two distinct local minima at $`Q=0.78`$ and $`Q=0.90`$, respectively. Both these minima correspond to the native overall topology. There is also a minimum at $`Q=0.50`$, which corresponds to the wrong topology. The $`Q=0.50`$ minimum is more narrow and slightly higher, so the native topology is the favored one. However, it should be stressed that it is difficult to discriminate between the two topologies using a pairwise additive potential (see Sec. 3.4). To be able to do that in a proper way, it is likely that one has to include multibody terms and/or more side-chain atoms in the model.

The main difference between the two minima at $`Q=0.78`$ and $`Q=0.90`$ lies in the shape and orientation of helix III, which comprises amino acids 41–55 in the native structure. At the $`Q=0.78`$ minimum, there tends to be a sharp bend in this segment, and the amino acids before the bend, 41–44, are disordered rather than helical. The remaining amino acids, 45–55, tend to make a helix, but its orientation differs from that in the native structure. Relative to the $`Q=0.90`$ minimum, where helix III is much more native-like, we find that the $`Q=0.78`$ minimum is entropically favored but energetically disfavored. The separation in energy between these minima is probably underestimated by our model. There is, for example, a stabilizing electrostatic interaction between helices I and III in the native structure (Glu16-Lys50), which should favor the $`Q=0.90`$ minimum but is missing in our model.

Also shown in Fig. 5a is the result for one of the random sequences. The probability of finding this sequence in the vicinity of the native structure is, not unexpectedly, very low. The same holds true for the other two random sequences too (data not shown).

To extract representative conformations for the collapsed state, we used simulated annealing followed by a conjugate-gradient minimization. Using this procedure, a large set of low-temperature Monte Carlo conformations were quenched to zero temperature. In Fig. 5b, we show the quenched conformations with lowest energy in a $`Q,E`$ scatter plot. Our minimum-energy structure is found at $`Q=0.44`$, corresponding to $`\delta =9.1`$ Å. However, our thermodynamic calculations show that this conformation is not very relevant, in spite of its low energy. If we restrict ourselves to conformations with the native-like and thermodynamically most relevant topology, then the lowest energy is at $`Q=0.97`$, corresponding to $`\delta =1.8`$ Å. This conformation is shown in Fig. 6 along with the native structure. It is worth noting that the $`Q=0.44`$ and $`Q=0.97`$ minima both were revisited in independent runs.

These results can be compared with those of Scheraga and coworkers , who tested an energy-based structure prediction method on the same sequence. With their energy function, the global minimum was found to have an rmsd of 3.8 Å from the native structure (calculated over Cα atoms).

### 3.2 Helix Stability

Having discussed the overall thermodynamic behavior, we now take a closer look at the stability of the secondary structure and how it varies along the chain. To this end, we monitored the hydrogen-bond energy between the CO group of amino acid $`i`$ and the NH group of amino acid $`i+4`$ \[see Eqs. (2,5)\], $`e_{\text{hb}}(i)`$, as a function of $`i`$. This was done not only for the protein A sequence, but also for the corresponding three one-helix segments, which are listed in Table 1. An experimental study  of essentially the same three segments found segment III to be the only one that shows some stability on its own.

The results of our calculations are shown in Fig. 7, from which we see that the difference between the full sequence and the one-helix segments is not large in the model. However, the segments I and II definitely make less stable helices on their own than as interacting parts of the full system; they are stabilized by interhelical interactions. Furthermore, among the three one-helix segments, the model correctly predicts segment III to be the most stable one. That this segment does not get more stable as part of the full system is probably related to the observation above that helix III is distorted at the $`Q=0.78`$ minimum.

A striking detail in Fig. 7 is that the beginning of segment II is quite unstable. This can be easily understood. This segment has a flexible glycine at position 30, and the amino acids before the glycine, 24–29, are all polar, so there are no hydrophobic interactions that can help to stabilize this part.

### 3.3 Kinetics

Using the semi-local update , we performed a set of 30 kinetic simulations at $`kT=0.54`$. The runs were started from random coils. There are big differences between these runs, partly because the system, as it should, sometimes spent a significant amount of time in the wrong topology. Nevertheless, the data show one stable and interesting trend, namely, that the formation of helices was never faster than the collapse. This is illustrated in Fig. 8, which shows the evolution of the similarity parameter $`Q_0`$, the hydrogen-bond energy $`E_{\text{hb}}`$ and the radius of gyration, $`R_\text{g}`$, in one of the runs. $`Q_0`$ is defined as $`Q`$ in Eq. (9), except that it measures similarity to the optimized model structure in Fig. 6 rather than the native structure. In Fig. 8, we see that $`E_{\text{hb}}`$ converges slowly, whereas the collapse occurs relatively early.

Now, at a first glance, it may seem easy to make the helix formation faster by simply increasing the strength of the hydrogen bonds. Therefore, it is important to note that the hydrogen bonds cannot be made much stronger without making the ground state non-compact and thus destroying the three-helix bundle . This means that the conclusion that the collapse is at least as fast as helix formation holds for any reasonable choice of parameters in this model.

It is interesting to compare these results to those of Zhou and Karplus , who studied the same protein using a Gō-type potential and observed fast folding when the Gō forces were strong. Under these conditions, the helix formation was found to be fast, whereas the collapse was the rate-limiting step.

However, a Gō-like model ignores a large fraction of the interactions that drive the collapse, which can make the collapse artificially slow. In a recent Gō model study , this problem was addressed by eliminating backbone terms from the potential until a reasonable helix stability was achieved. No such calibration was carried out in Ref. . This may explain why these authors find a behavior that our model cannot reproduce.

Let us finally mention that we also performed the same type of kinetic simulations for the designed sequence studied in Ref.  which, as discussed earlier, has a very abrupt collapse transition. It turns out that $`E_{\text{hb}}`$ and $`R_\text{g}`$ evolve in a strongly correlated manner in this case. So, the helix formation and collapse occur simultaneously for this sequence.

### 3.4 Fine-tuning?

In Sec. 3.1, we discussed the relative weights of the two possible overall topologies, which is a delicate issue. What changes are needed in order for the model to more strongly suppress the wrong topology? Is it necessary to change the form of the energy function, or would it be sufficient to fine-tune the interaction matrix $`\mathrm{\Delta }(s_i,s_j)`$ in Eq. (6)?

One way to do such a fine-tuning of $`\mathrm{\Delta }(s_i,s_j)`$ would be to maximize $`Q^{}`$, where $`Q`$ is the similarity parameter and $`^{}`$ denotes a thermodynamic average restricted to compact conformations ($`R_\text{g}<10`$ Å say). This is essentially the overlap method of Ref. . The gradient of the quantity $`Q^{}`$ can be written as

$$\frac{Q^{}}{\mathrm{\Delta }(s_i,s_j)}=\frac{ϵ_{\text{col}}}{kT}\left(QX^{}Q^{}X^{}\right),$$
(10)
where $`X`$ is a sum of Lennard-Jones terms, $`(\sigma _{\text{col}}/r_{ij})^{12}2(\sigma _{\text{col}}/r_{ij})^6`$, over all possible CβCβ pairs of type $`s_i,s_j`$.

We calculated the $`Q,X`$ correlation in Eq. (10) for all pairs $`s_i,s_j`$ with $`\mathrm{\Delta }(s_i,s_j)=1`$ at $`kT=0.54`$, and found that $`|Q^{}/\mathrm{\Delta }(s_i,s_j)|`$ was small ($`0.15`$) for all these pairs. Hence, there is no sign that a significant increase in $`Q^{}`$ can be achieved by fine-tuning $`\mathrm{\Delta }(s_i,s_j)`$; the contact patterns seem to be too similar in the two topologies. To include more side-chain atoms and/or multibody terms in the model is likely to be a more fruitful approach.

## 4 Conclusion

We have explored a five-letter protein model with five to six atoms per amino acid, where the formation of native structure is driven by hydrogen bonding and effective hydrophobicity forces. This model, which does not follow the Gō prescription, was tested on a small but real sequence, a three-helix-bundle fragment from protein A.

Using this model, the protein A sequence was found to collapse much more efficiently than random sequences with the same composition. In the collapsed phase, we found that the native topology dominates, although the suppression of the wrong three-helix-bundle topology is not strong. Energy minimization constrained to the thermodynamically favored topology gave a structure with an rmsd of 1.8 Å from the native structure.

In our kinetic simulations, the collapse was always at least as fast as helix formation, which is in sharp contrast with previous results for the same protein that were obtained using a Gō-like Cα model . A possible explanation for the conflicting conclusions is that the Gō approximation makes the collapse artificially slow by ignoring a large fraction of the interactions driving the collapse. In our model, the conclusion that the helix formation is not faster than collapse seems unavoidable; if one tries to speed up the helix formation by increasing the strength of the hydrogen bonds, then the chain does not fold into a compact helical bundle.

The force field of our model was deliberately kept simple. In particular, the hydrophobicity potential was taken to be pairwise additive, with a simple structure for the interaction matrix $`\mathrm{\Delta }(s_i,s_j)`$ \[see Eq. (7)\]. In the future, it would be very interesting to look into the behavior of the model in the presence of multibody terms. A simpler alternative is to stick to the pairwise additive potential and fine-tune the parameters $`\mathrm{\Delta }(s_i,s_j)`$. However, the calculations in this paper give no indication that there is much to be gained from such a fine-tuning.

### Acknowledgments

This work was in part supported by the Swedish Foundation for Strategic Research.


text:
# Sensitivity to perturbations in a quantum chaotic billiard

## Abstract

The Loschmidt echo (LE) measures the ability of a system to return to the initial state after a forward quantum evolution followed by a backward perturbed one. It has been conjectured that the echo of a classically chaotic system decays exponentially, with a decay rate given by the minimum between the width $`\mathrm{\Gamma }`$ of the local density of states and the Lyapunov exponent. As the perturbation strength is increased one obtains a cross-over between both regimes. These predictions are based on situations where the Fermi Golden Rule (FGR) is valid. By considering a paradigmatic fully chaotic system, the Bunimovich stadium billiard, with a perturbation in a regime for which the FGR manifestly does not work, we find a cross over from $`\mathrm{\Gamma }`$ to Lyapunov decay. We find that, challenging the analytic interpretation, these conjetures are valid even beyond the expected range.

Hypersensitivity to initial conditions is the key ingredient of classical chaos. In quantum mechanics, its absence led to the study of other features that could be associated with the chaos of the corresponding classical system. Celebrated examples are the Gutzwiller trace formula for the quantum spectral density, the description of the spectral fluctuations by the random matrix theory and the relation of spectral correlations to transport .

In an alternative point of view, Peres  suggested that quantum dynamics should distinguish regular and irregular classical dynamics if the time evolution of an initial state for slightly different Hamiltonians are compared. That is, the sensitivity of a quantum system should be searched not by changing the initial conditions but rather by perturbing the Hamiltonian. The natural quantity for this investigation is the ability of the system to return to the initial state $`|\varphi `$ after being evolved with a Hamiltonian $`_0`$ for a period $`t`$ followed by an identical period of unitary evolution with $`_1=(_0+\mathrm{\Sigma }).`$ This defines the quantum Loschmidt echo (LE)

$$M(t)=|\varphi \left|\mathrm{exp}[\mathrm{i}_1t/\mathrm{}]\mathrm{exp}[\mathrm{i}_0t/\mathrm{}]\right|\varphi |^2.$$
(1)
The perturbation $`\mathrm{\Sigma }`$ can represent the uncontrolled degrees of freedom of an environment. As in classical chaos, the LE is related to a ‘distance’ between a perturbed and an unperturbed evolution of the same initial state.

In recent years new hints were available due to the advances in Nuclear Magnetic Resonance. The LE was measured in a many-body system of interacting spins  in a range where it is known to have spectral signatures of chaos. A striking finding was that, when interaction with the environment and residual interactions are very weak, the decay of $`M(t)`$ becomes independent of the perturbation strength. In this situation, it depends on the dynamical scales of the systems, i.e. on $`_0`$. While the complexity of the experimental system did not allow for a derivation of the characteristic time for these specific system, Jalabert and Pastawski  studied the LE in a one-body classically chaotic Hamiltonian with a perturbation represented by a long range quenched disordered potential. They have showed analytically that $`M(t)`$ may decay exponentially with a rate given by the Lyapunov exponent of the classical system. As condition, the perturbation must be quantically strong to produce statistically unpredictable changes in the quantum phase but weak enough to leave the underlying classical dynamics undisturbed.

More recently, Jacquod, Silvestrov and Beenakker  predicted a cross-over from a perturbation dependent regime to the Lyapunov one. However, this prediction is based on the strong assumption that the perturbation lives in a FGR regime; i.e. the local density of states (LDOS) is a Breit-Wigner distribution whose width $`\mathrm{\Gamma }`$ varies quadratically with the perturbation strength. In this situation, $`M(t)`$ for a wave packet and the survival probability of an unperturbed eigenstate have a decay rate given by $`\mathrm{\Gamma }.`$ Both observables would describe the same physics if the correlation between states forming the wave packet could be neglected.

Our aim is to determine whether the perturbation independent Lyapunov regime and the cross-over from a $`\mathrm{\Gamma }`$ decay are possible in a fully chaotic system with a clear semiclassical description where the presence of the perturbation is not described by the FGR. This occurs when there are strong correlations that could be related to classical structures which prevents a description in terms of a random matrix theory. This perturbation is then said to be non generic  and the LDOS can be very different from the Lorentzian analyzed in Ref. . Our positive answer in such a case opens the question of a semiclassical interpretation for the weak perturbation regime.

We consider the paradigmatic disymmetrized Bunimovich stadium billiard . It consists of a free particle inside a 2-dimensional planar region whose boundary $`𝒞`$ is shown in Fig. 1. The radius $`r`$ is taken equal to unity and the enclosed area is $`1+\pi /4`$. This system not only has a great experimental relevance , but also it is fully chaotic one by oposition to the system considered in Ref.. Besides, it can rule out the diffusive effect of disorder suspected to affect the behavior of $`M(t)`$ in a Lorentz gas. The classical dynamics is completely defined once the boundary is given. On the other hand, to address the quantum mechanics, it is necessary to solve the Helmholtz equation, $`^2\varphi _\mu =k^2\varphi _\mu `$ with appropriate boundary conditions. $`k_\mu `$ is the wave number and by setting $`\mathrm{}=2m=1`$, $`k_\mu ^2`$ results the energy. The most commonly used boundary conditions are the Dirichlet (hard walls) and the Neumann (acoustics) conditions. However, we are interested in the possibility of perturbing the quantum system without breaking the orthogonal symmetry and leaving the classical motion undisturbed . This is possible using more generalized boundary conditions:

$$\varphi (q)+\xi g(q)\frac{\varphi }{𝐧}(q)=0,$$
(2)
where $`q`$ is a coordinate along the boundary of the billiard (see Fig. 1), and $`𝐧`$ is the unit vector normal to the boundary. $`g(q)`$ is a real function and $`\xi `$ the parameter controlling the strength of the perturbation. Dirichlet boundary conditions are recovered when $`\xi =0`$ while Neumann conditions are satisfied in the limit $`\xi \mathrm{}`$. The eigenfunctions and eigenenergies for the case $`\xi =0`$ are readily obtained by using the scaling method.

In order to compute the LE in this system, a relation between the eigenvalues and eigenfunctions for different values of the parameter $`\xi `$ is needed. Based on a recently developed Hamiltonian expansion for deformed billiards , it is easy to show that the eigenvalues and eigenfunctions for different values of the parameter $`\xi `$ can be obtained from the Hamiltonian $`_0+\mathrm{\Sigma }(\xi )`$ which is expressed in the basis of eigenstates at $`\xi =0`$ (from now on we will call $`\varphi _\mu `$ to these eigenstates),

$$\mathrm{\Sigma }_{\mu \nu }=\xi \times \mathrm{\Phi }_{\mu \nu }_𝒞g(q)\frac{\varphi _\mu }{𝐧}\frac{\varphi _\nu }{𝐧}dq.$$
(3)
The function $`g(q)`$ measures the strength of the change in the boundary condition along the contour. Within a perturbation theory it would represent the direction and strength of a distortion of the stadium . Here we use

$`g(q)=\{\begin{array}{ccc}\alpha & & 0q1,\\ (1+\alpha )\mathrm{sin}(q1)+\alpha & & 1<q1+\pi /2\end{array}`$

with $`\alpha =1/(2+\pi /2)`$ that could be assimilated to a dilation along the horizontal axis and a contraction along the perpendicular one. Notice that the integral above could be viewed as an inner product among the wave functions $`\frac{\varphi _\mu }{𝐧}`$ defined over $`𝒞`$. This relation defines an effective Hilbert space in a window $`\mathrm{\Delta }k`$ Perimeter/Area . The cut-off function $`\mathrm{\Phi }_{\mu \nu }=\mathrm{exp}\left[2(k_\mu ^2k_\nu ^2)^2/(k_0\mathrm{\Delta }k)^2\right]`$ restricts the effect of the perturbation to states in this energy shell of width $`Bk_0\mathrm{\Delta }k`$. It allows to deal with a basis of finite dimension with wave numbers around the mean value $`k_0`$ and restricting to a particular region $`\mathrm{\Delta }k`$ of interest.

Figure 1 shows the dependence of the energy levels on the perturbation. They exhibit many avoided crossings as $`\xi `$ is varied. While the energy levels show the typical behavior of a general system without constants of motion, we also recognize that some small avoided crossings are situated along parallel tilted lines. These energies correspond to the well known “bouncing ball” states which are highly localized in momentum. The selected perturbation does not modify substantially those states.

Since the LE is a classically motivated quantity, a Gaussian wave packet (with a mean value of momentum $`k_0`$ and velocity $`v_0`$) is a proper semiclassical selection for an initial condition. By evaluating its evolution in a system without perturbation ($`\xi =0`$) and other with perturbation strength $`\xi `$, we compute the LE (Eq. 1) as a function of time. At this point one must recognize that the choice of a semiclassical initial condition is very relevant in order to observe the ’Lyapunov’ regime .

While a global exponential decay of $`M(t)`$ can be clearly identified in almost any individual initial condition, the fluctuations for a system with $`k_0`$ not too large can introduce error in the estimation of the rate. Hence, we have taken an average over $`30`$ initial states. Fig. (2) (a) and (b) show typical sets of curves of $`M(t)`$ for $`k_0=50`$ and $`k_0=100`$ respectively. It can be seen that after a transient, $`M(t)`$ decays exponentially, $`\mathrm{exp}[t/\tau _\varphi ]`$. For $`\xi >\xi _c4.5/k`$ the decay rate $`\tau _\varphi `$ becomes independent of the perturbation and $`1/\tau _\varphi \lambda `$ with $`\lambda `$ the Lyapunov exponent of the classical system  in accordance with the conjecture. On the other hand, for large times $`M(t)`$ saturates to a finite value $`M_{\mathrm{}}1/N`$ with $`N`$ the effective dimension of the Hilbert space .

According to Ref.  the chaos controlled decay appears provided that $`\lambda >1/\stackrel{~}{\tau }`$ where $`\stackrel{~}{\mathrm{}}=v_o\stackrel{~}{\tau }`$ is the length over which the perturbation changes the quantum phase (mean free path) which, for a plane wave with wave number $`k`$ and velocity $`v_o`$. For a quenched disorder perturbation is evaluated from the FGR

$$\frac{1}{\stackrel{~}{\tau }_k}=\frac{2\pi }{\mathrm{}}\underset{\eta 0^+}{lim}\underset{k^{}}{}\left|\mathrm{\Sigma }_{k^{}k}\right|^2\frac{1}{\pi }\frac{\eta /2}{(E_k^{}E_k)^2+(\eta /2)^2}.$$
(4)
Ref.  realized that in the opposite regime of $`\lambda <1/\stackrel{~}{\tau },`$ the LE of an eigenstates $`\varphi _\mu `$ of $`_0`$ is just a survival probability an must decay exponentially under the action of the perturbation,

$$|\varphi _\mu \left|\mathrm{exp}[\mathrm{i}\left(_0+\mathrm{\Sigma }\right)t/\mathrm{}]\right|\varphi _\mu |^2\mathrm{exp}[t/\stackrel{~}{\tau }_\mu ],$$
(5)
given by Eq. 4 for$`\mathrm{}/B<t<\mathrm{}/\mathrm{\Delta }`$ ($`\mathrm{\Delta }`$ the mean level spacing) . The appearance of this FGR behavior requires that a typical matrix element $`U\left|\mathrm{\Sigma }_{\mu \nu }\left(\xi \right)\right|_{\mathrm{typ}.}`$of the perturbation to be $`U>\mathrm{\Delta }.`$ The Fourier transform of Eq. 5 is the LDOS which, although being discrete, would present a Lorentzian envelope  of width $`\mathrm{\Gamma }=1/\stackrel{~}{\tau }`$. In Ref.  it is conjectured that this decay can determine the LE decay with more general initial states. This is the regime controlled by the non-diagonal terms in the semiclassical expansion . Once the non-diagonal terms have decayed, one expects the chaos controlled decay of the diagonal ones will survive. This gives a cross-over criterion for the decay rate of the LE of $`1/\tau _\varphi =\mathrm{min}[\mathrm{\Gamma },\lambda ]`$ as the perturbative parameter $`\xi `$ changes.

The LDOS is shown in Fig. 3 for three different perturbation strengths. In contrast to the case of Ref.  our distribution is not Lorentzian. This is related to the fact that the used perturbation (the function $`g(q)`$) does not connect all different regions of phase space; for instance, the bouncing ball states are practically undisturbed by $`\mathrm{\Sigma }`$ determining the non-generic nature of the perturbation. In particular, we have evaluated the width $`\mathrm{\Gamma },`$ showing the spreading of the unperturbed eigenstates when expressed in terms of the new ones. The results show a linear dependence of $`\mathrm{\Gamma }`$ on $`\xi `$ shown in Fig. 3; that is, we obtain $`\mathrm{\Gamma }0.36\xi k^2`$. Moreover, taking into account that $`\lambda 0.86k`$, the critical value $`\xi _c`$ for the crossover from the $`\mathrm{\Gamma }`$ regime to the Lyapunov one is expected at $`\xi _c=2.4/k`$ (remember that from Fig. 2 it results $`\xi _c4.5/k`$). Then, for our system, the criterium works with a $`\mathrm{\Gamma }`$ given by the half width of the LDOS. This is shown in Fig. 4 where for perturbation strengths $`\xi <4.5k^1`$, the LE. decays as $`M(t)=\mathrm{exp}\left[t/\tau _\varphi \right]\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1}/\tau _\varphi =\mathrm{\Gamma }/2`$ for $`\lambda >\mathrm{\Gamma }/2`$.

These results contrast with the FGR dependence of $`1/\tau _\varphi \xi ^2`$ observed for weak perturbations. These are the Lorentz gas with a perturbed effective mass , the kicked top perturbed by a perpendicular delayed kick and general chaotic system perturbed by a quenched disorder where random matrix theory describes the $`\mathrm{\Gamma }`$ decay. In this context, the linear dependence of $`1/\stackrel{~}{\tau }`$ on $`\xi `$ may be considered as a further indicative that the physics of the LE decay can be very different from that described by Eq. 5 and that the result of Ref.  has more general validity than expected. In the non-perturbative regime, before the Lyapunov exponent takes over,the LE decays exponentially with a rate given by the perturbation dependent width of the LDOS. Another important feature is that $`\xi _c4.5/k0`$ when $`k\mathrm{}.`$ This confirms that in the classical limit Eq. (1) would decay with the Lyapunov regime regardless of the magnitude of $`\mathrm{\Sigma }`$, recovering the chaotic hypersensitivity to perturbations.

In summary, by studying one of the most important models in quantum chaos, a fully chaotic billiard system, we have shown that, for a wide range of parameters, the Loschmidt echo decays exponentially with rate given by the Lyapunov exponent of the classical system. Moreover, we have discussed the onset of this Lyapunov regime requires that $`\lambda >\mathrm{\Gamma }/2.`$ In the opposite situation, the presence of an exponential controlled by $`\mathrm{\Gamma }`$ even in absence of a generic perturbation described by the FGR, demands further studies to fully interpret the detailed mechanism controlling this regime. We finally remark that the $`M(t)`$ would behave much differently for intrisically quantum initial conditions. For an eigenstate of $``$ one finds a decay described by a FGR and it does not show a crossover into the Lyapunov decay . In the other quantum extreme, an initial state generated from the long time evolution of a semiclassical wave packet, we find a perturbation dependent Gaussian decay . These issues have begun to receive much attention  due to its strong connection with quantum computing stability, decoherence in waves, and quantum-classical transition. Furthermore, the dephasing time observed in transport experiments in mesoscopic devices shows a perturbation independent rate . So far, there is no consensus about the physical phenomenon causing it. Since the time scale $`\tau _\varphi `$ measured by the LE is a decoherence time and our methodology can obviously be adapted to treat the transport problem , our results open a rich field for exploration: the connection of both time scales.

We thank D. Cohen, R. Jalabert and M. Saraceno for very useful discussions and the support from SeCyT-UNC, CONICET, ANPCyT, ECOS-SeTCIP and Antorchas-Vitae. DAW received support from CONICET (Argentina) and AECI (Spain).


text:
# Discovery of the “missing” mode in HR 1217 by the Whole Earth Telescope

## 1 Introduction

After decades of trying, the search for solar-type oscillations in stars finally appears to have been successful (see, e.g., Bouchy & Carrier 2001 and Carrier et al. 2001). Although this led Gough (2001) to announce the “birth of asteroseismology”, for the past two decades asteroseismology has successfully investigated the interiors of many types of stars other than the solar-type stars. Remarkable success stories of observational and theoretical investigations of white dwarf stars and rapidly oscillating Ap stars have amply demonstrated the power of asteroseismology as a tool to advance our knowledge of the physics of stellar interiors and the details of stellar evolution (see, e.g., Kurtz et al 1989; Winget et al. 1991; Kawaler & Bradley 1994; Matthews et al. 1999).

With these successes, some mysteries have remained. In this Letter, we address apparently contradictory interpretations of the pulsation spectrum of the rapidly oscillating Ap star HR 1217. This star, discovered to be a pulsator by Kurtz (1982), was investigated with an extensive global campaign in 1986 (Kurtz et al. 1989). A key result from that data set was a list of six principal pulsation frequencies, reproduced in Table 1. As expected from the asymptotic theory of nonradial pulsations, five of the modes are nearly equally spaced in frequency.

The asymptotic frequency spacing, $`\nu _o`$, is a measure of the sound crossing time of the star, which in turn is determined by the star’s mean density and radius. With a typical mass of Ap stars of about 2M, $`\nu _o`$ reflects the radius of the star, with $`\nu _o`$ scaling as $`R^{3/2}`$. In the asymptotic limit, the number of nodes in the radial direction, $`n`$, is larger than the spherical degree $`\mathrm{}`$. Assuming adiabatic pulsations in spherically symmetric stars the pulsation frequencies are, to first order,

$$\nu _{n,\mathrm{}}=\nu _o(n+\mathrm{}/2+ϵ),$$

where $`ϵ`$ is a (small) constant (Tassoul 1980, 1990). Without precise identification of the degree ($`\mathrm{}`$) of the pulsation modes, asymptotic theory allows the frequency spacing to be uncertain by a factor of two, depending on whether modes of alternating even and odd $`\mathrm{}`$ are present (producing modes separated by $`\nu _o/2`$ in frequency), or only modes of consecutive $`n`$.

The results of the 1986 campaign were inconclusive as to whether $`\nu _o`$ was 68 $`\mu `$Hz or 34 $`\mu `$Hz. The prinicipal frequencies seen in the data are given in Table 1; they correspond to those found by Kurtz et al. (1989) for the 15 day stretch of best coverage (for comparison with the new data presented in Table 3). The highest frequency of HR 1217 in those data was 50 $`\mu `$Hz higher than the fifth mode, suggesting that $`\nu _o`$ was 34 $`\mu `$Hz. But the fine structure of the spacings was suggestive of alternating $`\mathrm{}`$ values. Fortunately, the two possible values could be assessed if the luminosity of the star were precisely known. If $`\nu _o`$ were 34 $`\mu `$Hz, then the radius of HR 1217 would be large enough that it would be far removed from the main sequence (i.e. more evolved) and therefore more luminous (Heller & Kawaler 1988). Matthews et al. (1999) used the hipparcos parallax measurement to place HR 1217 unambiguously close to the Main Sequence – meaning that $`\nu _o`$ is indeed 68 $`\mu `$Hz. This deepened the “mystery of the sixth frequency”, now $`\frac{3}{4}\nu _0`$ higher. No clear theoretical construct could explain it.

The asymptotic frequency spacing given in the equation above is valid only for linear adiabatic pulsations in spherically symmetric stars. However, the magnetic field, the chemical inhomogeneities, and rotation all contribute to break the spherical symmetry in roAp stars. Therefore, it is important to know the effects that these deviations from spherical symmetry have on the theoretical amplitude spectra of roAp stars, before comparing the latter with the observed amplitude spectra. The effects of the chemical inhomogeneities have been discussed recently by Balmforth et al. (2001), but those will not concern us further here. The effects of the magnetic field on the oscillations of roAp stars (Dziembowski & Goode 1996; Bigot et al. 2000; Cunha & Gough 2000), as well as the conjoined effect of rotation and magnetic field (Bigot 2002), have been determined by means of a singular perturbation approach. While generally the magnetic field effect on the oscillations is expected to be small, Cunha & Gough (2000) found that, at the frequencies of maximal magnetoacoustic coupling, the latter is expected to become significantly large, resulting in an abrupt drop of the separation between mode frequencies.

The observational consequence of the results of Cunha & Gough (2000) suggests that we should see equally spaced modes in roAp stars, with an occasional mode much closer to its lower frequency counterparts. More recently, Cunha (2001) suggested that the explanation of the strange separation between the last two modes observed in HR 1217 could rest on the occasional abrupt decrease of the large separations predicted by Cunha & Gough (2000). For this prediction to hold, she argued that the observations of Kurtz et al. (1989) must have missed detecting a mode at a frequency 34$`\mu `$Hz higher than that of the fifth mode they observed. She predicted that new, more precise measurements would find this “missing mode” if the Alfvénic losses were not large enough to stabilise it. Detailed re–examination of the data from 1986 shows no peak at the key position approximately 33 $`\mu `$Hz above f5 at the 0.1 mma level.

In 2000 November, we began an extensive, coordinated global photometry campaign on HR 1217 using the Whole Earth Telescope. A complete analysis of this extensive data set, which addresses many other aspects of roAp stars, is in preparation. In this Letter, we present a preliminary analysis of data from that run that clearly shows a previously unseen pulsation mode at a frequency about 36 $`\mu `$Hz above the fifth frequency, as predicted by Cunha (2001). In the next section, we describe the observational procedures and the data coverage and reduction. Section 3 presents the preliminary frequency analysis, and the results are discussed in Section 4, along with a brief discussion of the impact of this result on the theory of pulsations of roAp stars.

## 2 Observations

The WET run on HR 1217 began on 2000 November 6 at selected sites, and continued through early 2000 December. The bulk of the data, with the best global coverage, were obtained during 2000 November 14-30. A complete analysis of all of the available data is currently underway. For this Letter, we concentrate on the central portion of the WET run, with data from five sites. This subset of the full data set provides high signal–to–noise and a reasonable global coverage. It also extends over slightly more than one rotation cycle of HR 1217. Since the pulsation amplitude is modulated with the rotation period, this data subset is just long enough to begin to resolve rotational sidelobes of the main peaks.

Table 2 lists the individual observing runs in this data set. The telescopes used range in aperture from 0.6 m to 2.1 m. Data from all sites were obtained using photoelectric photometers, with 10 s individual integrations. At Beijing Astronomical Observatory, McDonald Observatory, Mauna Kea Observatory, and Observatorio del Teide, the observers used three-channel photometers that are functionally similar to the equipment described in Kleinman et al. (1996). The South African Astronomical Observatory observations were made with a single-channel photometer, and the observations at CTIO with a two–channel photometer. At all sites, observations were made through a Johnson $`B`$ filter, along with neutral density filters when needed to keep the count rates below $`10^6`$ s-1. Following the procedures described in Kleinman et al. (1996), the sky background was continuously monitored with the three-channel instruments. At sites using two-channel and single channel photometers, the sky was obtained several times during the night at irregular intervals, and then interpolated during reduction.

As can be seen in Table 2, we obtained 215.2 hr of observations during the interval from 2000 November 14–30, resulting in a duty cycle of 53%. Longitude coverage was adequate, though the longitudes around central Asia were not as well covered as the others.

## 3 Frequency analysis

The Fourier transform (FT) of the reduced data, in the frequency range where the pulsations are significant, is shown in in the top panel of Figure 1. The spectral window, shown in the middle panel of the figure, shows the response of the FT to a single, noise-free sinusoid sampled at the same times as the light curve of HR 1217. The side peaks correspond to aliases of 1 d-1 and 2 d-1. They are approximately 40% of the amplitude of the principal peak, and are caused by the (small) daily gaps present in the data causing incomplete global coverage.

The principal periodicities that we found in HR 1217 are listed in Table 3. Following initial identification of the main peaks in the FT, we did a successive least-squares fit to the light curve including all of the main peaks. We then included the rotational sidelobes in the fit, sequentially. Throughout this process we prewhitened the data by removing noise-free sinusoids at the fitted frequencies, amplitudes, and phases. We stopped when none of the remaining peaks was above the noise level. In all, we found 14 significant periodicities in this data set. In addition to the 7 principal frequencies, both rotational sidelobes of f3 and f4 were found. We also found the low-frequency rotational sidelobe of f2, f5, and f7. The frequencies listed in Table 3 are from the fit that included all 14 frequencies.

The bottom panel of Fig. 1 shows the FT of the residual light curve following the removal of 14 frequencies, on the same scale as the top panel. There are some residual peaks in this plot at interesting frequencies. Analysis of the full data set, including runs outside of the subset that we used, shows that some of these are real. They will be described in further detail in the full analysis of the data which is in preparation.

## 4 Results

### 4.1 Comparison with the 1989 data

Early in the run, it became clear that HR 1217 was pulsating with the same frequencies that were present in the 1986 data analysed by Kurtz et al. (1989). Tables 1 and 3 show that the principal frequencies from the 1986 study (f1 through f5 and f7) are highly consistent over a time span of 15 yr. Some of the amplitudes of these modes are higher in 2000 and some lower by small amounts than they were in 1986, but it is the frequencies (and presence or absence) of the modes that are of interest here.

The chief difference between the 2000 data and 1986 data is the presence of a frequency at 2791 $`\mu `$Hz listed as f6 in Table 3. That mode was not detected in the data of Kurtz et al. (1989 - Table 1) but was a clear signal in the WET Xcov20 2000 data. To ensure that this frequency is not an artefact of the data reduction algorithm, we repeated the frequency analysis of our data fitting just the large-amplitude peaks f2, f3, f4, and f5, and their rotational sidelobes (if present). We then removed those 10 frequencies. The results are illustrated in Fig. 2. This figure shows the original FT, and the FT of the data simulated by including f2-f5 and their rotational sidelobes. Clearly, there is excess signal at the positions of f1 and f7, but also at 2791 $`\mu `$Hz as well.

Thus we conclude that the “new” frequency, f6, is real. Table 3 shows that it lies at nearly $`\nu _o/2`$ above f5, as expected if it is a normal $`p`$mode and $`\nu _o68\mu `$Hz. It is much closer to f7 than $`\nu _o/2`$, as predicted by Cunha (2001).

### 4.2 Implications for roAp stars

Cunha (2001) speculated that the position of the f7 peak in the Kurtz et al. (1989) data is consistent with her model of the normal mode structure in Ap stars when magnetic fields are important to the pulsation dynamics. Since that peak was $`\frac{3}{4}\nu _0`$ above f5 (which is inexplicable in asymptotic theory), she suggested that there should be a peak at $`\nu _o/2`$ above f5. That is precisely what we see in the data from the WET run in November 2000 with the discovery of f6.

As discussed earlier, other explanations for the frequency spacing pattern from f5 to f7 are in direct conflict with the now well-determined hipparcos luminosity of HR 1217. We therefore conclude that the frequency pattern in HR 1217 suggests that the pulsations we see in this star are consistent with normal $`p`$mode pulsations whose frequencies are, in some cases, strongly affected by the magnetic field of the star.

Cunha (2001) suggested that large Alfvénic losses could help explain the missing f6 in the 1986 data, as these losses are maximal at the frequencies where the large separations experience the abrupt decrease. This energy loss could either stabilize the mode or contribute to decrease its amplitude (although it is not clear how the growth rates relate to the amplitude of the modes in roAp stars).

Since f6 is observed in the present data, the possibility that the Alfvénic losses are large enough to stabilise this mode can be ruled out, at least at the time of these observations. Whether at the time of the previous observations the efficacy of the magnetoacoustic coupling (which depends, among other things, on the exact frequency of the mode and on the characteristics of the magnetic field) was different, is something to which we do not have an answer. An attempt to monitor the amplitude of f6, as well as that of the other modes, in the future might, therefore, be worthwhile. However, the magnetic field does produce an important observable effect on the frequency of f7.

With the detection of f6, we move closer to a detailed understanding of the pulsation mechanism in roAp stars. Most intriguingly, this result for HR1217 suggests that, with appropriately detailed models, we may soon be able to probe the magnetic field structure below the surfaces of these stars through their pulsation frequencies – another application of asteroseismology to probing stellar interiors.

## Acknowledgments

We gratefully acknowledge support from the U.S. National Science Foundation through grant AST-9876655 to Iowa State University, and funding by UNESCO through the International Institute of Theoretical and Applied Physics at Iowa State. M.C is supported by FCT-Portugal through the grant PD/18893/98 and the grant POCTI/1999/FIS/34549 approved by FCT and POCTI, with funds from the European Community programme FEDER. PM is supported by KBN (Poland) through grant 5-P03D-012-20.