text:
# The case of the missing neutrality

## 1 Introduction

The influence of neutral networks in evolutionary processes was first elucidated by Peter Schuster’s group in Vienna in 1996. Put simply, two genotypes are considered neutrally equivalent if they map to the same phenotype. A neutral network is a set of genotypes connected by this neutrality relationship on links with Hamming distance 1 (ie each link of the network corresponds to a mutation at a single site of the genome).

These researchers noted that evolution tended to proceed by diffusion along these neutral networks, punctuated occasionally by rapid changes to phenotypes as an adaptive feature is discovered. The similarity of these dynamics with the theory of Punctuated Equilibria was noted by Barnett. It was also noted that if a giant network existed that came within a hop or two of every possible genotype, evolution will be particularly efficient at discovering solutions.

Most work on neutrality in evolution uses the genotype-phenotype mapping defined by folding of RNA. This mapping is implemented in the open source Vienna RNA package111http://www.tbi.univie.ac.at/~ivo/RNA, so is a convenient and well known testbed for ideas of neutrality in evolution.

Also in 1996, I developed a definition of the genotype-phenotype mapping for Tierra, which was first published in 1997. I noticed the strong presence of neutrality in this mapping at that time, which was later exploited to develop a measure of complexity of the Tierran organism. In 2002, I started a programme to visualise Tierra’s phylogenetic trees and neutral networks in order to “discover the unexpected”. Two key findings came out of this: the first being that Tierra’s genebanker data did not provide clean phylogenetic trees, but had loops, and consisted of many discontinuous pieces. This later turned out to be due to Tierra’s habit of reusing genotype labels if those genotypes were not saved in the genebanker database. This might happen if the population count of that genotype failed to cross a threshold. This is all very well, except that a reference to that genotype exists in the parent field of successor genotypes. The second big surprise was the paucity of neutral mutations in the phylogenetic tree. We expect most mutations to an organism to be deleterious, and so expect that neutral mutations will have disproportionately greater representation amongst successor genotypes than one would expect if each mutation was equally likely.

## 2 Neutrality in Tierra

Tierra is a well known artificial life system in which small self-replicating computer programs are executed in specially constructed simulator. These computer programs (called digital organisms, or sometimes “critters”) undergo mutation, and radically novel behaviour is discovered, such as parasitism and hyperparasitism.

It is clear what the genotype is in Tierra, it is just the listing of the program code of the organism. The phenotype is a more diffuse thing, however. It is the resultant effect of running the computer program, in all possible environments. Christoph Adami defined this notion of phenotype for a similar artificial life system called Avida. In Avida, things are particularly simple, in that organisms either reproduce themselves at a fixed replication rate, or don’t as the case may be, and optionally perform range of arithmetic operations on special registers (defined by the experimenter).

In Tierra, organisms do interact with each other via a template matching mechanism. For example, with a branching instruction like `jmpo`, if there is a sequence of `nop0` and `nop1` instructions (which are nooperations) following the branch, this sequence of 1s and 0s is used as a template for determining where to branch to. In this case the CPU will search outwards through memory for a complementary sequence of `nop0`s and `nop1`s. If the nearest complementary sequence happens to lie in the code of a different organism, the organisms interact.

To precisely determine the phenotype of a Tierran organism, one would need to execute the soup containing the organism and all possible combinations of other genotypes. Whilst this is a finite task, it is clearly astronomically difficult. One means of approximation is to consider just interaction of pairs of genotypes (called a tournament). Most Tierran organisms interact pairwise — very few triple or higher order interactions exist. Similarly, rather than running tournaments with all possible genotypes, we can approximate matters by using the genotypes stored in a genebanker database after a Tierra run. In practice, it turns out that various measures, such as the number of neutral neighbours, or the total complexity of an organism are fairly robust with respect to the exact set of organism used for the tournaments.

So the procedure is to pit pairwise all organisms in the genebanker against themselves, and record the outcome in a table (there is a small number of possible outcomes, which is detailed in ). A row of this table is a phenotypic signature for the genotype labeling that row. We can then eliminate those genotypes with identical signatures in favour of one canonical genotype. This list of unique phenotypes can be used to define pragmatically a test for neutrality of two different organisms. Pit each organism against the list of unique phenotypes, and if the signatures match, we have neutrality. The source code for this experiment is available from the author’s website.222http://parallel.hpc.unsw.edu.au/getaegisdist.cgi/getsource/eco-tierra.3, version 3.D3

Tierra has three different modes of mutation:

A site of the soup is randomly chosen and mutated;

Data is mutated during the copy operation;

Instructions occasionally produce erroneous results

Furthermore, in the case of cosmic ray and copy mutations, a certain proportion of mutations involve bitflips, rather than opcodes being substituted uniformly. This proportion is set as a parameter in the soup\_in file (`MutBitProp`) — in these experiments, this parameter is set to zero.

In order to study the issue of whether neutrality is greater or less than expected in Tierra, I generated three datasets with each of the 3 modes of mutation operating in isolation. The sizes of each data set was 69,139, 87,003 and 198,982 genotypes respectively, generated over a time period of about 1000 million executed instructions. Genebanker’s threshold was set to zero, so all genotypes were captured. This led to a proper phylogenetic tree. After performing a neutrality analysis, a set of 83, 86 and 158 unique phenotypes was extracted as the test set for the tournaments.

Since the neighbourhood size increases exponentially with neighbourhood diameter, I restrict analysis to single site, or point mutations. In each data set, around 7% of these genotypes were created were created by a mutation at a single site and were neutrally equivalent to its parent. For each of these, I compute the number of neutral neighbours $`n_i`$ existing in the 1 hop neighbourhood of the parent genotype $`i`$, which is of size $`32^\mathrm{}_i`$, where $`\mathrm{}_i`$ is the length of the genome. For a given parent $`i`$, the ratio

$$r_i=\frac{\nu _i32^\mathrm{}_i}{o_in_i}$$
(1)
gives the proportion of neutral links actually followed relative to the number of neutral links available (neutrality excess), where $`\nu _i`$ is the number of neutrally equivalent offspring, and $`o_i`$ the total number of offspring and the size of the 1 hop neutral neighbourhood. Fig. 1 shows the running average of this quantity over these transitions, with the genotypes numbered in size order.

In this analysis, no selection is operating, so one would expect that the neutrality excess should be identical to 1. However, in the case of instruction flaws, it is rather unpredicatable what the effect is. In the case of cosmic ray mutations, 50% of time one would expect the parent to be mutated, rather than the daughter. In the case of a mutation affecting a crucial gene of a parent genotype, the organism may not be able to reproduce at all, thus favouring neutral mutations. Only copy mutations should affect all sites of the genome equally, leading to a neutrality excess equal to one. The measured value, however is about 1.3, substantially greater than one. The reason for this is not known at this point in time.

The the datasets were further subsetted to include just those transitions whose daughter organism successfully reproduced, ie with a maximum population count greater than 1. The neutrality excess in this case is substantially less than 1, so something in Tierran evolution is favouring nonneutral evolution.

## 3 Vienna RNA folding experiments

It is quite well known that evolution using the RNA folding map exhibits a great deal of neutrality, at least for a standard genetic algorithm optimising a well defined fitness function. Tierra is a coevolutionary system, and does not have a well defined fitness function. Rather, the chance of an organism surviving at each time step depends on what other organisms are in the environment at the time, and has a significant component of contingency.

One possible cause for the repression of neutrality in Tierra is this coevolutionary nature. The reasoning follows from the idea of the Red Queen effect. This says that organisms must continuously evolve just to remain adaptive.333Like the Red Queen in Lewis Caroll’s Through the Looking Glass, who had to keep running, just to stay where she was. In such a circumstance, neutral evolution is maladaptive, and surely be suppressed.

A convincing argument in favour of this hypothesis would be the demonstration of a coevolutionary system based on the RNA folding map that exhibited this repression of neutrality. Whilst I haven’t achieved this goal, I will report on a couple of attempts.

The first attempt is an instantiation of the simplest possible coevolutionary system. It consists of two populations: a tracker population $`T`$ which attempts to be as similar to the other population as possible, and an evader population $`E`$ that attempts to be as different from the tracker population as possible. The Vienna RNA folding library is used, and fitness functions $`f_T(x,E)`$, $`f_E(x,T)`$ defined for the trackers and evaders respectively, based on the average distance between their phenotypes:

$`f_T(x,E)`$ $`=`$ $`{\displaystyle \frac{1}{N_E}}{\displaystyle \underset{yE}{}}d(x,y),xT`$
$`f_E(x,T)`$ $`=`$ $`{\displaystyle \frac{1}{N_T}}{\displaystyle \underset{yT}{}}d(x,y),xE`$ (2)
where $`N_T`$ and $`N_E`$ are the population counts of trackers and evaders respectively, and $`d(x,y)`$ is the string edit distance444Please consult the Vienna RNA package documentation for a precise definition of string edit distance between the folded structure of $`x`$ and $`y`$.

I implemented a simple genetic algorithm with just a point mutation operator. All RNA strings have the same length. During reproduction, each RNA string copies itself, possibly with a mutation. During the selection step, the least fit 50% of organisms are culled, bringing the population count back to the starting value. For the results presented in Fig. 2, the GA parameters are shown in Table 1. Source code for this experiment is available from the author’s website.555http://parallel.hpc.unsw.edu.au/getaegisdist.cgi/getsource/rnafold/, version D1

Fig. 2 shows the neutrality excess for trackers and evaders, defined in the same way as eq (1). A third baseline run of a single population tracking a fixed target is also shown (corresponding to a classic fixed fitness function genetic algorithm). The most obvious thing about these results is that neutrality is highly adaptive to evaders, who clearly are trying to make themselves occupy as small a footprint in phenotype space as possible. Given the evaders propensity to stay still, trackers will tend to behave like their fixed target counterparts.

This indicated that trackers were dominating over the evaders. Another experiment I performed was where all coevolving populations were symmetric. In this case, each population would track a second and evade a third population. The simplest arrangement of these has 3 coevolving populations in a “Rock, Scissors, Paper” configuration, so no one population dominates. I also tried other combinations up to 10 separately coevolving populations, wired up randomly to each other (but fixed at the start of the experiment). In all cases, the results were much the same — the neutrality excess was less than for the trackers in Fig. 2, but still just slightly greater than 1.

So how do can we introduce the Red Queen effect to this system? One possibility I haven’t explored as yet is to somehow give evaders room to move, perhaps by including a genome lengthening operator as part of the GA.

## 4 Conclusion

The suppression of neutrality in Tierran evolution is a real effect. It is quite likely that this is a Red Queen effect, with organisms needing to change to remain adaptive. Experiments with using the RNA folding map to try to reproduce this effect have proven inconclusive. However, it was noted that there is significant evolutionary pressure to increase neutrality in evading populations.

## Acknowledgments

I would like to thank the Australian Centre for Advanced Computing and Communications for a grant of computing time used in this project.


format: text

text:
# A note on transition, turbulent length scales and transport in differentially rotating flows

## 1 Introduction

If it is widely accepted that plane shear flows are subject to non-linear instability, the situation concerning their rotating counterpart is more controversial. Indeed, rotation alone is known to introduce constraints that could result in the stabilization of the flow. Differentially rotating flows are ubiquitous in astrophysics. From stellar interiors to accretion disks, from galactic rotation to planetary disks, they are met everywhere.While in many communities it is generally accepted that differential rotation could give rise to turbulence, the debate continues in the accretion disk community. Since the dawn of accretion disk modeling in the early 1970’s, and in the absence of quantified transport properties supported by an identified physical mechanism, turbulence transport has been modeled with ad hoc prescriptions for the anomalous viscosity (the most popular being the so-called ”alpha prescription”). Differential rotation has been proposed early as the possible source of this anomalous transport (Shakura & Sunyaev 1973). Recent numerical simulations fail to show evidence of such instability in accretion disks. However, whether or not these simulations are adequate or sufficient to give a definitive answer to this issue remains open to question (Longaretti 2002; Richard 2001, 2003). The difficulty of studying turbulence and instabilities has been reflected in the lack of a detailed physical mechanism along with dynamical and transport properties to apply in Astrophysical models. Several relevant instabilities are currently being investigated in the accretion disk community (Klahr & Bodenheimer (2003), Longaretti (2002), Balbus & Hawley (1991)), the most studied being the Magneto-Rotational Instability (MRI) (Balbus & Hawley 1991). This instability, linear in nature, is an excellent candidate as long as the disk is ionized enough, but is more unlikely to be efficient in cold disks such as proto-planetary nebulae (leading to the introduction of the ”dead zone” layered accretion disk model by Gammie (1996)). The issue of differential rotation has been recently addressed with rather different approaches in several publication such as Chagelishvili et al. (2003), Longaretti (2002), Longaretti (2003), Richard (2003) and Tevzadze et al. (2003). It is our purpose to try to present these different results in a synthetic manner (section 2) and to propose in their light a possible basic coherent scenario for the development of instabilities and the properties of the resulting turbulent state (section 3).

## 2 Flow stability

In this section we briefly comment on various published results concerning the stability of differentially rotating flows.

### 2.1 The effect of rotation

If plane shear flows are known to bear non-linear instabilities, the picture is quite different where the constraint of rotation is added. In a widely referred-to analysis, Balbus et al. (1996) argued that hydrodynamic turbulence would be difficult to sustain in Keplerian flows without taking into account hydro-magnetic effects. Based on the equations of evolution of the turbulent velocity fluctuations energy, they deduced that rotating flows could not experience the turbulence due to the energy sink introduced by the Coriolis acceleration. Taking a closer look, one should note that the Coriolis force embodying the rotation influence on the fluid motion does not have any effect on the total energy equation. Its effect is to redistribute energy between the radial and azimuthal components of the motion, vanishing when these equations are summed. This is the natural consequence of the fictitious character of the Coriolis force. Therefore, the rotation merely stiffens the system, add constraints, but does not introduce any energy sink. The rotation and curvature of the flow introduces fundamental differences in comparison with plane shear flow. It does not necessary mean that the instability mechanisms present in plane flows disappear, but more likely that the stability and the characteristics of the resulting turbulence are modified (One can refer for example to Longaretti (2003) for a comparative discussion on the role of rotation, curvature and viscous dissipation in plane and rotating shear flows).

### 2.2 Finite amplitude perturbations and transition to turbulence

In contrast with their linear counterpart, non-linear instabilities require the presence of finite (non-infinitesimal) amplitude perturbations to be triggered. Chagelishvili et al. (2003) (and sequel Tevzadze et al. (2003)) proposed the mechanism of transient growth and bypass mechanism (used for many years in the aerodynamics community) to provide the required amplitudes. Even in a linearly stable flow, the linear operator can provide transient amplification of perturbations. The results from Chagelishvili et al. (2003) show that the initial perturbation amplitude can be greatly amplified. This amplification can introduce into the flow perturbations of amplitudes sufficient so that their non-linear interactions can no longer be neglected.

Richard (2003) recently proposed a simple model describing the necessary conditions for self-sustained turbulence in differentially rotating flows, where it can be seen that the energy extraction is directly proportional to the shear present in the base flow (a classical result for shear flows). It also implies that there is a critical Rossby number ($`Ro=r_r\mathrm{\Omega }/2\mathrm{\Omega }`$, where $`r`$ is the local radius and $`\mathrm{\Omega }`$ is the mean rotation.) above which the energy extraction is sufficient to compensate for the stiffness introduced in the system by the mean rotation. It has to be pointed (as Chagelishvili et al. (2003) also suggested) that the non-linear interactions (also referred to as turbulent diffusion) do not participate in the energy extraction, but only redistribute it and counteract the effects of rotation. In this scenario, once the critical Rossby number has been reach (Richard 2003) and the critical amplitude is present, the flow can then undergo a transition from its laminar state to a state where non-linear shear turbulence is developed.

## 3 Turbulent flow properties

Richard (2003) focused on the stability properties of differentially rotating flows. Using a similar approach, we wish here to address the properties of the turbulent state of the bifurcated flow, in particular the turbulent transport and scale properties .

### 3.1 Rossby number and turbulent scaling

Rotation has a known tendency to constrain larger scales in a two-dimensional state while efficient turbulent diffusion is achieved by three-dimensional motions. Typical geophysical rotating flows exhibit both two-dimensional (at large scales) and three-dimensional (at smaller scales) structures. Within this picture, the relevant scale determining the turbulent diffusion should be the larger scale for which the turbulence can achieve to be three-dimensional.Baroud et al. (2003) performed laboratory experiments on a rotating annulus, and have shown that low Rossby numbers are associated with two-dimensional turbulence whereas higher Rossby numbers (of order unity in their experiments) are associated with three-dimensional turbulent structures. The Rossby number for the mean flow can be re-written as

$$Ro=\frac{1}{2}\frac{r}{(_r\mathrm{ln}\mathrm{\Omega })^1},$$
(1)
where it appears as the ratio of local radius over the characteristic length scale of the shear. It can also be linked to a quantity more often referred to in astrophysics, the epicyclic frequency $`\omega `$, through the relation

$$\frac{\omega ^2}{4\mathrm{\Omega }^2}=(1+Ro).$$
(2)
The Rossby number of a turbulent structure of characteristic length scale $`\lambda `$ and velocity $`u`$ (”turbulent Rossby number” hereafter), rotating with mean flow velocity, can be approximated by $`Ro_\lambda ru/2\lambda ^2\mathrm{\Omega }`$ : The denominator of the Rossby number is twice the rotation experienced by the turbulent structure, i.e. $`2(\mathrm{\Omega }\pm u/r)2\mathrm{\Omega }`$ ; The numerator is the derivative with respect to the radius of the turnover time, i.e. $`r_r\mathrm{\Omega }+r_r(u/\lambda )r_r\mathrm{\Omega }+ru/\lambda ^2ru/\lambda ^2`$, where we have used the relation

$$u/\lambda r_r\mathrm{\Omega },$$
(3)
from Richard (2003), and the condition $`\lambda <<r`$. We can then write, from Eq.(3) that

$$Ro_\lambda \frac{r}{\lambda }Ro.$$
(4)
Hence, the ratio of the characteristic turbulent length scale over the local radius writes as

$$\frac{\lambda }{r}\frac{Ro}{Ro_\lambda }.$$
(5)
From Eq. (4), it follows that (for a given Rossby number associated with the mean flow) the turbulent Rossby number increases when going to smaller scales along the turbulent cascade. Remembering that high (resp. low) Rossby number are associated with three- (resp. two-) dimensional motions, this result is coherent with the classic picture of a turbulent spectrum showing two-dimensional structures at large scales and three-dimensional ones at smaller scales.

### 3.2 Turbulent diffusion

From an Astrophysical point of view, the main interest in studying turbulence might well be to characterize the turbulent transport of angular momentum or mass, by deriving the relevant model for the turbulent viscosity, which is classically written as

$$\nu _tu\lambda _\nu ,$$
(6)
where $`\lambda _\nu `$ stand for the largest turbulent scale participating in the diffusion process and $`u`$ its associated velocity. From Eqs. (3) and (5), it derives that

$$\nu _t\left(\frac{Ro}{Ro_{\lambda _\nu }}\right)^2r^3_r\mathrm{\Omega }.$$
(7)
Two of the results from Richard & Zahn (1999), derived from experimental data from Wendt (1933) and Taylor (1936), are :

\- The ratio $`\lambda _\nu /r`$ is constant. From Eq.(5) it follows that $`(Ro/Ro_{\lambda _\nu })`$ is a constant.

\- The turbulent viscosity scales as $`\nu _t=\beta r^3_r\mathrm{\Omega }`$, where $`\beta `$ is a numerical constant. This result combined with Eq.(7), implies that $`\beta (Ro/Ro_{\lambda _\nu })^2`$.

Experiments from Taylor (1936) used by Richard & Zahn (1999) to derive the $`\beta `$ formulation of turbulent viscosity in rotating flows, all had the same maximum $`Ro`$. The reason is that the Rossby number within a laminar Couette-Taylor flow when the inner cylinder is at rest is given by $`Ro=R_i/r`$, where $`R_i`$ is the radius of the inner cylinder and r is the local radius. It follows that the absolute value of the Rossby number is maximum at the inner boundary of the flow and decreasing outward. In his classic experiments, Taylor modified the aspect ratio of his apparatus by changing the radius of the outer cylinder, hence did not modify the maximum Rossby number value in the flow. It follows that in Eq. (7), $`Ro`$ is to be treated as a constant for this set of experiments. We conclude that both the apparent scaling of $`\lambda _\nu `$ and the constant value of the $`\beta `$ parameter derived in Richard & Zahn (1999), can be explained if we postulate that there exists a critical (transitional) $`Ro_{\lambda _\nu }`$ in the turbulent cascade above which the turbulence is three-dimensional and under which it is two-dimensional, and that its value is a constant that does not depend on the flow itself. Concluding on the form of the turbulent viscosity, in the general case (from Eq.(7)) it writes

$$\nu _t\frac{1}{Ro_{\lambda _\nu }}\frac{(r^2_r\mathrm{\Omega })^3}{4r\mathrm{\Omega }^2}.$$
(8)
In the particular case of Keplerian flows, it reduces to

$$\nu _t\mathrm{\Omega }_Kr^2,$$
(9)
where $`\mathrm{\Omega }_K=\sqrt{GM/r^3}`$ is the Keplerian angular velocity. This form is similar to the $`\beta `$-viscosity prescription, first introduced by Lynden-Bell & Pringle (1974) and more recently reintroduced by Richard & Zahn (1999) (based on laboratory fluid experiment data) and Duschl et al. (2000), and applied by Huré et al. (2001), Davis (2003), Wehrstedt & Gail (2003) and Lachaume et al. (2003).

## 4 Discussion

Based on recent developments, we have proposed a simple scheme for the transition and self-sustenance of turbulence in non-magnetic differentially rotating flows. In this scenario, initial velocity fluctuations are amplified by their interaction with the mean shear (provided by the linear operator in the equations of motion). While this coupling between perturbations and averaged flow can not by itself destabilize the flow (due to its linearly stable nature) it can provide transient growth introducing finite amplitude fluctuations. These fluctuations in turn feed the non-linear interactions and trigger the instability, under the condition that the Rossby number of the mean flow has reached a critical value sufficient for the non-linear energy transfer to overcome the stiffness introduced by the mean rotation. We conjecture that once the flow has become unstable, two- and three-dimensional turbulence coexist at different scales. The bi- or three-dimensional character of the motions is dictated by the value of the Rossby number at a given scale. Finally we have introduced a new formulation for the turbulent viscosity of such turbulence, which reduces to the $`\beta `$ prescription in the case of constant Rossby number flows, including Keplerian rotation.

###### Acknowledgements.

The authors wish to thank Dr. Jeffrey Cuzzi and Dr. Robert Hogan for their comments and support. D.Richard is supported by a Research Associateship from the National Research Council / National Academy of Sciences. This research has made use of NASA’s Astrophysics Data System.


text:
# Combined BIMA and OVRO observations of comet C/1999 S4 (LINEAR)

## 1 Introduction

The volatile component of comets provides an excellent tracer of the chemical composition present during the formation of the Solar System (Mumma et al., 1993; Irvine & Bergin, 2000; Irvine et al., 2000). Recent years have seen remarkable progress in this area with the appearances of two especially active comets, Hyakutake (C/1996 B2) and Hale-Bopp (C/1995 O1), that passed close to Earth. Observations at (sub) millimeter wavelengths of these spectacular comets, in addition to numerous less prominent but no less interesting comets over the last decade, indicate that the composition of comets is heterogeneous. This holds not only from comet to comet (e.g., A’Hearn et al. 1995), but also within the same comet as was shown with millimeter-interferometers that resolve the emission from the coma (e.g., Blake et al. 1999). The latter work illustrated the information that high spatial resolution observations can offer about cometary composition and evaporation processes.

In the Summer of 2000, the comet C/1999 S4 (comet LINEAR from here on) passed close to the Sun (0.765 AU) and at 0.374 AU to the Earth (Durig et al., 1999). Its orbit identified it as a dynamically young comet, and it showed a high level of activity at large heliocentric distance. As is not uncommon in dynamically young comets, the increase in its activity lagged significantly behind what was expected after its discovery, reaching only a medium level of activity at perihelion passage. It generated some excitement when, at closest approach to the Sun, the nucleus of the comet disintegrated and the residual material faded away in a few days, first reported by Kidger (2000). The comet’s characteristics and its breakup are discussed in detail by (Bockelée-Morvan et al. 2001; Farnham et al. 2001; Lisse et al. 2001; Mäkinen et al. 2001; Mumma et al. 2001; Weaver et al. 2001; Bonev et al. 2002; Altenhoff et al. 2002; Tozzi & Licandro 2002; Schulz & Stüwe 2002; Cochran & Cochran 2001; Jockers et al. 2001; Rodríguez et al. 2002; Kidger 2002, and Magee-Sauer et al. 2001). Of most relevance to our discussion here are: the conclusion that the comet’s nucleus may have broken up already several weeks before it fully disintegrated in late July (July 23–24); the molecular-line observations reported by Bockelée-Morvan et al. including the HCN-production rate; and the report by Magee-Sauer et al. (2001) of infrared HCN line observations.

This Paper describes the results of an observing campaign of the volatile composition of comet LINEAR, carried out jointly by the millimeter arrays at Hat Creek (operated by the Berkeley-Illinois-Maryland Association) and Owens Valley (operated by the California Institute of Technology), both in California (§2). We present a weak detection of HCN $`J`$=1–0 in auto-correlation mode and upper limits on the emission in cross-correlation mode, and place limits on the spatial extent of the emission (§3). Using models of HCN excitation in comets from Crovisier (1987) (see also Bockelée-Morvan et al. 1984) we investigate whether the spatial extent matches expectations (§4). Our conclusions are summarized in §5. The Appendix discusses the methods used to combine data from the two millimeter arrays.

## 2 Observations

Observations of comet LINEAR were obtained between 2000 July 17 and 2000 July 27 at the millimeter interferometers of the Berkeley-Illinois-Maryland Association (BIMA)111The BIMA Array is operated by the Berkeley-Illinois-Maryland Association under funding from the National Science Foundation. at Hat Creek, California ($`121^{}\mathrm{\hspace{0.17em}28}^{}\mathrm{\hspace{0.17em}18}\stackrel{}{\mathrm{.}}49`$ West, $`40^{}\mathrm{\hspace{0.17em}49}^{}2\stackrel{}{\mathrm{.}}5`$ North; altitude 1043 m), and the Caltech Millimeter Array at the Owens Valley Radio Observatory (OVRO)222The Owens Valley Millimeter Array is operated by the California Institute of Technology under funding from the U.S. National Science Foundation (#AST96–13717). near Big Pine, California ($`118^{}\mathrm{\hspace{0.17em}16}^{}\mathrm{\hspace{0.17em}55}\stackrel{}{\mathrm{.}}92`$ West, $`37^{}\mathrm{\hspace{0.17em}14}^{}\mathrm{\hspace{0.17em}2}\stackrel{}{\mathrm{.}}04`$ North; altitude 1222 m). The six 10-meter antennas of the OVRO array were in ‘compact’ configuration, the ten 6-meter antennas of the BIMA array in ‘C’ configuration: both configurations result in comparable projected baselines: 2–20 k$`\lambda `$ at BIMA and 3–15 k$`\lambda `$ at OVRO. Since these were summer-time observations, when atmospheric stability can be poor for millimeter-interferometry, these compact configurations were adopted to avoid decorrelation by atmospheric phase fluctuations on long baselines.

OVRO observed the comet in crosscorrelation mode from rise in the early morning local time until set in the late afternoon local time. BIMA operated in crosscorrelation mode from rise until mid-morning local time. After mid-morning the phase stability at Hat Creek deteriorated significantly during the observing period, and the observing mode was switched to autocorrelation. This observing mode1 does not suffer from phase-decorrelation but has a large beam-dilution in the $`143^{\prime \prime }`$ full-width at half-maximum (FWHM) primary beam of the six-meter antennas. The phase stability at Owens Valley, approximately 500 km south-east of Hat Creek and on the opposite side of the mountain range of the Sierra Nevada, was significantly better compared to Hat Creek, and successful crosscorrelation data were obtained throughout each day. In addition, the digital correlator at OVRO cannot operate in autocorrelation mode, which fortunately was not a limitation.

Receivers and digital correlators of both instruments were configured to record the lines of HCN $`J`$=1–0 (88.6318470 GHz) and HCO+ $`J`$=1–0 (89.1885230 GHz). At OVRO each of the HCN and HCO+ line frequencies were recorded in a single band of 64 channels with a resolution of 125 kHz (0.42 km s-1). At BIMA the lines were recorded in two bands per line, centered on the line frequency, with 128 and 256 channels each and with respective channels widths of 391 and 98 kHz (1.3 and 0.33 km s-1). The autocorrelation mode of the BIMA correlator offers twice the number of channels at half the width for each of the bands, giving velocity resolutions of 0.65 and 0.17 km s-1.

Great care was taken to ensure that both arrays were using the same comet ephemeris at all times. On each of the observing days the most recent orbit solution from the Jet Propulsion Laboratory Solar System Dynamics Group333Available at http://ssd.jpl.nasa.gov/horizons.html was used (Table 1). Changes of the predicted comet’s position between subsequent orbit solutions were less than a few arcseconds, smaller than the synthesized beam of $`12^{\prime \prime }`$ and much smaller than the primary beam of $`143^{\prime \prime }`$ for the BIMA autocorrelation observations. The small day-to-day uncertainty in the comet’s position therefore has no appreciable effect on the data. Because OVRO and BIMA use different coordinate systems (B1950.0 vs. J2000, respectively) and use, respectively, geocentric and topocentric representations for the coordinates of Solar-System objects, much effort was spent to ensure that both instruments actually pointed at the same position on the sky. Note that the parallax of the comet at a geocentric distance of 0.3 AU amounts to $`2^{\prime \prime }`$ over the distance between the two arrays. The topocentric velocity of the comet was also Doppler tracked at both sites, keeping the cometocentric velocity centered at 0 km s-1 in the spectral bands.

The complex gains of crosscorrelation data were calibrated separately for the BIMA and OVRO data, although the same flux and phase calibrators were used for both data. Depending on the position on the sky of the comet, the quasars 0716+714, 0923+392, and 1156+295 were used, with respective fluxes of 3.0, 1.2, and 1.2 Jy determined from contemporaneous measurements of planets at OVRO (Table 1). Using the same calibrator at each observatory ensures that the flux scaling of both data sets is identical. The OVRO data were calibrated using the ‘MMA’ package specific to Owens Valley (Scoville et al., 1993); the BIMA data were calibrated using the ‘MIRIAD’ software package (Sault et al., 1995). The OVRO data, originally stored in FITS format were subsequently converted to MIRIAD format, and both crosscorrelation data sets were added together after regridding to a common spectral resolution. The Appendix describes the data-combination procedure in more detail, including a test case validating the process. The combined data were then inverted and imaged using the MIRIAD software package.

The autocorrelation observations were corrected for atmospheric opacity and fluctuations by switching to an ‘off’ position at $`+20^{}`$ in Right Ascension every 30 seconds. The data were recorded on the main-beam brightness temperature scale in K; subsequent observations of the Galactic SiO maser R Leo in cross- and autocorrelation mode confirmed the flux scaling. A number of bad channels were flagged due to spikes, and second- to fourth-order polynomial baselines were subtracted from the autocorrelation spectra from the individual antennas. Finally the spectra from the individual antennas were added to yield a single autocorrelation spectrum for each of the spectral bands. All processing of the autocorrelation data was carried out within the MIRIAD software package.

## 3 Results

### 3.1 Autocorrelation spectra

The autocorrelation spectra of HCN 1–0 (Fig. Combined BIMA and OVRO observations of comet C/1999 S4 (LINEAR)) and HCO+ 1–0 were inspected using the average data of single, individual days and averaging over two, three, and six consecutive days. In HCN a weak signal was detected in two- and three-day averages between 2000 July 21.85 and 2000 July 24.16. These are the three days leading up to the disruption of the nucleus between 2000 July 24.9 and 25.9 (Kidger, 2000) and coincide with a reported brightening at visible wavelengths, likely indicating a minor outburst event. The three-day average main-beam brightness temperature of HCN 1–0 is $`14\pm 4`$ mK km s-1 in the $`143^{\prime \prime }`$ FWHM beam of the six-meter dishes. The emission is centered around a velocity of $`1`$ km s-1, close to the expected velocity of the comet (nominally at zero), and has a width of $`2`$ km s-1 FWHM. The latter value is uncertain because of the low signal-to-noise of the emission; the hyperfine components at +4 and $`7`$ km s-1 were likewise below the noise level. The reported intensity is a total for all three components, obtained by multiplying the integrated intensity of the single detected component by the expected relative intensities in thermodynamic equilibrium: $`0.2+1.0+0.6=1.8`$. Although the overall signal-to-noise of the detection is low, the three-day averages leading up to (2000 July 21.85–24.16; $`14\pm 4`$ mK km s-1) and following (2000 July 24.95–27.17; $`<10`$ mK km s-1 2$`\sigma `$ upper limit) the comet’s final breakup around 2000 July 23–24 clearly show a dearth of the HCN 1-0 emission. The decrease in HCN $`J`$=3–2 emission observed with the IRAM 30-meter telescope by Bockelée-Morvan et al. (2001) after 2000 July 23.8 supports our detected fall-off of HCN. No emission was detected in any of the HCO+ spectra, to a 2$`\sigma `$ upper limit of 4 mK km s-1 in three-day averages.

### 3.2 Crosscorrelation images

The crosscorrelation data were inspected in a similar way as the autocorrelation data. The data were integrated over a single day, over multiple consecutive days, and over 4-hour time intervals to search for transient features. This was done for the data from the two arrays individually as well as for the combined data set. In all cases the emission was integrated over 2 km s-1 around a cometocentric velocity that was varied from $`10`$ to $`+10`$ km s-1 to account for possible small tracking offsets in the comet’s velocity. The reported flux (limits) refer to the sum of the hyperfine components, added with the appropriate weights (1:5:3).

None of the images obtained in this way revealed any positive detection of HCN 1–0 or HCO+ 1–0 emission. Figure Combined BIMA and OVRO observations of comet C/1999 S4 (LINEAR) illustrates this for HCN 1–0 at an adopted cometocentric velocity of 0 km s-1 and including all OVRO and BIMA data. Using the three-day average over 2000 July 21.85–24.16, over which emission was detected in autocorrelation mode, an upper limit to the HCN 1–0 emission is found of 144 mJy beam-1 km s-1 (2$`\sigma `$) in the $`12^{\prime \prime }\times 9^{\prime \prime }`$ synthesized (naturally weighted) beam, corresponding to 203 mK km s-1 . The corresponding values for HCO+ 1–0 are smaller, at 80 mJy beam-1 km s-1 and 117 mK km s-1, because HCO+ has no hyper-fine splitting to sum over.

## 4 Comparison to model calculations

The previous section showed that HCN 1–0 emission was detected, at low signal-to-noise, in autocorrelation mode but not in crosscorrelation mode. If the emission is unresolved ($`<10^{\prime \prime }`$), the beam-averaged value of $`14\pm 4`$ mK km s-1 detected in the $`143^{\prime \prime }`$ primary beam corresponds to an intensity of $`2.7\pm 0.8`$ K km s-1 in the $`12^{\prime \prime }\times 9^{\prime \prime }`$ synthesized beam of the crosscorrelation data. Since we obtained an upper limit of 203 mK km s-1, it follows that the emission is significantly resolved out in the interferometer beam. Alternative explanations for reduced flux in the crosscorrelation data are decorrelation due to uncorrected phase variations and positional uncertainties of the comet larger than the synthesized beam. These processes are unlikely to affect the data: the measured decorrelation on the phase calibrator is typically less than 20%, and, as mentioned in §2, the position of the nucleus was known to within a few arcseconds from day to day, which is much smaller than the synthesized beam.

Extended emission is expected from molecules in the expanding coma. For a parent species like HCN, the density in the coma follows a Haser distribution,

$$n(r)=\frac{Q}{4\pi V_0r^2}\mathrm{exp}\left(\frac{rr_\mathrm{n}}{r_\lambda }\right),$$
(1)
where $`n`$ is the number density of molecules, $`Q`$ the production rate of that molecule by the nucleus, $`V_0`$ is the expansion velocity in the coma, $`r_\mathrm{n}`$ is the size of the nucleus, and $`r_\lambda `$ is characteristic photodissociation scale length. Bockelée-Morvan et al. (2001) quote $`Q`$(HCN)=$`3\times 10^{25}`$ s-1; we adopt $`V_0=0.9`$ km s-1 and $`r_\lambda =V_0/\beta =3.7\times 10^4`$ km, where $`\beta =2.45\times 10^5`$ s-1 is the HCN photodissociation rate for the Solar spectrum at the time of the observations (Bockelée-Morvan 2001, private communication).

Is a distribution of HCN following eq. (1) consistent with our detection in the $`143^{\prime \prime }`$ autocorrelation beam and our upper limit in the $`12^{\prime \prime }\times 9^{\prime \prime }`$ crosscorrelation beam? Or is an additional source of HCN required, distributed over $`>10^{\prime \prime }`$, e.g., evaporating icy fragments that have broken off the nucleus? Fragments exceeding $`1`$ m, non-uniformly illuminated by the Sun, are expected to undergo non-gravitational forces and accelerate away from the inner ($`10^{\prime \prime }`$) coma (e.g., Desvoivres et al. 1999). Our data uniquely probe the spatial distribution of the HCN 1–0 emission. Bockelée-Morvan (2001, private communication) report HCN 1–0 emission in the $`28^{\prime \prime }`$ IRAM 30-meter beam of $`93\pm 13`$ mK km s-1, averaged over 2000 July 21–23; we include this data point in our analysis. We do not consider the HCN 3–2 emission reported by Bockelée-Morvan et al. (2001) here, because we limit our discussion to the 1–0 transition.

The emission distribution of HCN 1–0 depends on the spatial distribution of the molecules following eq. (1) and on their excitation as function of position. The latter is given by the balance of (de-) excitation through collisions with water and electrons, excitation into vibrationally excited levels through absorption of Solar infrared photons, and de-excitation through spontaneous emission. Bockelée-Morvan et al. (1984) and Crovisier (1987) discuss the excitation of HCN in comets in detail. Instead of carrying out a full excitation calculation, we will use the HCN level populations modeled by Crovisier (1987): this paper presents three cases: water production rate $`Q`$(H2O)=$`2\times 10^{29}`$ s-1 and heliocentric distance $`r_h=1.0`$ AU (model 1); $`Q`$(H2O)=$`2\times 10^{28}`$ s-1 and $`r_h=1.6`$ AU (model 2); and $`Q`$(H2O)=$`1\times 10^{30}`$ s-1 and $`r_h=0.8`$ AU (model 3). None of the parameters of these models are exactly identical to those for comet LINEAR ($`Q`$(H2O)=$`3\times 10^{28}`$ s-1, Bockelée-Morvan et al. 2001; and $`r_h=0.76`$ AU), but we show that especially the *relative* line intensities in the $`12^{\prime \prime }\times 9^{\prime \prime }`$, $`28^{\prime \prime }`$, and $`143^{\prime \prime }`$ beams are comparatively insensitive to the model details and are largely determined by the underlying column density distribution. Although the exact values of the of the relative line intensities in the threee beams depend on the details of the excitation distribution in the coma that results from the balance between fluorescent and collisional processes, values within the same range are found in the extreme cases were either of these processes dominate. The models described above can therefore be used to test if our data agree with a Haser distribution for the column density.

Using a Haser distribution following eq. (1) with the parameters listed above, and the $`J`$=0 and $`J`$=1 HCN level populations of the three models from Crovisier (1987), we calculate the expected HCN 1–0 emission distributions at a distance of 0.374 AU from Earth. We convolve the emission distributions with Gaussians of $`143^{\prime \prime }`$ and $`28^{\prime \prime }`$ FWHM to get the main-beam intensities in these beams. For the crosscorrelation data, we construct synthetic visibilities based on the predicted emission distribution and the actual antenna positions of the data. This explicitly takes into account the spatial filtering of the interferometers and their insensitivity to large-scale structure, as well as the data combination procedure outlined in the Appendix. We process the synthetic visibilities in the same way as the data to obtained integrated-intensity images. We find that no more than 10–30% of the flux is resolved out by the interferometer.

Table 2 lists the observational results in the $`143^{\prime \prime }`$, $`28^{\prime \prime }`$, and $`12^{\prime \prime }\times 9^{\prime \prime }`$ beams, alongside the predicted intensities in these beams for the three models. In terms of absolute intensities the observed values are within a factor of two of the model results, indicating a general agreement between the models and the observations. Table 2 also lists the observed and modeled intensities normalized to the value in the $`28^{\prime \prime }`$ beam: $`F_{\mathrm{auto}}/F_{28}`$ and $`F_{\mathrm{cross}}/F_{28}`$. For the three models, $`F_{\mathrm{auto}}/F_{28}`$ ranges from 0.21 to 0.28. Our observed value, $`0.15\pm 0.05`$, is consistent with the models at the $`1\sigma `$ level. Given the possibility of calibration differences between the IRAM 30-meter data and the BIMA data, we conclude that this indicates that the line emission over $`143^{\prime \prime }`$ is consistent with the expectation of a Haser distribution when compared to a $`28^{\prime \prime }`$ region. Bockelée-Morvan et al. (2001) estimate that only fragments with sizes $`1`$ m are large enough to survive outside the inner few thousand km of the coma. Our $`143^{\prime \prime }`$ result shows that such fragments, if they exist, do not carry a significant contribution of the volatile material ($`50`$%). This is consistent with the results of Magee-Sauer et al. (2001), who report that spatially resolved infrared HCN line observations are consistent with a nuclear source, and of Bockelée-Morvan et al. (2001), who conclude that the evolution of the HCN production rate suggests that HCN is mainly released by icy fragments of small size.

The ratios $`F_{\mathrm{cross}}/F_{28}`$ of the three models show a larger range, 1.57–3.65, compared to $`F_{\mathrm{auto}}/F_{28}`$. This indicates that on these smaller scales, excitation differences exist depending on the heliocentric distance and evaporation rate (= density of collision partners). Our non-detection of emission in the crosscorrelation beam yields a strict upper limit $`F_{\mathrm{cross}}/F_{28}<2.5`$. This is consistent with the range predicted by the models, especially when also taking into account possible decorrelation of our crosscorrelation data up to 20% (see above) and, again, possible calibration differences between the IRAM 30-meter and OVRO+BIMA data. In other words, the non-detection in our combined crosscorrelation data is expected for a Haser-model HCN distribution. This indicates, that to obtain high-spatial resolution (sub) millimeter observations of moderately active comets such as LINEAR, higher sensitivity is required. The Combined Array for Research in Millimeter Astronomy (CARMA), planned to be operational in 2005, will already provide the increase in sensitivity required for a marginal detection of a comet like LINEAR; with the advent of the Atacama Large Millimeter Array (ALMA) such objects will become easily accesible. Needless to say, detailed excitation calculations along the lines of Bockelée-Morvan et al. (1984), Crovisier (1987), and Biver et al. (1999) will be needed to interpret these future observations.

## 5 Conclusions

1. We detected HCN 1-0 emission in autocorrelation mode with the BIMA array at $`14\pm 4`$ mK km s-1 in the $`143^{\prime \prime }`$ beam, and find a 2$`\sigma `$ upper limit of $`144`$ mJy beam-1 km s-1 (203 mK km s-1) in the $`12^{\prime \prime }\times 9^{\prime \prime }`$ synthesized beam of combined crosscorrelation data from OVRO and BIMA. Only upper limits are obtained for HCO+ of $`<4`$ mK km s-1 in autocorrelation and $`<80`$ mJy beam-1 km s-1 in crosscorrelation. HCN emission was detected in the three days preceding (2000 July 21.9–24.2) the breakup of the comet’s nucleus around 2000 July 23–24; no emission in HCN or HCO+ was detected after the breakup.
2. Compared to models of HCN excitation in cometary comae, we find that the relative and absolute integrated intensities in $`143^{\prime \prime }`$, $`28^{\prime \prime }`$, and $`12^{\prime \prime }\times 9^{\prime \prime }`$ are consistent with expectations. This indicates a nuclear source for HCN, and no significant contribution from a spread-out ($`>1300`$ km) distribution of evaporating icy fragments. Fragments exceeding 1 m in size are thought to be able to resist evaporation over such distances, and our data limit their contribution to the volatiles to $`50`$% during the days of 2000 July 21–23 immediately preceding the full disruption of comet LINEAR.
3. We showed that data obtained at the BIMA and OVRO interferometers can be successfully combined, creating a ‘virtual’ CARMA array that increases the signal-to-noise by approximately a factor of $`\sqrt{2}`$ due to the increased collecting array. While for these observations the array configurations were chosen to provide overlapping, short $`uv`$ spacings because of atmospheric phase-stability concerns, virtual ‘CARMA’ observations have the possibility to offer complementary array configurations increasing imaging capabilities. These, however, fall short of what a truly combined array with baselines between the individual BIMA and OVRO antennas can offer: a combined CARMA array is foreseen to be operational by 2005. With the arrival of the Atacama Large Millimeter Array, spatially resolved comet observations at (sub) millimeter wavelengths will become increasingly accesible.

We are indebted to the staffs of the Owens Valley and Hat Creek Radio Observatories for providing excellent support during the observations. Wilson Hoffman and Steve Scott are especially thanked for creating and maintaining the software required to track fast-moving targets such as this comet. The research of M. R. H. in Berkeley was supported by the Miller Institute for Basic Research in Science. This work was partially funded by: NASA NAG5-4292, NAG5-4080, NAG5-8708, and NGT5-0083; NSF AST 96-13998, AST96-13999, AST96-13716, AST96-15608, and AST99-81363; Taiwanese grants NSC 86-2112-M-003-T and 89-2112-M-003-004; and the Universities of Illinois, Maryland, and California, Berkeley. We greatly thank the referee, D. Bockelée-Morvan, for providing constructive criticism of our paper and pointing out ways in which to significantly improve our analysis.

## Appendix A Combining data from different interferometer arrays

The results presented in this Paper include images that are constructed based on data obtained simultaneously at two different millimeter interferometers, BIMA and OVRO. The process of combining such data contains a number of steps where special care needs to be taken to ensure that the resulting images correctly represent the data obtained at the individual arrays. This Appendix describes in detail the data combination process using the MIRIAD software package, and presents simulations that validate the procedure (Fig. Combined BIMA and OVRO observations of comet C/1999 S4 (LINEAR)).

As a first step, even before the data are taken, one needs to ensure that both arrays are using the same source position as phase tracking system, and that the receivers and correlators are configured to record the same line frequency using the same systemic velocity of the object at all times. Especially for a moving target like a comet, these requirements are not trivial when dealing with two separate observatories.

After the data have been taken, the time variation of the complex gains (amplitude and phase) has to be calibrated for the two data sets separately since there is no correlation between the gain variations at the two sites. It is not necessary that both arrays use the same objects, usually quasars, as phase and amplitude calibrators. For the absolute calibration of the flux, however, it can be advantageous to include the same secondary calibrator (a bright quasar), the flux of which is derived from contemporaneous observations of a primary calibrator (planet), preferably with both arrays. This ensures that there are no offsets in flux calibration between the arrays. Especially if the arrays have a different response to spatial scales, such flux offset would affect the resulting combined image.

Before the calibrated data sets can be combined, they need to be converted to a common data format (if not already the same). The OVRO data are usually exported in FITS format, while BIMA data are written in a format particular to the MIRIAD software package. The MIRIAD task ‘FITS’ can convert FITS-files to MIRIAD-files. To successfully convert OVRO data using this task, use ‘options=varwt’ so that the weights of the visibility data written by MMA (the OVRO data calibration package) are properly interpreted as the reciprocal of the noise variance. Further processing in the MIRIAD package does not require any additional conversions: many MIRIAD tasks accept input data sets that, for example, have different spectral gridding as do our OVRO and BIMA data.

Figure 3 illustrates the data combination procedure. To test our combination procedure, this figure also shows the results of combining data with a relative spatial offset: this results in an image with two ‘sources’, each with approximately 50% of the flux because in our observations the sensitivity of the OVRO and BIMA data was about equal.


text:
# 1 Introduction

## 1 Introduction

Linear electron-positron colliders for physics in the Higgs and supersymmetric particle regime have to reach centre-of-mass energies in excess of 500 GeV, and hence high accelerating fields in the accelerating structures are needed to limit the size of the accelerator complex. TESLA  is the only collider project based on superconducting cavities. In the 500 GeV baseline design (TESLA-500) the accelerating field, often called ”gradient” hereafter, amounts to 23.4 MV/m. The 1.3 GHz nine-cell niobium cavities for the TESLA Test Facility (TTF) linac have achieved an average gradient of 26.1 $`\pm `$ 2.3 MV/m at a quality factor $`Q_0110^{10}`$ in the most recent industrial production of 24 cavities. The technology developed for TTF is hence adequate for TESLA-500 but considerable improvements are needed for an upgrade of the collider to 800 GeV (TESLA-800). A detailed description of the present status of the nine-cell cavity layout, fabrication, preparation and tests can be found in .

The TESLA cavities (Figure 1) are made from 2.8 mm thick niobium sheets by deep drawing and electron beam welding. A damage layer of about 100 $`\mu `$m thickness has to be removed from the inner surface to obtain an optimum performance in the superconducting state. For the TTF cavities this has been done so far by chemical etching, see below. Using an alternative preparation method, electrolytic polishing (or electropolishing for short), scientists at the KEK laboratory in Tsukuba (Japan) have achieved gradients of up to 40 MV/m in 1.3 GHz single-cell cavities . This was the motivation for starting an R&D program on electropolishing of single-cell cavities of the TESLA shape. This program has been carried out in collaboration between CERN, DESY and Saclay.

## 2 Niobium Chemistry

Niobium metal has a natural Nb2O5-layer with a thickness of about 5 nm. Below this layer other oxides and sub-oxides can be found . Nb2O5 is chemically rather inert and can be dissolved only with hydrofluoric acid (HF). The strong acids used in the niobium surface treatment like HF, HNO3 or H2SO4 play also an essential role for removing defects from the inner surface of the cavity either by dissolving the foreign material itself or the surrounding niobium. Such defects may be abrasive particles from grinding, imprints from the deep drawing process, niobium protrusions from scratches or dirt particles sticking to the surface. Cleaning the niobium with chemical methods is the most practical way to achieve a high quality superconducting surface. Chemical or electrolytic processes can be applied to niobium resonators if the safety and environmental standards regarding HF-containing acid mixtures are obeyed. Nitrous gases, oxygen and hydrogen are produced during etching or electropolishing. A review of the preparation methods of niobium cavities is given in .

### 2.1 Chemical etching

The sheet rolling of niobium produces a damage layer of about 100 $`\mu `$m thickness which has to be removed in order to obtain a surface with excellent superconducting properties. The standard procedure is chemical etching which consists of two basic steps: dissolution of the Nb2O5 layer by HF and re-oxidation of the niobium by a strongly oxidizing acid . Then the new oxide layer will be dissolved by the HF again. The most frequently used oxidizing agent for niobium is nitric acid (HNO3). The material removal rate in a mixture of HF (40%) and HNO3 (65%) is large, of the order of 30 $`\mu `$m per minute . The reactions are strongly exothermic and can cause a thermal runaway situation. Additionally, large quantities of gases (hydrogen, nitrous gases, HF) are produced. For these reasons the mixture of HF and HNO3 alone is not suitable for the etching of cavities.

To obtain a better process control a buffer substance like phosphoric acid H3PO4 (concentration of 85%) may be added . Furthermore the mixture is cooled below 15C which also reduces the migration of hydrogen into the niobium lattice . The standard procedure with a removal rate of about 1 $`\mu `$m per minute is called buffered chemical polishing333The term “polishing” is not really adequate since the resulting surface roughness is in the $`\mu `$m range and grain boundaries are enhanced.  (BCP) with an acid mixture containing 1 part HF, 1 part HNO3 and 2 parts H3PO4 in volume. Another possibility is to use concentrated sulphuric acid (H2SO4) as a buffer . At TTF, a closed-circuit chemistry system is used in which the acid is pumped from a storage tank through a cooling system and a filter into the cavity and then back to the storage. The gases produced are not released into the environment without prior neutralization and cleaning. The cavities are rinsed with low pressure ultrapure water immediately after the chemical treatment.

By now there exists compelling evidence that the BCP process limits the attainable accelerating fields of multi-cell cavities to about 30 MV/m even if niobium of excellent thermal conductivity is used. Typical accelerating gradients achieved with BCP in defect-free cavities are shown in figure 2. Measurements on etched one-cell cavities in other laboratories yielded similar results . Etched cavities exceeded 30 MV/m only in rare cases .

The etching process is accompanied with undesirable effects such as migration of hydrogen into the bulk niobium and strong grain boundary etching. The first effect can be reduced by cooling the acid below 15C . In principle, grain boundary etching could be suppressed by using an acid mixture with a high etching rate. This, however, would not be advisable in the nine-cell TESLA cavity with an inner surface of about 1 m2 since the large amount of produced heat would speed up the reaction even more and preclude a well-controlled etching process.

### 2.2 Electrolytic polishing

An alternative method to etching is electrolytic polishing or electropolishing (EP) in which the material is removed in an acid mixture under the flow of an electric current. Sharp edges and burrs are smoothed out and a very glossy surface can be obtained. The electric field is high at protrusions so these will be dissolved first. On the other hand, the field is low in the grain boundaries and little material will be removed here while in the BCP process strong etching is observed in the boundaries between grains.

Electropolishing of niobium cavities has been known for 30 years. The most widely used electrolyte is a mixture of concentrated HF and concentrated H2SO4 in volume ratio of 1:9 . A pulsed electric current was used in a horizontal EP setup for superconducting niobium cavities prepared at CERN in collaboration with Karlsruhe in 1979 . A continuous method for horizontal EP has been developed at KEK in 1989 .

The chemical processes are as follows :

$`2\text{Nb}+5\text{SO}_4^{}+5\text{H}_2\text{O}`$ $``$ $`\text{Nb}_2\text{O}_5+10\text{H}^+`$
$`+5\text{SO}_4^{}+10\text{e}^{}`$
$`\text{Nb}_2\text{O}_5+6\text{HF}`$ $``$ $`\text{H}_2\text{NbOF}_5`$
$`+\text{NbO}_2\text{F}0.5\text{H}_2\text{O}+1.5\text{H}_2\text{O}`$
$`\text{NbO}_2\text{F}0.5\text{H}_2\text{O}+4\text{HF}`$ $``$ $`\text{H}_2\text{NbF}_5+1.5\text{H}_2\text{O}`$
Most systems for the EP of accelerator cavities are horizontal as shown in figure 3. The advantage is that the produced gases (mainly hydrogen) are rapidly removed from the wetted niobium surface . Gas bubbles sticking on the surface could lead to etching pits. In a vertical setup these bubbles would move slowly upwards and create axial wells. A drawback of the horizontal arrangement is that the cavity has to be rotated. There is some difficulty to achieve a leak tight rotary sleeve for the acid mixture. In addition, the removal rate is reduced by a factor of two since the surface is immersed only half of the time in the acid to allow the hydrogen gas to escape through the upper part of the beam tube.

### 2.3 Comparison of etched and electropolished surfaces

Micro-graphs of BCP and EP treated niobium samples are compared in figure 4. One can see that EP smoothes out the grain boundaries far better than BCP. The average roughness of chemically etched niobium surfaces is of the order of $`1\mu `$m  while EP surfaces are at least one order of magnitude smoother, see figure 5. It should be noted, however, that the measured unevenness of the surface depends on the scan length of the atomic force microscope. On a nanometer scale the roughness of BCP and EP treated samples is comparable. This can be understood because at such small lengths the measurement takes place on a single niobium grain. The main difference between EP and BCP is the smoothening of the ridges at grain boundaries which is measured only when the scan length is sufficiently large.

On etched surfaces the step height at grain boundaries can be a few $`\mu `$m. Particularly large steps may occur at the electron-beam welds at the cavity equator. For non-optimized weld parameters steps as high as 30 $`\mu `$m have been observed . The steps at grain boundaries and in the weld seams may lead to a magnetic field enhancement and a premature breakdown of superconductivity . This interpretation is corroborated by the observation that several chemically etched cavities exhibited a quench in the region of high magnetic field near the equator weld without any indication of contamination by foreign materials. It was possible to shift the quench position along the equator region by applying a new chemical etching, but the maximum achievable field was not improved .

An important question for electropolishing is of course how much material has to be removed to achieve a smooth surface. On electropolished samples the roughness drops below 1 $`\mu `$m after 150 $`\mu `$m of material have been removed. This corresponds to the thickness of the damage layer which has to be removed anyway. Therefore an electropolishing of at least 100 $`\mu `$m is a reasonable choice both for surface smoothening and damage layer removal.

### 2.4 Description of the EP systems

#### 2.4.1 Half-cell electropolishing

In the framework of the DESY-CERN-Saclay collaboration a setup for the EP of half-cells was developed at CERN. The goal was to investigate the viability of a novel cavity fabrication and preparation concept: removal of most of the damage layer by EP on simple subunits and assembly of the multi-cell cavity from pre-polished half-cells. An electrolyte developed at CERN  was used (see table 1) with a removal rate of 1 - 2 $`\mu `$m per minute. The electrode was made from copper. An electropolishing of about 100 $`\mu `$m was applied to the half cells. Then 4 one-cell cavities were electron-beam welded at an industrial company with a standard electron-beam welding apparatus (vacuum $`10^5`$ mbar), and 11 cavities were welded at CERN under much better vacuum conditions (pressure in the welding chamber below $`10^6`$mbar).

Two of the cavities were tested immediately after welding and high-pressure water rinsing. Not unexpectedly, the results were poor: the cavities reached a maximum accelerating field of just 8 - 10 MV/m with a strong decrease in quality factor. A slight chemical etching which would have been sufficient to dissolve a possible copper deposition on the niobium surface cured the Q-degradation, however electron field emission was still present at low excitation implying that another contaminant besides copper was present. After a chemical etching (BCP) of 20 - 50 $`\mu `$m both cavities reached the standard performance (25 MV/m) of BCP-treated cavities. The obvious explanation is that welding deposits on the inner surface of the cavity must be removed by a BCP or EP of 20 - 50 $`\mu `$m to arrive at a defect-free surface.

The conclusion to be drawn from these observations is that the fabrication of a multi-cell cavity from prepolished half-cells is indeed a viable concept which would substantially reduce the amount of EP needed in the finished resonator.

### 2.5 Electropolishing of complete cavities

A setup for the electropolishing of single-cell cavities was built at CERN. The parameters were chosen similar to those of the successful KEK system  and are summarized in table 1. The essential difference is the cathode material: pure aluminum at KEK and copper at CERN. Both materials have their respective advantages and disadvantages. The copper electrode requires some care, see below. The cathode is surrounded with tube made from a porous PTFE444Polytetrafluoroethylene, for example Teflon® cloth to keep the electrolytically produced hydrogen gas away from the niobium surface.

Except for the cathode all components of the EP system are made from plastic material which is inert against the aggressive acid mixture. All parts in contact with the HF should be made from fluoroplastics such as PFA555Polyperfluoroalkoxyethylene, PVDF666Polyvinylidene Fluoride or PTFE. The EP system is placed in a vented area where the exhaust gases are pumped through a neutralization system to avoid environmental hazards. The acid mixture is contained in a closed-circuit. The electrolyte is stored in containers cladded with Teflon which can be cooled by water flowing through Teflon-covered piping. The volume above the electrolyte is filled with dry nitrogen to avoid water vapour absorption by the strongly hygroscopic H2SO4. The acid mixture is pumped with a membrane pump through a cooler and a filter with 1 $`\mu `$m pore size before it reaches the inlet of the cathode. Then the electrolyte flows through the cathode to the center of the cell (see figure 3(b)). The acid returns to the storage tank via an overflow.

The cavity is installed horizontally together with the cathode and then the lower half is filled with the electrolyte which attacks the niobium only very slowly when no voltage is applied (etch rate less than 1 nm per hour). After the equilibrium filling level has been reached, the rotation is switched on and a leak check is done. Then the current is switched on and the current-voltage characteristic is measured. At a certain voltage (see table 1) current oscillations set in which are an indication that two alternating processes take place: dissolution of Nb2O5 by HF and re-oxidation by H2SO4. The best polishing results are obtained for an oscillation amplitude of 10 - 15% around the mean value. The temperature of the acid has to be around 30 - 35C during the EP. Temperatures above $`40^{}`$C must be avoided as they result in etching pits on the surface.

When the desired amount of material has been removed, the current is switched off. The rotation is stopped and the cavity is put into vertical position to release the acid mixture. After rinsing with pure water the cavity and the electrode are dismounted. Another low-pressure water rinsing follows and the wet cavity is taken into a glove-box with nitrogen atmosphere. Here the cavity is rinsed with an ultrapure high-pressure water jet to remove remaining chemical residues. A rinsing with filtered ethanol follows to speed up the drying process. The cavity is then stored overnight in a vacuum of $`10^3`$ mbar. Then the cavity is either rinsed with high-pressure water for the performance tests at CERN or filled with nitrogen gas and sent to Saclay or DESY for high-pressure water rinsing and tests.

It should be mentioned that the copper electrode used in the CERN system has a disadvantage. While the current is flowing, the copper becomes passivated but after the current has been switched off copper ions can be dissolved and precipitate on the cavity surface. It has been observed in a test that this leads to an increased residual resistance and to strong field emission at low field ($`<10`$MV/m) . To remove the copper deposition the cavities are rinsed in sequence with HNO3 (60 %), water, HF and water again. These steps are repeated twice. The HNO3 dissolves the copper and oxidises the niobium. Hydrofluoric acid removes the oxide layer and several contaminants as indicated by our surface studies. The overall material removal by this procedure is only about 10 nm and does not deteriorate the smoothness of the electropolished surface.

## 3 Measurements on electropolished cavities

### 3.1 First tests

In their first tests the cavities which were electropolished in the CERN system showed an unexpected performance limitation: the excitation curves exhibited a strong degradation in quality factor at high field as can be seen in figure 6. Field emission of electrons could be excluded as an explanation for the performance degradation since neither X rays nor secondary electrons were observed. Temperature mapping revealed a global heating in the areas of high magnetic fields around the equator indicating a rapid increase of the microwave surface resistance towards high microwave fields.

In contrast to this observation, the excitation curves of the EP-treated cavities at KEK showed only a moderate drop of the quality factor and reached much higher accelerating fields, see figure 7 .

At Saclay a similar strong $`Q`$ degradation was observed for chemically etched cavities , see figure 8. The temperature mapping revealed a global heating of the surface similar to that in figure 6b.

It was discovered that the unsatisfactory performance could be considerably improved by applying a moderate thermal treatment to the finished cavity . This will be called bakeout (or in situ bakeout) hereafter. The procedure at Saclay was as follows. After the last high-pressure water rinsing the cavities were evacuated and then heated up to 170C for 70 hours. The remarkable observation was that this low-temperature baking improved the quality factor at the highest field by nearly a factor of 3. This result was confirmed in other tests on BCP-treated cavities. However, tests revealed also that the baking procedure did not increase the maximum field level in chemically etched cavities .

### 3.2 Application of low-temperature bakeout to EP cavities

Building on the experience with BCP cavities at Saclay and with EP cavities at KEK (where the bakeout at 85–100 C had been part of the standard preparation) it was decided to apply the bakeout to the EP cavities of the CERN-DESY-Saclay collaboration. Figure 9 shows that a dramatic improvement in performance is achieved: several single cell cavities reach now accelerating gradients of up to 40 MV/m with quality factors above $`510^9`$. This behaviour is very similar to the observations made at KEK.

The heat treatment of niobium cavities is a well-known method to improve the performance. The temperature is in general quite high though: the TTF cavities are heat treated in an UHV furnace at 800C for hydrogen degassing and stress annealing and afterwards at 1400C for increasing the thermal conductivity of the bulk niobium777During this furnace treatment a titanium getter layer is evaporated onto the niobium surface which protects the niobium from being oxidized by the residual oxygen in the furnace atmosphere.. The surprising observation made with the bakeout effect is that a thermal treatment at such a low temperature, where the diffusion of any gases dissolved in the niobium lattice is extremely slow, has a major influence on the high-gradient performance. It it obvious that only the thin surface layer which is essential for the microwave superconductivity can be modified by the bakeout.

### 3.3 High-field performance at different helium temperatures

Figure 10 shows the excitation curves of an electropolished cavity before and after the bakeout at helium bath temperatures between 1.5 K and 2.2 K. At low field the quality factor exhibits the well-known temperature dependence which is caused by the exponential temperature dependence of the BCS surface resistance (see below). However, the maximum achieved gradient and the corresponding $`Q_0`$ value are almost independent of the bath temperature as long as the helium coolant is in the superfluid state. These results are consistent with thermal model calculations . Only when the temperature is in the vicinity of the Lambda-point (2.17 K) of liquid helium, where the transition from the superfluid to the normal-fluid phase takes place, a degradation of the quality factor due to insufficient cooling can already be seen at low field.

### 3.4 High temperature heat treatments of electropolished cavities

Electrolytic polishing generates hydrogen which can penetrate into the niobium lattice. If too much hydrogen is dissolved the danger exists that niobium hydride compounds are formed if the material is exposed to temperatures around 100 K for an extended period. These Nb-H compounds have very high microwave losses and may lead to a reduction of the quality factor by one or two orders of magnitude. This unfortunate effect has been named Q-disease.

The typical hydrogen concentration in etched niobium is in the order 100 - 400 atomic ppm hydrogen in the bulk material and 4 atomic percent in a 50 nm surface layer . A high temperature treatment at 800C is very effective in reducing the hydrogen concentration to 3 atomic ppm in the bulk and $``$ 1 atomic% in the surface layer.

The hydride formation depends on the dwell time the cavity spends in the dangerous region around 100 K. With a fast cooldown from 300 K to 4.2 K within 1 hour, hydride formation will be reduced to an acceptable level.

The cavities studied in this work have been made from electropolished half-cells. In the half-cell EP system the niobium is not protected from the hydrogen gas by a membrane. Therefore the danger of quality factor degradation was to be expected. A dedicated experiment was carried out to this end, the results are shown in figure 11. The electropolished and baked cavity was tested after a fast cooldown and showed the usual good performance. Then the cavity was exposed to a temperature of 100 K for 2 days. In the new test the quality factor was very low and showed the typical strong field dependence described in . In the next step the cavity was heated in a vacuum furnace to 800C to remove the hydrogen from the bulk material. A short electropolishing of about 20 $`\mu `$m was done to clean the surface from dirt which might have been introduced during the furnace treatment. Again the cavity was tested after fast cooldown and showed the initial good performance. The exposure of the cavity to the dangerous temperature of 100 K for 2 days was repeated. This time, however, no $`Q`$ degradation was observed, indicating that the formation of hydrides was strongly suppressed by the furnace treatment. This is an important and encouraging result which implies that no hydrogen will be introduced into the niobium during a short EP if precautions are taken (PTFE-cloth around the cathode) to keep the gas away from the niobium surface.

Seven cavities received the furnace treatment at 800C after electropolishing. Figure 12 shows the maximum gradients before and after the 800C treatment. The average gradient rises from 35 MV/m to 35.4 MV/m.

Finally, one electropolished cavity was subjected to the $`1400^{}`$C furnace treatment with titanium getter which is part of the standard preparation of the TTF cavities (see ). The idea was to check whether the increased thermal conductivity would have a similar beneficial effect on EP cavities as on BCP cavities. The residual resistivity ratio RRR increased from 300 to 500. The cavity reached the same maximum field of 35 MV/m as before the 1400C treatment, but at a slightly higher $`Q_0`$. More tests are needed to decide whether the $`1400^{}`$C furnace treatment would constitute a significant improvement of the electropolished multi-cell cavities of the TESLA-800 collider.

### 3.5 Comparison of etching and electropolishing

The test results on etched single-cell cavities are consistent with the performance achieved in the TESLA nine-cell cavities. The average gradient achieved in the one-cell cavities without any heat treatment is 24 MV/m. If one evaluates the single-cell performance of those etched nine-cell cavities which were tested after the 800C heat treatment an average gradient of 23.5 MV/m is obtained. The 1400C heat treatment which is routinely applied to the BCP-treated TTF cavities raises the average gradient to more than 25 MV/m .

A summary of the performance after etching and electropolishing is given in figure 13. The average gradient of the electropolished cavities is 35.7 MV/m. This is convincing evidence that electropolishing leads to substantially higher accelerating fields. One reason is certainly that magnetic field enhancements at grain boundaries are avoided. Moreover, the smoother EP surface contains fewer defects and is easier to clean from chemical residues.

### 3.6 Influence of the electron-beam welding on the cavity performance

The TESLA cavities are made from deep-drawn niobium half-cells which are joined by electron-beam (EB) welding. The quality of the EB weld seam at the equator is crucial for the maximum field achievable in the cavities. The one-cell production was made in two different EB welding machines. The first batch of 4 cavities was produced by an industrial company with experience in niobium cavity welding (ACCEL Instruments, Bergisch-Gladbach). The second batch of 11 cavities was welded at CERN with an EB welding machine equipped with a stainless steel vacuum chamber conforming to ultra high vacuum standards. Whereas the vacuum in the industrial production is of the order of a few times $`10^5`$ mbar, the pressure in the CERN system could be lowered down to $`310^7`$ mbar by using a cryopump and a longer time for pumping down.

As the residual gas in the machine is partly gettered by the molten niobium in the weld seam, the residual resistivity ratio RRR in this region can be degraded. For an initial RRR $``$ 300 the reduction can be of the order of 20% for a vacuum of $`10^5`$ mbar 888For TIG (Tungsten-Inert-Gas) welds the reduction can be very large (RRR $``$ 400 is reduced to RRR $``$ 120) . In addition, defects can be introduced using this welding technique.. It can be expected that the better vacuum conditions lead to an improved cavity performance. The average gradient of the CERN batch is in fact 3.8 MV/m higher as compared to the cavities prepared under standard vacuum conditions (see figure 14) 999All cavities were treated in a vacuum furnace at 800C to remove hydrogen from the niobium..

The CERN EB machine is equipped with a 60 kV electron gun without the possibility to raster the electron beam in a rhombic pattern across the weld which produces the smoothest underbead of the weld seam . The results in figure 14 indicate that the electron-beam welding facility should combine a state-of-the-art electron gun and a UHV vacuum system. Such an electron-beam welding facility is presently being commissioned at DESY.

### 3.7 Influence of bakeout on the microwave surface resistance

In figure 15 the microwave surface resistance before and after the low-temperature bakeout is plotted against the inverse temperature. According to the BCS (Bardeen-Cooper-Schrieffer) theory the surface resistance should drop exponentially when plotted against $`T_c/T`$. For temperatures $`T<\frac{T_c}{2}`$ and an energy of the microwave photons $`hf`$ much less than the energy gap $`\mathrm{\Delta }`$, the surface resistance can be approximated by:

$$R_{BCS}(T,f)=A\frac{f^2}{T}\mathrm{exp}\left(\frac{\mathrm{\Delta }}{k_BT}\right)$$
(1)
The factor $`A`$ depends on material parameters like the coherence length $`\xi `$, the electron mean free path $`\mathrm{}`$, the Fermi velocity $`v_F`$ and the London penetration depth $`\lambda _L`$. The BCS resistance of niobium at 1.3 GHz is about 600 n$`\mathrm{\Omega }`$ at 4.2 K and 10 n$`\mathrm{\Omega }`$ at 2 K.

A refined expression, derived from the two-fluid model of superconductors, is :

$$R_{BCS}(T,f)\sigma _{nc}\lambda _{eff}^3\frac{f^2}{T}\mathrm{exp}\left(\frac{\mathrm{\Delta }}{k_BT}\right)$$
(2)
where $`\sigma _{nc}`$ is the conductivity due to the normal-conducting component and $`\lambda _{eff}=\lambda _L\sqrt{1+\xi /\mathrm{}}`$. The BCS resistance is hence proportional to the mean free path $`\mathrm{}`$ of the normal-conducting electrons and does not assume its minimum for extremely pure niobium ($`\mathrm{}`$ very large) but rather for somewhat “dirty” Nb with $`\mathrm{}\xi `$. If we make the assumption that the bakeout creates a “dirty” surface layer with a shorter mean free path than in the bulk niobium, then the residual resistance will be larger and the BCS resistance will be lower than in the unbaked cavity. Figure 15 shows that this is indeed the case. The BCS part of the surface resistance is reduced by a factor of $`1.5`$ due to the bakeout while the residual resistance increases. Similar observations were made by Kneisel .

It should be noted that the resistance data are derived from measurements of the cavity quality factor at low field (a few MV/m). The peculiar high-field behaviour of unbaked cavities, namely the strong $`Q`$ degradation, cannot be explained in terms of the BCS resistance because this quantity is independent of the magnitude of the electromagnetic field in the cavity. On the contrary, a strongly field-dependent resistance would be needed to account for the rapid decrease of $`Q_0`$ towards large gradients. This effect is not yet understood.

## 4 Conclusions and outlook

Electropolished bulk niobium cavities offer accelerating fields of more than 35 MV/m and are hence suitable for the upgrade of the TESLA collider to 800 GeV. The surface preparation by electropolishing has a definite advantage in comparison with chemical etching: a much smoother surface is obtained and the field enhancement at grain boundaries is avoided. This magnetic field enhancement may be the origin of a local premature breakdown of superconductivity in etched cavities. The physical and physico-chemical principles of the bakeout effect are not well understood yet and require further research. Meanwhile, the first multi-cell cavities have been electropolished in a collaboration between KEK and DESY. The results will be published in a separate paper.

## 5 Acknowledgments

The authors would like to express their gratitude to O. Aberle, S. Ehmele, S. Forel, A. Lava, A. Insomby, R. Guerin, M. Kubly, the DESY group MKS3, K. Twarowski, J.-P. Poupeau, Y. Gasser for their invaluable help in measurements, cavity fabrication, chemical treatments and handling of cavities.

We are particularly grateful to K. Saito and E. Kako for their advice and information on electropolishing. We would like to thank E. Haebel for the discussions on RF superconductivity and help with the measurement system at CERN. We would like to thank P. Kneisel for stimulating discussions on the bakeout effect.