text:
# PROGRESS IN PERTURBATIVE COLOR TRANSPARENCY

## I Overview

Several experiments indicate that color transparency  and nuclear filtering  have been observed at large nuclear number $`A`$. The first color transparency experiment of Carroll et al  convincingly showed that interference effects in proton-proton scattering were filtered away in nuclear targets. Fits to the attenuation cross section in nuclear targets show values significantly below the Glauber theory values . The FNAL E-665 experiment  also proved consistent with filtering effects . The observation of increasing longitudinal final state polarization in $`\gamma ^{}A\rho A`$ as a function of $`Q^2`$ is noteworthy. We still await confirmation of predicted longitudinal polarization increasing as a function of $`A`$ .

Electron beam experiments remain controversial, with few signals of interesting $`Q^2`$ dependence . A basic feature of $`\gamma ^{}`$-initiated reactions is that most events are knocked out from the back side of the nucleus. This minimizes the resolving power of such experiments to measure the size of propagating states. The $`A`$ dependence is a particularly useful tool  to measure effective attenuation cross sections. O’Neill et al  showed that effective attenuation cross sections extracted from $`A(e,e^{}p)`$ SLAC data were smaller than Glauber theory calculations by a statistically significant amount. However, the precision of the data  was insufficient to establish a large effect, and model dependence in the choice of the normalization of hard scattering is another complication. Reports on new $`(e,e^{}p)`$ beam experiments from CEBAF are expected shortly.

Progress on the theory front has come from looking deeper at the basic factorization methods , and doing the work of labor-intensive calculations . The asymptotic factorization scheme of Lepage and Brodsky  is inadequate. An integration over the transverse separation of quarks is needed in the description. We call this “impact parameter factorization”, which is needed to describe the interactions with the nuclear target, which otherwise vanish prematurely in the pure short-distance scheme. The impact parameter method was originally found necessary to regulate Landshoff and Sudakov effects in $`pp`$ scattering . We adapted it to describe color transparency and nuclear filtering . Impact-parameter factorization has subsequently become very popular for the description of free-space form factors , which remain controversial  at laboratory $`Q^2`$ values.

## II The Calculations

Elsewhere  we report details of calculations of hard-pion knockout, $`\gamma ^{}A\pi A^{}`$, and hard nucleon knockout, $`\gamma ^{}ApA^{}.`$ These are exploratory concept studies, designed to see how $`pQCD`$ predicts color transparency and filtering with few parameters. Since all the details except for experimental acceptances have been incorporated, the calculations are also fully quantitative predictions of the type needed to compare to experiments.

The case of pionic transparency deserves special mention. First, the pion is the cleanest theoretical laboratory one would desire. A short-distance wave function is known from experiments on pion decay , without relying on the sometimes circular logic of schemes such as $`QCD`$ sum rules. Second, a pion is ultra-relativistic at energies as low as a few $`GeV`$. This helps strengthen the approximations made in $`pQCD`$. Finally, the pion has only two quarks in its valence state, and one transverse separation $`b`$, reducing the complexity of the calculations.

Working in configuration (impact-parameter $`b`$) space the expression for a $`\gamma ^{}meson`$ form-factor becomes:

$$F_\pi (Q^2)=𝑑x_1𝑑x_2\frac{d^2\stackrel{}{b}}{(2\pi )^2}𝒫(x_2,b,P_2,\mu )\stackrel{~}{H}(x_1,x_2,Q^2,\stackrel{}{b},\mu )𝒫(x_1,b,P_1,\mu ).$$
(1)
Here $`𝒫(x,b,P,\mu )`$ represent the Fourier transforms of the wave functions, including Sudakov factors; $`\stackrel{~}{H}(x_1,x_2,Q^2,\stackrel{}{b},\mu )`$ represents the hard scattering kernel from perturbation theory. The impact parameter $`\stackrel{}{b}`$ is conjugate to $`\stackrel{}{k}_{T1}\stackrel{}{k}_{T2}`$, $`\mu `$ is the renormalization scale, and $`P_1`$, $`P_2`$ are the initial and final momenta of the meson.

The nuclear medium modifies the quark wave function by an interaction kernel $`f_A`$, which is called the nuclear filtering amplitude. An eikonal form  appropriate for $`f_A`$ is: $`f_A(b;B)=exp(_z^{\mathrm{}}𝑑z^{}\sigma (b)\rho (B,z^{})/2)`$. Here $`\rho (B,z^{})`$ is the nuclear number density at longitudinal distance $`z^{}`$ and impact parameter $`B`$ relative to the nuclear center. We parametrize $`\sigma (b)`$ as $`kb^2`$ for our calculations. Finally, we must include the probability to find a target at position $`B,z`$ inside the nucleus. Then the wave functions $`𝒫_A`$ appropriate for the nuclear target are

$$𝒫_A(x,b,P,\mu )=f_A(b;B)𝒫(x,b,P,\mu ).$$

Putting together the factors, the process of knocking out a hadron from inside a nuclear target has an amplitude $`M`$ given by

$`M`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}d^2B(\mathrm{\Pi }dx_id^2b_i){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑z\rho (B,z)\times F_\pi (x_1,x_2,b,Q^2)\times f_A(b,B)`$ (2)
The analysis for the proton is similar but vastly more complicated. A 9 dimensional integration over the various $`x_i,b_i`$ coordinates is performed by Monte Carlo. Sudakov effects are set to depend on the maximum of the three quark separation distances, $`b_{max}=max(b_1,b_2,b_3)`$.

We find that the physics is not described by a free-space hard scattering, followed by some model of propagation with or without “expansion”, which is the ansatz of most competing groups. The integrations over the transverse quark variable extend over the whole volume of the nucleus. There is no easy decoupling into a simple product of “hard” times “nuclear” effects. Color transparency truly probes the internal structure of hadrons.

We found uncertainties in the nuclear correlations at the $`10\%`$ level to be a major concern, in some cases exceeding the theoretical uncertainties from the rest of the calculation . Primary new results include a discovery that the $`\mathrm{𝑠𝑙𝑜𝑝𝑒}`$ of the transparency ratio with $`Q^2`$ depends strongly on the $`x`$ dependence of wave functions (distribution amplitudes). This mysterious effect was traced to the fact that central wave functions are more effective in maintaining short distance. End-point dominated distributions tend to exacerbate long-distance effects, which are found not to produce transparency. Nuclear filtering was observed to depend on the choice of wave functions as consistent. Thus both the slope of transparency ratios, and the magnitude of effective attenuation cross sections extracted from data, are probes of $`x`$ dependence. Extensive details and nearly a dozen plots are given in Ref. .

Acknowledgments: This work was supported by BRNS grant No. DAE/PHY/96152, the Crafoord Foundation, the Helge Ax:son Johnson Foundation, DOE Grant number DE-FGO2-98ER41079, the KU General Research Fund and NSF-K\*STAR Program under the Kansas Institute for Theoretical and Computational Science.


format: text

text:
# Evolution of Stellar Collision Products in Globular Clusters – II. Off-axis Collisions.

## 1 INTRODUCTION

Globular clusters have, until recently, been studied under the auspices of either stellar evolution or stellar dynamics. Stellar evolutionists often treat clusters as collections of isolated stars all with the same age and primordial chemical abundances, while stellar dynamicists often treat clusters as a system of point masses moving under the influence of their mutual gravities. However, it is becoming clear that the dynamics of a cluster can affect its populations, and vice versa (Bailyn, 1995; Portegies Zwart et al., 1999; Joshi, Nave & Rasio, 2000). Therefore, in order to model globular clusters effectively, we need to consider both stellar evolution and stellar dynamics, and especially their influence on each other.

One situation of particular interest in which stellar dynamics and stellar evolution interact involves direct stellar collisions. These collisions, especially those involving binary stars, can have a significant impact on the energy budget of a cluster. In particular, they can destroy and change the properties of hard primordial binaries. Dynamical interactions involving hard primordial binaries are thought to be the most important source of heating that can support a cluster against core collapse (see, e.g., Goodman & Hut, 1989; Gao et al., 1991; Rasio, 2000). In addition, it has been suggested that, through the indirect heating by mass loss associated with the stellar evolution of collision products in the cluster core, collisions could provide a mechanism for supporting a cluster against gravothermal collapse after exhaustion of its primordial binaries (Goodman & Hernquist, 1991). The products of stellar collisions between main-sequence stars have the additional useful property that they are potentially observable in the color-magnitude diagrams of clusters as blue stragglers. Since the population of blue stragglers is determined by the dynamics of the cluster, we can use the blue stragglers to investigate the dynamical state of the cluster. For example, Ferraro et al. (1999) have suggested that the large population of blue stragglers in M80 implies that the cluster is close to core collapse. In order to use blue stragglers to study globular clusters, however, we need to understand clearly the evolution and observable properties of collision products.

In previous work (Sills et al., 1997, hereafter Paper I), we used the results of Smoothed Particle Hydrodynamics (SPH) simulations of stellar collisions as starting models for stellar evolution calculations, and produced detailed evolutionary models for head-on (i.e., zero impact parameter) collision products in globular clusters. However, collisions between stars are, of course, never exactly head-on. Lombardi, Rasio & Shapiro (1996) calculated collisions with non-zero impact parameters, and some of the resulting off-axis collision products were rotating close to their break-up speeds when they returned to dynamical equilibrium. The effects of rotation on the subsequent evolution of these collision products will be substantial. These effects have never before been included in any theoretical study of blue straggler evolution.

In this paper, we present four new SPH simulations of off-axis collisions, and use the results to initiate stellar evolution calculations of the rotating collision products. In §2, we outline the method and results of the SPH calculations. We summarize our technique for combining the results of SPH simulations of stellar collisions with detailed evolutionary models, and discuss the modifications we have made for rotating stars. §3 presents the implementation of rotation used in the Yale Rotational Evolution Code (YREC), and the results of the evolution calculations. We discuss the implications of our results for blue straggler creation mechanisms and for using blue stragglers as tracers of the dynamical evolution of globular clusters in §4.

## 2 HYDRODYNAMIC SIMULATIONS OF OFF-AXIS STELLAR COLLISIONS

### 2.1 Implementation of the SPH Code

Our hydrodynamic simulations were performed using a new, parallel version of the SPH code described in Paper I. This code was originally developed by Rasio (1991) specifically for the study of hydrodynamic stellar interactions such as collisions and binary mergers (see, e.g., Rasio & Shapiro, 1991, 1992, 1994). SPH is a Lagrangian hydrodynamics technique in which pressure gradient forces are calculated by kernel estimation directly from particle positions (for a review, see Monaghan, 1992). In SPH simulations we solve the equations of motion of a large number $`N`$ of Lagrangian fluid particles moving under the influence of hydrodynamic forces as well as the fluid self-gravity (in this paper $`N`$ is as large as $`1.05\times 10^5`$). For a basic overview of the SPH equations, see Rasio & Lombardi (1999).

Local densities and hydrodynamic forces at each particle position are calculated by smoothing over $`N_N`$ nearest neighbors. The size of each particle’s smoothing kernel is evolved in time to keep $`N_N`$ close to a predetermined optimal value. For the high-resolution calculations in this paper the optimal number of neighbors was set at $`N_N=100`$. Neighbor lists for each particle are recomputed at every iteration using a linked-list, grid-based parallel algorithm. Note that neighborhood is not a symmetric property, as one particle can be included in another’s neighbor list but not vice versa. Pressure and artificial viscosity (AV) forces are calculated by a “gather-scatter” method, looping over all particles in turn, gathering the force contribution on each particle from its neighbors, and scattering out the equal and opposite force contribution on the neighbors from each particle.

The fluid self-gravity in our code is calculated by an FFT-based convolution method. The density field is placed on a 3D grid by a cloud-in-cell method, and convolved with a kernel function calculated once during the initialization of the simulation. We use zero-padding of our grids to obtain correct boundary conditions for an isolated system, at the expense of memory storage. Gravitational forces are calculated from the gravitational potential by finite differencing on the grid, and then interpolated for each particle using the same cloud-in-cell assignment.

A number of variables are associated with each particle $`i`$, including its mass $`m_i`$, position $`𝐫_i`$, velocity $`𝐯_i`$, entropic variable $`A_i`$ and numerical smoothing length $`h_i`$. The entropic variable $`A`$ is a measure of the fluid’s compressibility and is closely related (but not equal) to specific entropy: for example, both $`A`$ and specific entropy are conserved in the absence of shocks. For convenience, we will refer to $`A`$ as entropy. In this paper we adopt an equation of state appropriate for a monatomic ideal gas: that is, the adiabatic index $`\gamma =5/3`$ so that $`P_i=A_i\rho _i^{5/3}`$, where $`P_i`$ and $`\rho _i`$ are the density and pressure of particle $`i`$.

The rate of increase, due to shocks, for the entropy $`A_i`$ of particle $`i`$ is given by

$$\frac{dA_i}{dt}=\frac{\gamma 1}{2\rho _i^{\gamma 1}}\underset{j}{}m_j\mathrm{\Pi }_{ij}(𝐯_i𝐯_j)_iW_{ij},$$
(1)
where the summation is over neighbors and $`W_{ij}`$ is a symmetrized smoothing kernel. We adopt the form of AV proposed by Balsara (1995):

$$\mathrm{\Pi }_{ij}=\left(\frac{p_i}{\rho _i^2}+\frac{p_j}{\rho _j^2}\right)\left(\alpha \mu _{ij}+\beta \mu _{ij}^2\right),$$
(2)
where

$$\mu _{ij}=\{\begin{array}{cc}\frac{(𝐯_i𝐯_j)(𝐫_i𝐫_j)}{h_{ij}(|𝐫_i𝐫_j|^2/h_{ij}^2+\eta ^2)}\frac{f_i+f_j}{2c_{ij}}\hfill & \text{if }(𝐯_i𝐯_j)(𝐫_i𝐫_j)<0\hfill \\ 0\hfill & \text{if }(𝐯_i𝐯_j)(𝐫_i𝐫_j)0\hfill \end{array}.$$
(3)
The terms $`h_{ij}`$ and $`c_{ij}`$ are, respectively, the average smoothing length and sound speed associated with particles $`i`$ and $`j`$. Here $`f_i`$ is the so-called form function for particle $`i`$, defined by

$$f_i=\frac{|𝐯|_i}{|𝐯|_i+|\times 𝐯|_i+\eta ^{}c_i/h_i},$$
(4)
where

$$(𝐯)_i=\frac{1}{\rho _i}\underset{j}{}m_j(𝐯_j𝐯_i)_iW_{ij}$$
(5)
and

$$(\times 𝐯)_i=\frac{1}{\rho _i}\underset{j}{}m_j(𝐯_i𝐯_j)\times _iW_{ij}.$$
(6)
The function $`f_i`$ acts as a switch, approaching unity in regions of strong compression ($`|𝐯|_i>>|\times 𝐯|_i`$) and vanishing in regions of large vorticity ($`|\times 𝐯|_i>>|𝐯|_i`$). Consequently, this AV has the advantage that it is suppressed in shear layers. In this paper we use $`\eta ^2=0.01`$, $`\eta ^{}=10^5`$ and $`\alpha =\beta =\gamma /2`$. This choice of AV treats shocks well, while introducing only relatively small amounts of numerical viscosity (Lombardi et al., 1999).

The code was parallelized using MPI (the message passing interface), to run efficiently on multiple processors. Simpler methods, such as those in HPF, prove impossible to use for SPH loops, since the “gather-scatter” technique, which enforces Newton’s Second Law, disallows the use of DO INDEPENDENT loops. In our code, the lists of neighbors are split and divided up evenly between the processors. All loops over SPH particles containing sums over quantities involving neighbors, including the AV, are distributed among the processors, summed in part, and collected at the end of the relevant subroutine by an MPI\_ALLGATHER command. The gravitational force calculations have been parallelized using the rfftwnd\_mpi package of the FFTW library (Frigo & Johnson, 1997). All 3D grids used by gravity routines are split up and distributed evenly among processors in the $`z`$-direction, which saves a considerable amount of computer memory. Force contributions from grid cells contained within each process are calculated in parallel, and gathered at the end of the subroutine by an MPI\_ALLGATHER command. In benchmarking tests, we have found that our parallel code scales very well when using up to 32 processors on a distributed shared-memory supercomputer. For example, doubling the number of processors up to 32 results in an increase of total CPU time by no more than $`10\%`$ on the SGI/Cray Origin2000 supercomputer. Both the hydrodynamics and gravitational subroutines show similar speed-up when increasing the number of processors.

Our simulations employ equal mass SPH particles, both to keep the resolution high in the stellar cores and to minimize spurious mixing during the simulation; see Lombardi et al. (1999) for a discussion of spurious transport induced by unequal mass particles. Due to extreme central densities in the parent stars, the Courant stability condition requires exceedingly small timesteps. Consequently, an elapsed physical time of one hour in our high resolution calculations requires roughly 1000 iterations and 100 CPU hours on the SGI/Cray Origin2000 supercomputer at NCSA. The three high resolution calculations presented in this paper each required a total of $`1000`$ CPU hours.

### 2.2 Initial Conditions

All of our hydrodynamic simulations use parent star models whose characteristics are based on the results of YREC calculations, as discussed in Sills & Lombardi (1997). In particular, we evolved (non-rotating) main-sequence stars of total mass $`M=0.6`$ and $`0.8M_{}`$ with a primordial helium abundance $`Y=0.25`$ and metallicity $`Z=0.001`$ for 15 Gyr, the amount of time needed to exhaust the hydrogen in the center of the $`0.8M_{}`$ star. The thermodynamic and structural profiles of the parent stars are shown in Figure 1. The total helium mass fractions for the $`0.6`$ and $`0.8M_{}`$ parent stars are 0.286 and 0.395, and their radii are $`0.52R_{}`$ and $`0.96R_{}`$, respectively. From the pressure and density profiles of these models, we compute the entropy profile and assign values of $`A`$ to SPH particles accordingly. In addition, the chemical abundance profiles of 15 different elements are used to set the composition of the SPH particles. To minimize numerical noise, each parent star’s SPH model is relaxed to equilibrium using an artificial drag force on the particles, and then these relaxed models are used to initiate the collision calculations.

The stars are initially non-rotating and separated by 5 $`R_{TO}`$, where $`R_{TO}=0.96R_{}`$ is the radius of a turnoff star. The initial velocities are calculated by approximating the stars as point masses on an orbit with zero orbital energy (i.e., we assume the parabolic limit, appropriate for globular clusters) and a pericenter separation $`r_p`$. A Cartesian coordinate system is chosen such that these hypothetical point masses of mass $`M_1`$ and $`M_2`$ would reach pericenter at positions $`x_i=(1)^i(1M_i/(M_1+M_2))r_p`$, $`y_i=z_i=0`$, where $`i=1,2`$ and $`i=1`$ refers to the more massive star. The orbital plane is chosen to be $`z=0`$. With these choices, the center of mass resides at the the origin.

### 2.3 Hydrodynamic Results and Conversion to YREC Format

Table 1 summarizes the parameters and results of all of our SPH simulations. The first column gives the name by which the calculation is referred to in this paper; we use lower case letters to distinguish the present calculations from the corresponding calculations involving polytropic parent stars presented in Lombardi, Rasio & Shapiro (1996). The second and third columns give the masses $`M_1`$ and $`M_2`$ of the parent stars. Column (4) gives the ratio $`r_p/(R_1+R_2)`$, where $`r_p`$ is the pericenter separation for the initial orbit and $`R_1+R_2`$ is the sum of the two (unperturbed) stellar radii. This ratio has the value $`0`$ for a head-on collision, and $`1`$ for a grazing encounter. Column (5) gives the number of SPH particles. Column (6) gives the final time $`t_f`$ at which the calculation was terminated. Column (7) gives the total angular momentum of the merger remnant. Column (8) gives the ratio $`T/|W|`$ of rotational kinetic energy to gravitational binding energy of the (bound) merger remnant in its center-of-mass frame at time $`t_f`$. Column (9) gives the mass of the merger remnant (the gravitationally bound fluid). Columns (10) and (11) give the velocity components $`V_x`$ and $`V_y`$ for the merger remnant’s center of mass at time $`t_f`$ in the system’s center-of-mass frame. Since the amount of mass ejected during a parabolic collision is small, the merger remnant never acquires a large recoil velocity. The case e, f and k simulations implemented $`N_N=100`$ neighbors per particle, while case j implemented $`N_N=32`$ neighbors. Case j is a low resolution calculation of a nearly head-on collision.

Figure 2 illustrates the dynamical evolution for case e: a turnoff main-sequence star ($`M_1=0.8M_{}`$) collides with a slightly less massive star ($`M_2=0.6M_{}`$). The parabolic trajectory has a pericenter separation $`r_p=0.25(R_1+R_2)`$. The first collision at time $`t2`$ hours disrupts the outer layers of the two stars, but leaves their inner cores essentially undisturbed. The two components withdraw to apocenter at $`t3.2`$ hours before colliding for the second, and final, time. The merger remnant is rapidly and differentially rotating, and shear makes the bound fluid quickly approach axisymmetry within a few dynamical timescales.

Figure 3 shows contours of the density $`\rho `$, the entropy $`A`$, and the $`z`$component of the specific angular momentum $`j=\omega r_{cyl}^2`$ for each collision product in a slice containing its rotation axis. Here $`r_{cyl}`$ is the cylindrical radius measured from the rotation axis. The contours shown represent an average over the azimuthal angle, as the product approaches axisymmetry at the end of the simulation. By definition, the collision products are not barotropes since the specific angular momentum $`j`$ (and hence $`\omega `$) has a distinct dependence on $`z`$. For a rotating star, stable thermal equilibrium requires $`d\omega /dz=0`$ in chemically homogeneous regions, as discussed in §3.2. From the set of specific angular momentum contours presented in Figure 3, it is therefore evident that the outer layers of the collision products (where the composition gradient is small) are not in thermal equilibrium. The shape of the specific angular momentum contours is due to the tendency of dynamical shear to equalize the angular velocity $`\omega `$ on surfaces of constant density (Pinsonneault et al., 1989), so that $`j`$ increases with $`z`$ at fixed $`r_{cyl}`$.

To convert our three-dimensional hydrodynamic results into the one-dimensional YREC format, we begin by averaging the entropy and specific angular momentum values in $`25`$ bins in enclosed mass fraction. The resulting profiles are extrapolated over the outermost $``$5% of the remnant mass, since in practice the simulations must be terminated before all bound fluid has fallen back to the remnant’s surface (for a discussion of how the values of the final time $`t_f`$ are chosen, see §2.3 of Lombardi, Rasio & Shapiro, 1996). The specific angular momentum profile is then renormalized (by a factor of close to unity), to ensure that the extrapolation does not modify the total angular momentum of the remnant. With the $`A`$ and $`j`$ profiles given, the structure of the rotating remnant is uniquely determined in the formalism of Endal & Sofia (1976), by integrating the general form of the equation of hydrostatic equilibrium \[see their eq. (9)\]. To do so, we implement an iterative procedure in which initial guesses at the central pressure and angular velocity are refined until a self-consistent YREC model is converged upon. The temperature profile is calculated from the ideal gas equation of state, while the luminosity profile is calculated according to the equation of radiative transport. This first approximation for the structure of the star was used as a starting model for the stellar evolution code.

Although the pressure, density, and angular velocity profiles of the resulting YREC model cannot also be simultaneously constrained to equal exactly the corresponding (averaged) SPH profiles, the differences are slight (see, for example, the $`\omega `$ profiles in Figure 4). Indeed, the subsequent evolutionary tracks are not significantly different for remnant models that use the form of the $`\omega `$ profile taken directly from SPH. The decision to constrain the $`j`$ profile, as opposed to the angular velocity $`\omega `$ profile, guarantees that the total angular momentum of the initial YREC model agrees with that from the SPH results. Furthermore, the SPH $`\omega `$ profile can be slightly less accurate than the $`j`$ profile, since the latter is less sensitive to the resolution of the calculation.

## 3 EVOLUTION OF THE ROTATING COLLISION PRODUCTS

### 3.1 Implementation of the Stellar Evolution Code

Our stellar evolution calculations are performed with YREC. YREC is a one-dimensional code, in which the star is divided into shells along surfaces of constant gravitational plus rotational potential. The code solves the equations of stellar structure with the Henyey technique, and follows the rotational evolution with the formalism of Endal & Sofia (1976). Guenther et al. (1992) give a detailed description of the physics implemented in the evolution code. We used the same opacities, equation of state and model atmospheres as in Paper I. The free parameters in the code (the mixing length and parameters that set the efficiency of angular momentum transport and rotational chemical mixing) are set by calibrating a solar mass and solar metallicity model to the sun.

Rotation is treated by evaluating physical quantities on equipotential surfaces rather than the spherical surfaces usually used in stellar models. The hydrostatic equilibrium and radiative transport equations contain terms that account for the lack of spherical symmetry in the rotating star. A number of rotational instabilities that transport angular momentum and material within the star are followed, including dynamical shear (Pinsonneault et al., 1989), meridional circulation (von Zeipel, 1924), secular shear (Zahn, 1974), and the Goldreich-Schubert-Fricke (GSF) instability (Goldreich & Schubert, 1967; Fricke, 1968). Angular momentum transport and the associated chemical mixing are treated as diffusion processes, with diffusion coefficients that account for each active mechanism within unstable regions of the star. The diffusion coefficients are proportional to circulation velocities, which have been estimated by Endal & Sofia (1987). In addition to the internal rearrangement of angular momentum, angular momentum also can be drained from an outer convection zone through a magnetic wind, using the formalism given by Chaboyer, Demarque & Pinsonneault (1995). The Endal & Sofia scheme is valid across a wide range of rotation rates, for a restricted class of angular momentum distributions. This scheme requires that the angular velocity is constant on equipotential surfaces, which does not allow for modeling of latitude-dependent angular velocity profiles. See Meynet & Maeder (1997) for a detailed discussion of the validity of this approach. For a detailed description of the implementation of rotation in YREC, see Pinsonneault et al. (1989).

### 3.2 Evolutionary Results

Although in dynamical equilibrium, the nascent collision products are initially far from thermal equilibrium and contract on a thermal timescale. During this stage, their luminosity is predominately provided by the conversion of gravitational potential energy, analogous to the situation in pre-main-sequence stars. After the contraction, central hydrogen burning quickly becomes significant, and we refer to this point as the zero age main sequence (ZAMS). Note, however, that unlike normal ZAMS stars, these stars are far from being chemically homogeneous. Many qualitative features of the subsequent evolution are similar to that of normal intermediate mass stars; for example, it is possible to define the subgiant and giant stages of evolution. Nevertheless, the timescales for this evolution, and even the actual evolutionary tracks, can be substantially different from those of normal stars, because of the initial structure of the collision product and because of the effects of the high rotation rate.

The initial total angular momentum of the collision products (Table 1) is as much as 10 times larger than that of normal stars of comparable mass at the pre-main-sequence birth line. The stellar collision products differ substantially from normal pre-main-sequence stars in another significant way. Normal pre-main-sequence stars have large convective envelopes, and therefore can lose substantial amounts of angular momentum via a magnetic wind. The sun has lost about 99% of its initial angular momentum, for example (Pinsonneault et al., 1989). However, the stellar collision products are not convective at the end of the SPH simulation, nor do they develop convection zones during the thermal readjustment and main-sequence stages. Therefore, the collision products do not lose any angular momentum by a magnetic wind until the subgiant branch.

Rotational instabilities transport angular momentum from the center of the star to the surface approximately on an Eddington-Sweet timescale (the thermal timescale divided by the ratio of centrifugal acceleration to gravitational acceleration). For the rapidly rotating stellar collision products, this timescale is on the order of only $`10^7`$ years for the stars as a whole, and is even shorter in the outer layers. The dominant mixing mechanism in these stars is the GSF instability, which addresses axisymmetric perturbations: A structure in which the angular velocity $`\omega `$ is constant on cylinders is the only thermally stable configuration in regions of homogeneous chemical composition. Since dynamical shear acts to remove any latitude dependence of the angular velocity (see Pinsonneault et al., 1989), regions that are subject to the GSF instability evolve towards rigid rotation. In the case of our stellar collision products, the angular velocity profile is initially a steeply decreasing function of enclosed mass fraction (Figure 4). Therefore, the GSF instability serves to raise the surface rotation rate and decrease the interior rotation rate, flattening the $`\omega `$ profile. Since the timescale for this to occur is shortest at the surface, the $`\omega `$ profile begins to flatten there, and then the flattened profile penetrates inwards. The star has a very large reservoir of angular momentum, and this redistribution of angular momentum can cause the surface of the star to rotate faster than its breakup velocity.

An obvious solution here would be simply to remove the mass from the calculation as it becomes unbound to the star. Unfortunately, under this scheme, rapidly rotating stellar collision products could be completely disrupted. The total angular momentum of the star is so large that if the star were rotating as a solid body, its rotational velocity would be more than its breakup velocity. This is true immediately following the collision, when the star has its largest radius, and is even more obvious as the star collapses to the main sequence. The breakup velocity is proportional to $`R^{1/2}`$, so it increases as the radius $`R`$ decreases. However, the surface rotation rate is proportional to $`R^1`$, and increases faster than the breakup velocity as $`R`$ decreases. As the GSF instability works to flatten the rotation profile further and further into the center of the star, more of the star becomes unbound, and eventually the entire star can be spun apart. This poses a problem for the theory that blue stragglers are created via stellar collisions: either blue stragglers are not created through physical off-axis stellar collisions, or some mechanism(s) can remove angular momentum from the star on short timescales. Possible mechanisms will be discussed in §4.

For now we will simply postulate that the angular momentum of the collision products can be reduced. In order to study the effects of rotation on stellar collision products, we performed some experiments on case k ($`0.6M_{\mathrm{}}+0.6M_{\mathrm{}},r_p=0.25(R_1+R_2)`$). We artificially decreased the initial rotation velocity of each shell by a large factor before evolving the model. In particular, we divided the rotation velocities by 1000, 100, 10, 5 and 2, which reduced their total angular momentum by almost the same factors. In this last case, the star still had enough angular momentum that it began to rotate faster than its breakup velocity. The evolutionary tracks for the other angular momentum values are shown in Figure 5. The tracks with the least angular momentum are not noticeably different from the evolutionary track of the corresponding non-rotating model. As the total angular momentum of the star increases, the star spends more time on the main sequence due to two effects. First, rotation provides an additional source of support for the star, so that it does not have to produce as much pressure through nuclear burning to prevent its collapse. Second, rotationally induced mixing moves hydrogen into the core of the star, increasing the amount of available fuel.

We also evolved the product of the low resolution SPH simulation of a barely off-axis collision (case j). In this situation, two 0.6 $`M_{\mathrm{}}`$ parent stars collide with an impact parameter of only $`r_p=0.02R`$, where $`R=0.52R_{}`$. The SPH results are used directly as the initial conditions for the evolution code as described in §2.3, without rescaling the angular velocity profile. This collision has sufficiently little angular momentum that even if the collision product rotated as a solid body, it would rotate slower than its breakup velocity. We used this case to confirm that the results of SPH simulations could indeed be evolved directly, and to see what the effect of even a barely off-axis collision was on the subsequent evolution of the collision product. The results are shown in Figure 6. For comparison, a non-rotating version of the same star is also shown. The non-rotating star was evolved using the same starting model as the rotating case, but with zero angular momentum. The effects of rotation, even for a collision that is barely off-axis, are remarkable. This collision product has an initial total angular momentum of $`2.1\times 10^{50}`$ g cm<sup>2</sup> s<sup>-1</sup>, comparable to that of a pre-main-sequence star of the same mass. Since the collision product does not lose angular momentum during its evolution through the main-sequence stage, its rotational velocity remains quite large (272 km s<sup>-1</sup> at the turnoff). Substantial rotational mixing occurs, prolonging its main-sequence lifetime to 3.7 Gyr (compared to 1.8 Gyr for the equivalent non-rotating star).

The evolutionary tracks for the case e, f and k collision products are shown in Figure 7. Here we have reduced the initial angular velocities by a factor of 5, in order to prevent the stars from rotating faster than breakup. None of the other profiles are changed initially, although YREC converges on a solution which satisfies the rotational equations of stellar structure. Some details of the various evolutionary tracks are given in Table 2. Cases e and f involve the same parent stars (0.6 $`M_{}`$ and 0.8 $`M_{}`$) but have different pericenter separations, namely 0.25 and 0.5 $`(R_1+R_2)`$, respectively. Since case f has significantly more angular momentum than case e, much more mixing occurs, and the evolutionary tracks are substantially different. The track of the case f collision product also shows the most striking differences from the track of a normal, slowly rotating star. Although not even the case f collision product is fully mixed, hydrogen is carried into the core nearly as fast as it is depleted by nuclear burning, extending the main-sequence lifetime in case f by a factor of 5 over that of case e.

In addition to transporting angular momentum and carrying hydrogen into the core, rotational instabilities also carry helium to the surface, affecting the surface opacity and causing the star to become continuously brighter and bluer while on the main sequence. The most dramatic example of this effect can be seem in the case f track of Figure 7, which reaches a luminosity of $`300L_{}`$ and a surface temperature of $`2.5\times 10^4`$ K while still on the main sequence. In this case, the helium mass fraction at the surface reaches an incredible $`Y=0.85`$. Figure 8 shows the surface abundance of helium as a function of time for the collision products featured in this paper.

The differences between the case e and case f tracks, which correspond to collision products of nearly the same mass, help to emphasize that the total mass of a collision product is not the primary factor in determining its relative position on a Hertzsprung-Russell (HR) diagram, unlike the situation for normal main-sequence stars. As a second example, consider the collision product of case k, which has a mass of only $`1.16M_{}`$, but which is nevertheless brighter and bluer at turnoff than the $`1.36M_{}`$ product of case e. As discussed in Paper I and Sills & Lombardi (1997), a collision that involves a turnoff mass star (such as in cases e and f) will produce a blue straggler with a hydrogen-depleted core, since the core of the turnoff star sinks to the center of the collision product. Therefore, the main-sequence lifetime for the higher mass products is generally shorter. In lower mass collisional blue stragglers, there is then more time for helium to be mixed to the surface, and the surface helium abundance can grow larger than in the higher mass collision products (compare cases e and k in Figure 8). As a result, a lower mass collision product can become brighter and bluer than one of higher mass. Consequently, blue straggler masses that have been determined from color-magnitude diagram position could be in error if the blue stragglers are rapidly rotating.

Figure 9 shows the surface rotational velocity as a function of time for the four cases featured in this paper. The ages of the ZAMS, the turnoff, and the base of the giant branch are marked in each of the panels. All four cases show similar behavior: the rotational velocity adjusts rapidly as the star contracts to the main sequence. Once on the main sequence, the surface rotational velocity remains approximately constant. As the star becomes a subgiant, it develops a surface convection zone, and the standard angular momentum loss mechanism (via a magnetic wind) begins to operate. The loss rate is proportional to the surface angular velocity to the third power, so these rapidly rotating stars very efficiently shed most of their angular momentum before reaching the giant branch.

## 4 SUMMARY AND DISCUSSION

We have modeled blue stragglers created in off-axis collisions between cluster main-sequence stars. Since even a minutely off-axis collision can produce a rapidly rotating star, most collision products cannot remain gravitationally bound as they contract to the main sequence unless they somehow lose angular momentum. The standard angular momentum loss mechanism (via a magnetic wind) does not operate in collisional blue stragglers, since magnetic braking is effective only in stars with surface convection zones deep enough to support a dynamo. The stellar collision products lack a surface convection zone, even during their collapse to the main sequence. Therefore, either blue stragglers cannot be created by stellar collisions, or some other mechanism for removing angular momentum must occur.

Angular momentum can be lost from stars in a variety of ways. Stars more massive than these collision products have substantial mass loss rates during their main-sequence evolution, which can be on the order of $`10^9M_{\mathrm{}}`$ yr<sup>-1</sup> for stars of 3 $`M_{\mathrm{}}`$ (Lamers, 1981). Such a high mass loss rate will remove a substantial amount of angular momentum from a star. The mass loss in these stars takes the form of a line radiation-driven wind (Castor, Abbott & Klein, 1975), which requires high luminosities in order to remove a large amount of mass. This mechanism may be of importance in extremely bright collision products (such as in case f), although low luminosity blue stragglers are unlikely to be significantly affected.

If the collision product is a member of a binary system, it can lose angular momentum through the orbital evolution of the system. Angular momentum loss is usually discussed in the context of contact or near-contact binary systems, and requires that one of the stars has a stellar wind and a co-rotating magnetic field (Stepien, 1995). It is possible to have a stellar collision in which the product has a binary companion (through an encounter of a single star with a binary system (Davies, 1995), or two binary systems with each other (Bacon, Sigurdsson, & Davies, 1996). However, it is not clear that the resulting binary system will be in contact, and it is highly unlikely that all blue stragglers formed in such a fashion.

The most plausible mechanism for removing angular momentum from these stellar collision products is disk-locking. In low mass pre-main-sequence stars, the torque placed on a star through a magnetic field coupled to a disk of material at a few stellar radii is enough to slow the rotation of these stars by as much as an order of magnitude (Sills, Pinsonneault & Terndrup, 2000). This mechanism requires the presence of both a magnetic field and a disk of material around the star. The SPH simulations of the stellar collisions show that typically a few percent of the mass of the system is not bound to the final collision product. However, this mass will leave the system on a dynamical time (on the order of hours) and cannot become the disk. A more reasonable source of disk mass is the outer layers of the star that are spun off during the initial contraction on a thermal timescale. Since the amount of mass necessary for such a disk is on the order of a few hundredths of a solar mass, the disk-locking could begin very soon after the collision, and remove the angular momentum quite quickly. A simple calculation shows that this scenario is not unreasonable. We assume that angular momentum is carried away by mass which is coupled to the surface rotation rate of the star down to the Alfven radius (Königl, 1991). In order to lose 80% of the star’s angular momentum (in case f), we need to lose only $`0.1M_{\mathrm{}}`$. The other ingredient in this scenario is the presence of a magnetic field. We know that low mass stars have substantial magnetic fields. It is not clear what happens to those fields after the stellar collision, and whether they will be coherent enough to lock the star to a disk. More work is necessary to determine if this or any other angular momentum loss mechanism occurs in stellar collision products.

We assumed that the stellar collision products would lose a large fraction (up to $``$80%) of their initial angular momentum early in their evolution. We then continued their evolution through the main sequence to the base of the giant branch. The rotating stars are bluer and brighter than their non-rotating counterparts, and have longer main-sequence lifetimes. Both of these changes to the evolutionary tracks will have an effect on the derived blue straggler distributions in the color-magnitude diagram, and therefore will change our conclusions about the dynamical state of globular clusters.

Since rotating blue stragglers live on the main sequence longer, we expect that their distributions will be more populated near the main sequence than those presented in Sills & Bailyn (1999). Since the observed distributions of blue stragglers in the color-magnitude diagram are already extended to cooler temperatures than the predicted distributions, the increase in blue straggler main-sequence lifetime will exacerbate the lack of agreement between observations and theory. It is possible that the width of the observed distribution is a selection effect, caused by the blended light of chance superpositions or true binary stars. More careful observations of blue straggler populations in globular clusters may clear up this particular mystery.

While none of the stars investigated in this work are fully mixed, rotation does cause helium to be mixed to the surface. Rotating blue stragglers tend to be much brighter than their non-rotating counterparts, and much hotter. Therefore, rotation could be responsible for anomalously bright blue stragglers, such as the central 6 blue stragglers in NGC 6397 (Lauzeral et al., 1992) or the brightest three blue stragglers in M3 (Ferraro et al., 1997). Detailed comparisons between observed distributions and theoretical distributions based on rotational tracks are necessary to determine if there are two distinct populations of collisional blue stragglers in globular clusters, or if all blue stragglers can be explained by a population with a range of collisional impact parameters.

The only blue straggler in a globular cluster to have its rotation rate measured is BSS 19 in 47 Tucanae (Shara, Saffer, & Livio, 1997). This star is rotating at $`v\mathrm{sin}i=155\pm 55`$ km s<sup>-1</sup>, which is consistent with the rotation rates of the collision products presented in this work. The position of this star in the HR diagram is shown in Figure 7. The rotation velocities of the collision products in the vicinity of this star range from about 230 km s<sup>-1</sup> to 400 km s<sup>-1</sup>, which is consistent with the BSS 19 data point and a reasonable inclination angle $`i`$ between $`20^{}`$ and $`40^{}`$. Observations of more rotation rates of blue stragglers in globular clusters are necessary to resolve some of the issues involved in studying off-axis collision products. However, we can say that the conclusion of Shara, Saffer, & Livio (1997) that such a rapid rotation rate necessarily indicates a binary merger formation for BSS 19, rather than a collisional origin, is currently unwarranted. This conclusion was based on the assumption of Leonard & Livio (1995) that collisional blue stragglers would develop convective envelopes shortly after the collision, and would then lose angular momentum through a magnetic wind. As shown in Paper I and confirmed in this work, collision products do not develop significant surface convection zones, and therefore do not suffer from standard magnetic braking.

The other blue stragglers with measured rotation rates are found in the old open cluster M67 (Peterson, Carney & Latham, 1984). These stars are rotating more slowly than the blue straggler in 47 Tuc, and much more slowly than the models presented here. A large fraction of these blue stragglers are in binary systems. Could these two circumstances be related? Blue stragglers are thought to be formed through mass transfer in a binary system as well as through direct stellar collisions. Detailed calculations will be necessary to determine if binary merger products have different rotational properties than stellar collision products, and if we can distinguish between the two blue straggler creation mechanisms on the basis of their current rotation rates.

This research was supported in part by NSF Grants AST-9618116 and PHY-0070918, and NASA ATP Grant NAG5-8460 at MIT, and by an award from Research Corporation and NSF Grant AST-0071165 at Vassar. This work was also supported by the National Computational Science Alliance under grant AST980014N and utilized the NCSA SGI/Cray Origin2000. F.A.R. was supported in part by an Alfred P. Sloan Research Fellowship. A.S. gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada.


text:
# Radio Jet-blown Neutral Hydrogen Supershells in Spiral Galaxies?

## 1 Introduction

An important manifestation of the activity inside the disks of gas-rich galaxies is their highly structured HI distribution, marked by cavities, shells and supershells. First discovered in the Milky Way (Heiles 1979, 1984), such features are now known to exist in a number of spiral galaxies (e.g. Lehnert & Heckman 1996, Irwin & Seaquist 1990, Puche et al. 1992, Brinks & Bajaja 1986). Exceptionally huge HI arcs and loops extending across several kiloparsecs have been identified with greater clarity in the HI images of a number of edge-on spirals, such as NGC 5775 (Irwin 1994), NGC 4631 (Rand & van der Hulst 1993), NGC 3044 (Lee & Irwin 1997, hereafter LI97) and NGC 3556 (M 108, King & Irwin 1997, hereafter KI97). These have been interpreted as expanding supershells because of a loop-like or circular appearance in projection and either a persistence over a wide velocity range or, in a few cases, as some evidence for expansion in position-velocity space. Two main classes of explanations for the supershells posit the source of their kinetic energy to be, respectively, internal and external to the parent galaxy.

The internal source model involves starbursts, driving stellar winds (or superwinds) and subsequent supernova explosions (e.g. Lehnert & Heckman 1996). The “chimney model” (Norman & Ikeuchi 1989), for example, attempts to explain disk-halo features and other halo gas via processes related to underlying star formation. The association between extraplanar H$`\alpha `$ filaments and star forming regions in the disk of NGC 891 and other correlations between halo emission and in-disk tracers of star formation (Dahlem et al. 1995; Rand 1997) argue in favour of such models. If the presence of HI supershells is found to correlate with the existence of other halo gas, as might be expected in the chimney model, then stellar winds and supernovae are expected to be responsible for the HI supershells as well.

The main difficulty with the starburst model for HI supershells lies in the required input energies for the largest shells. Using standard assumptions that the expanding supershells are in the post-Sedov phase following an ‘instantaneous’ injection of energy (cf. Chevalier 1974), HI supershells often require energy input from staggering numbers of spatially correlated supernova events. This was realized early on for our own Galaxy (Heiles 1979, 1984). For external edge-on galaxies, since we are selectively observing only the largest shells, the energy deficit problem is exacerbated. In some cases, hundreds of thousands of clustered supernovae are required (e.g. KI97, LI97), a conclusion which is not changed significantly if the energy is injected continuously over the lifetime of the shells. Other evidence against star formation processes creating the HI shells is also emerging. Rhode et al. (1999) find no optical evidence for recent star formation in the numerous lower energy HI holes of Holmberg II and note that X-ray and FUV emission are also absent. They conclude that supernovae have not played a part in the formation of the HI shells. Efremov et al. (1998) outline numerous other examples in which there appears to be no relation between HI shells and star formation. They, as well as Loeb & Perna (1998), propose that the HI shells are produced, instead, by gamma ray bursts.

The alternative external source hypothesis invokes infall of massive gas clouds on to the galactic plane, as a result of gravitational interaction with neighbouring galaxies (see Tenorio-Tagle & Bodenheimer 1988). This resolves the energy problem since input energy is then a function of the mass and velocity of the infalling cloud. Evidence in favour of this hypothesis comes from observations of high velocity clouds (HVCs) around our own Milky Way and the signatures of interaction in M 101 (van der Hulst & Sancisi 1988) and NGC 4631 (Rand & Stone 1996). It does, however, require that the galaxy be in some way interacting with a companion or, at least, that sufficiently massive clouds be in the vicinity.

Recent observations are revealing galaxies which are apparently isolated, yet harbour extremely large HI supershells. Two striking examples are the nearby, SB(s)cd galaxy, NGC 3556 (KI97) and the SBc galaxy, NGC 3044 (LI97). Both of these galaxies exhibit radio continuum halos extending to $`5`$ kpc from the galactic plane and have a number of supershells requiring energies up to a few $`\times `$ 1056 ergs. These supershells are too large and energetic to have been produced by conventional clustered supernovae. At the same time, there appears to be no evidence for interaction or nearby companions, either.

We propose here a new explanation for HI supershells. That is, that they have been formed by radio jets which plow through the interstellar medium (ISM), accreting ISM gas and sometimes inflating bubbles. This allows for an internal energy source for the HI shells, provides a natural explanation for any spatial symmetries seen in the HI features, and also resolves the energy problem. In Sect. 2, we provide arguments in favour of jet inflated HI bubbles, Sect. 3 presents the model, and Sect. 4 discusses the implications of this scenario.

## 2 The Case for a Radio Jet Origin for HI Supershells

### 2.1 Radio Jets and AGN in Spiral Galaxies

Seyferts are one class of disk galaxy for which several examples of the nucleus ejecting a radio jet pair have been found (e.g. Ulvestad & Wilson 1984a, 1984b, Kukula et al. 1995, Aoki et al. 1999). Likewise, several cases of jets occurring in normal spiral galaxies have been reported (e.g. Hummel et al. 1983). Prominent examples include NGC 3079 (de Bruyn 1977), NGC 5548 (Ulvestad et al. 1999) and Circinus (Elmouttie et al. 1998).

The total energy output from such nuclear activity can approach $`10^{60}`$ erg, assuming that the nuclear activity lasts $`10^7`$ years, and arises from accretion onto a supermassive black hole of $`10^8M_{}`$ at 10% of the Eddington limit with the canonical 10% efficiency. In “normal” spiral galaxies, the central mass may be lower; for example, the bipolar outflow from the SAB(rs)cd galaxy, M 101 may be produced by a central $``$ 106 M black hole (Moody et al. 1995). While jets are not always observed directly, the growing number of supermassive black holes inferred to be present in the nuclei of normal disk galaxies (Kormendy & Richstone 1995, van der Marel 1999), including late-type disks (Ho et al. 1997), the Milky Way itself (e.g. Genzel et al. 1997) and gas-rich, large low-surface-brightness galaxies (Schombert 1998), lend weight to the idea that many such galaxies may have undergone phases of nuclear activity accompanied with a bi-directional ejection of relativistic plasma jets, before entering a dormant nuclear phase.

Recent studies have also revealed that the ejection of jets can take place at small angles to the large-scale galactic disk (e.g. Nagar & Wilson 1999; Kinney et al 2000), plausibly leading to clear signatures of their dynamical interaction with the ISM. A classic case of such an interaction is the bow-shaped morphology of the radio lobe in the well known Seyfert galaxy, NGC 1068, which strongly suggests that the kiloparsec-scale radio jets are ploughing through the disk medium (Wilson & Ulvestad 1987, Axon et al. 1998). A dynamical interaction of the jets with the disk is also evident in the case of the spiral galaxy M 51 which is at the low end of nuclear activity scale. Here, the jets seem to have created a pair of ‘radio lobes’ on opposite sides of the nucleus; whereas one of the lobes is bow-shaped, the other one is actually resolved into a ring which is also detected in H$`\alpha `$ emission (Ford et al. 1985). Another example is NGC 4258 whose VLBI jet has been ejected close to the galaxy disk (see Cecil et al. 1995, Herrnstein et al. 1997, Cecil et al. 1998). Large scale (14 kpc length) “anomalous arms” are also seen within the galactic disk and have been interpreted as manifestations of a larger-scale jet (e.g. Martin et al. 1989, but see Cox & Downes 1996). In M 101, the bipolar outflow is also roughly confined to the disk, with the outflow not expected to extend beyond a height of 400 pc (Moody et al. 1995). Additional support for the jet-disk interaction hypothesis in these disk galaxies comes from the detection of shock-excited optical emission lines associated with the radio lobes, as discussed in some of the references cited above.

### 2.2 Jet Inflated Bubbles in the Lobes of Radio Galaxies

It is further interesting to note that radio bubbles and shells of synchrotron plasma have been discovered within the lobes of a few radio galaxies. Two spectacular examples are 3C310 (van Breugel & Fomalont, 1984) and Hercules A (Dreher & Feigelson 1984). In Hercules A, in addition, peculiar dark shells, about 3 kpc in size, have been found straddling the nucleus along the radio outflow axes (Baum et al. 1996). For the dark shells, Baum et al. prefer an explanation in terms of expanding bubbles of hot gas ejected from the active nucleus along the radio axis, in agreement with plasmon-like ejection from the core, though they do not rule out other possibilities.

These cases provide a clue that the radio bubbles/shells may have been ‘puffed’ up due to localized instabilities in the radio jets, occurring at distances of a few kiloparsecs from the active nucleus (e.g. Morrison & Sadun 1996). Should such instabilities arise also in the jets ejected within disk galaxies (Sect. 2.1), it is quite plausible that the resulting increase in pressure at those locations would inflate bubbles and shells out of the ambient ISM material which is HI rich, in the case of spirals. Such features, puffed out of the disk during the brief phase of jet instability, would remain visible past the radiative lifetime of the synchrotron plasma against the radiative and expansion losses, which is typically of order $`10^610^7`$ years.

### 2.3 Candidate Galaxies

#### 2.3.1 NGC 3556

In NGC 3556, the two most prominent supershells are located symmetrically near the opposite extremities of the galaxy. KI97 have carried out a detailed mapping of the HI and radio continuum emission from this nearly edge-on galaxy located 11.6 Mpc away. The superimposition of 2 channels of HI emission, symmetrically placed in velocity about systemic, are shown in Fig. 1 a. The most obvious features are two loops of HI, producing extensions to the east and west along the major axis. Fig. 1b shows the HI emission from these velocities superimposed on the radio continuum map, the latter taken from entirely independent observations. The abrupt density drop-off in the radio continuum map occurs roughly where the optical galaxy ends (Fig. 1a inset). The salient features of this system most relevant to the present study (see Fig. 1; KI97) are:

(a) Two symmetrically placed giant HI loops are seen originating at a projected galactocentric radius of about 12 kpc on the east and west major axis. They appear to originate right at the opposite edges of the optical disk and extend, not perpendicular to the disk plane, but rather parallel to it and outwards.

(b) The giant HI loops or supershells are slightly bent into a Z-symmetry.

(c) These HI loops are most obvious in velocity channels which are equally spaced about the systemic velocity of the galaxy (though they extend over a larger velocity range, see KI97).

(d) From their velocity dependent morphology, each of the loops shows evidence for shell expansion, but with only half of the shell present. In both cases, the receding side (with respect to galaxy rotation) is open.

(e) Other smaller HI loops and extensions do exist at smaller galactocentric radii.

(f) The parameters of these two HI supershells are as follows: The eastern supershell has a mass, M = 5.9 $`\times `$ 107 M, a radius, R = 3.2 kpc, an expansion velocity, $`\mathrm{\Delta }v`$ = 51.7 km s-1, a kinematic age, $`\tau `$ = 6.0 $`\times `$ 107 yr, a kinetic energy, Ekin = 1.6 $`\times `$ 1054 ergs, and an implied input energy (assuming instantaneous input from supernovae) of E = 2.6 $`\times `$ 1056 ergs. The parameters of the western supershell are: M = 2.3 $`\times `$ 107 M, R = 1.85 kpc, $`\mathrm{\Delta }v`$ = 41.4 km s-1, $`\tau `$ = 4.4 $`\times `$ 107 yr, Ekin = 3.9 $`\times `$ 1053 ergs, and E = 3.4 $`\times `$ 1055 ergs.

(g) The energy needed to create the eastern supershell alone would require a star cluster populated by $`>1.7\times 10^5`$ OB stars, if supernovae and stellar winds are the drivers.

(h) The brighter and more spectacular eastern supershell has associated HI gas which extends as far as 30 kpc beyond the eastern edge of the optical disk and possibly much farther (see Fig. 1a inset). The huge eastern HI extension appears to open up from one side of the eastern shell and extends fairly straight outwards in the radial, rather than the vertical direction as if formed in a short period of time in comparison to galactic rotation.

(i) The galaxy contains a large nonthermal radio halo with a scale height of $`5`$ kpc above the plane (see also de Bruyn & Hummel 1979, Bloemen et al. 1993).

(j) Radio continuum emission appears to be associated with the HI loops, but is not spatially coincident with them. The radio continuum emission extends farther out than the HI supershells, as if “funnelling” along or through the HI features (KI97; Fig. 1b).

(k) NGC 3556 has no obvious interacting companions. The brightest galaxy in the vicinity is the 16th magnitude, MCG+09-19-001, 25′′ in size, undisturbed in appearance and $``$ 12 to the east of NGC 3556. This is most likely a background source.

#### 2.3.2 Other Examples

NGC 5775, an interacting galaxy, shows symmetrically placed HI features (Lee et al. 2000) in the sense that HI extensions occur on opposite sides of the major axis in which there appears to be an underlying disturbance; these features occur at similar galactocentric radii, in projection (see also Fig. 3 of Irwin 1994 which shows 3 of the extensions). NGC 2613 (Chaves & Irwin 2000) also shows six HI extensions which occur in pairs symmetrically on either side of the major axis; two of the pairs occur near the ends of the optical disk.

Another potentially very interesting example is M 31. Blitz et al. (1999) show that two massive HI clouds exist symmetrically placed on opposite sides of this galaxy. They bear a remarkable resemblance to extragalactic radio jets. The two clouds are redshifted by $`200`$ km s-1 with respect to the systemic velocity of M 31. Conceivably, this could be explained by ram pressure sweeping as M 31 falls through the IGM towards the Milky Way.

## 3 A New Model for the Formation of HI Supershells

### 3.1 Need for a New Mechanism

As mentioned above, the galaxy NGC 3556 puts to a severe test both of the standard explanations for the supershell phenomenon, namely: (i) localized starburst (generating intense stellar winds, followed by multiple supernova explosions), and (ii) infall of massive gas clouds on to the galactic disk. The major difficulties faced by these models, as highlighted in KI97, are:

(a) Localized starbursts: The energy needed to create the eastern supershell alone would require a star cluster populated by $`>1.7\times 10^5`$ OB stars. Even for the recently discovered ‘super star clusters (SSC)’ in some starburst galaxies, not more than $`10^4`$ OB stars are implied (Meurer et al. 1995). There is no indication of an OB association even remotely approaching the level of SSCs at the locations of the supershells in NGC 3556. Since the kinematic ages of the supershells are $``$ 5 $`\times `$ 107 yr, the starburst would have to have occurred within this time frame. Yet starburst durations and OB association ages are typically also of this order, so some evidence of the starburst might be expected to have survived. No such evidence is presently found, however. NGC 3556 does not show a markedly high supernova rate globally (Irwin et al. 1999) and independent low and high resolution radio continuum observations show no evidence for a source of energy at the bases of the supershells (KI97, Irwin et al. 2000).

(b) Infalling clouds: Several arguments have been advanced against this possibility (KI97). Firstly, this galaxy appears isolated, making the source of the putative clouds puzzlesome. Secondly, to impart sufficient kinetic energy for the creation of the eastern supershell, an infalling cloud of mass $`>10^8M_{}`$ would be required. A cloud this massive would easily have been detected on the sensitive HI maps, yet no such cloud or clouds are seen. (An exception would be if the clouds had very narrow line widths, low fractal dimensions and were optically thick, in which case they could have been missed due to low filling factors.) Thirdly, if an encounter with intergalactic clouds has occurred in the past, it would be necessary to postulate that two very massive clouds just happened to fall in parallel to the major axis at opposite ends of the galaxy, roughly at the same time, a scenario which seems implausible. Fourthly, we could instead suggest that a “rain” of HI clouds is infalling (since other high latitude HI structures are also seen in this galaxy), including two which fell in along the major axis at either end. However, in order for the rain to be no longer visible, the infall would have to be completed over a timescale which is less than or equal to the age of the supershells, i.e. a few 107 years. This is unlikely since infall timescales are of order a dynamical timescale which is typically a few 108 yr.

While these arguments have been applied to NGC 3556 alone, the difficulty with internal energy from starbursting applies to a number of galaxies, as outlined in Sect. 1. Infalling clouds are certainly plausible, but run into difficulties for isolated galaxies such as NGC 3556 and NGC 3044. Barring new information coming to light on these galaxies, we therefore argue that a new model should be considered.

### 3.2 Outline of the Model

As mentioned earlier, the duality of the expanding HI supershells in NGC 3556, together with their locations at the opposite edges of the disk marked by steep density gradients, their Z-symmetric deviations from spherical symmetry, as well as their anomalously large energy requirements (Sect. 2.3.1), together lead us to argue that the origin of these two supershells is linked to the jet phenomenon. In Sect. 2 we have noted several manifestations of jet instabilities giving rise to bubble-shaped features, or shells, in the form of optical nebulosities or nonthermal radio emission in normal spirals and radio galaxies. In the present context of HI supershells, we examine the possibility of the supershell being inflated due to localised instabilities in the radio jets as they plough through the ISM of the disk of NGC 3556. The present discussion should only be viewed as a feasibility check based on energy considerations, and not as a detailed quantitative model for the supershells. Even in the case of NGC 3556, the standard mechanisms for the supershells (Sect. 1) may well have contributed at some level.

According to our proposal, the two radio jets ejected close to the plane of the disk, during the active phase of the nucleus, undergo a rapid flaring as they encounter the region of large scale density decline near the outskirts of the galactic disk. Consequently, both synchrotron plasma and the entrained thermal plasma deposited by the jets in these two regions get heated and the resulting high-pressure bubble of hot plasma undergoes a rapid expansion, sweeping the gas-rich ambient medium into the shape of the HI supershells. An idealization to such a situation is the flaring of a jet crossing an ISM/IGM interface, which was first discussed in the context of radio galaxies analytically by Gopal-Krishna & Wiita (1987), followed by two-dimensional (Norman et al. 1988, Wiita et al. 1990, Zhang et al. 1999) and three-dimensional (Loken et al. 1995, Hooda & Wiita 1996, 1998) numerical simulations. The simulations by Loken et al. showed that intermediate power radio jets associated with wide-angle-tail (WAT) sources undergo a rapid flaring upon crossing the ISM/IGM density discontinuity where a moderately supersonic jet becomes subsonic. As shown by the numerical simulations, a density drop of a factor of just a few can cause such a flaring (e.g. Hooda & Wiita 1996). Observational evidence for the jet flaring comes from Wide-Angle-Tail (WAT) radio sources which are associated with the dominant members of rich clusters of galaxies (e.g. O’Donoghue et al. 1993). Further evidence is provided, e.g., by the recent radio/optical/X-ray study of the WAT radio galaxy 3C465 in the Abell cluster A2634 (Sakelliou & Merrifield 1999). Particularly instructive for the present study is the case of the active disk galaxy IRAS 04210+0400; the two $`10`$ kpc long radio jets associated with this disk galaxy are seen to flare up near the opposite optical boundaries of the galaxy (Holloway et al 1996). It may be noted that during their passage through the inner parts of the galactic disk, the jets are likely to encounter large ISM density fluctuations. However, these are less prone to disrupt the jets because of the higher jet thrust there, combined with the expected small spatial scale of the density fluctuations compared with the jets’ cross-section.

Since the postulated radio jets in NGC 3556 are currently too weak for detection, a point which we re-address below, we shall adopt here the jet parameters inferred for another massive disk galaxy, NGC 4258 (Sect. 2.1), and apply them to known conditions in NGC 3556. The jet velocity, as measured from the emission line gas, is at least vj = 2000 km s-1 which, together with a density of $`\rho _\mathrm{j}`$ = 0.02 $`m_\mathrm{p}`$ g cm-3, where $`m_\mathrm{p}`$ is the mass of the proton, correspond to a kinetic luminosity of Lj $``$ 1042 erg s-1 (Falcke & Biermann 1999). Assuming that the ram pressure in the jet dominates over the internal thermal pressure, then across the shock front we require $`\rho _\mathrm{j}(v_\mathrm{j}v_\mathrm{s})^2=\rho _0v_\mathrm{s}^2`$, where $`\rho _0`$ is the ambient gas density, and $`v_\mathrm{s}`$ is the shock velocity. Taking $`\rho _0`$ = 0.26 $`m_\mathrm{p}`$ g cm-3 at the positions of the supershells (KI97), we obtain $`v_\mathrm{s}`$ = 435 km s-1. This is slightly higher than the value of $``$ 300 km s-1 estimated for NGC 4258. The shock velocity puts a minimum timescale on the duration of the jet in this model, since there must be sufficient time for the jet to propagate out to the locations of the supershells at 12 kpc galactocentric radius. Thus, the minimum duration of the jets is 2.7 $`\times `$ 107 yr for the parameters in this illustration.

The mechanical power of the jet would then be $`L_\mathrm{j}`$ = $`1/2\rho _0\pi R^2v_\mathrm{s}^3`$ = 5 $`\times `$ 1041 erg s-1, where the effective jet radius, R has a typical value of 1 kpc (e.g. Cecil et al 1995). This is similar to the jet power estimated for the well known spiral galaxy M 51, which is at the low end of the scale of nuclear activity (cf. Ford et al 1985). In NGC 3556, the integrated mechanical luminosity of the jet over its minimum lifetime is then $``$ 4 $`\times `$ 1056 erg. This lower limit is adequate to account for the observed kinetic energy associated with the larger eastern supershell (Sect. 2.3.1). Here we have assumed that the efficiency for converting the input energy into the kinetic energy of the shell is of order 1% as usually taken for the multiple supernova model. For a higher efficiency, the required jet power can be lower.

The expansion of a typical bubble in the ISM should then proceed similarly to previously modelled scenarios, the difference being the source of input energy. Since the input energy is much higher than conventional supernovae, the bubbles are more likely to achieve blowout, providing a natural conduit through which cosmic rays (including those supplied by the jets) can escape into the halo. Thus, the presence of supershells is expected to correlate with the existence of a nonthermal radio halo.

The locations of the anchor points of the two supershells in NGC 3556 suggest that the postulated blasting of each supershell would have occurred at a galactocentric distance which is just past the peak of the HI rotation curve (see Fig. 6 of KI97). The radially outward expansion of each supershell would then expose it to the velocity shear in the medium, exerting a side-ways wind pressure counter to the galactic rotation. Due to this, the half of the supershell towards the direction of the galactic rotation would be compressed and hence brightened, whilst the other half would be dragged out due to the velocity shear in the external medium and, consequently, dimmed. Such a deformation of the two supershells from spherical symmetry, in course of their expansion, would give rise to a Z-symmetry, which is consistent with the observations (Sect. 2.3.1; KI97).

The cooling time for the X-ray emitting sheaths around the jets of NGC 4258, assuming unity filling factor, is only 4 $`\times `$ 106 yr (Cecil et al. 1995) and therefore, once the jets have turned off, such a signature of their existence would disappear rapidly. The lifetime of the synchrotron emitting particles in the jets is longer, typically from 106 to 107 yr. For NGC 4258, it has been estimated to be between 1 to 5 $`\times `$ 107 yr (Martin et al. 1989). However, shorter timescales are possible and will depend on local magnetic field strength and spectral index. For instance, the magnetic field in the north-east radio jet in NGC 1068, which extends to $``$ 450 pc from the nucleus, is $``$ 5 $`\times `$ 10-4 Gauss and the spectral index is -1 (Wilson & Ulvestad 1987). This gives a synchrotron lifetime of only 1.5 $`\times `$ 105 yr. The particle lifetime is expected to be shorter in shocks where the magnetic field is compressed. Thus, we expect the radio jet in NGC 3556 to fade over a timescale less than the ages of the supershells.

Similar arguments apply to the radio core which is expected to decay faster than the jets, once the nuclear activity has switched off, given the likelihood of a stronger magnetic field and flatter spectral index in the core. Even if some radio emission from a core were to persist after the activity ceases, it may be below the rms noise level of the available maps. Chary & Becklin (1997), for example, estimate the radio luminosity of the core of NGC 4258, which is known to have a supermassive central object of mass 3.6 $`\times `$ 107 M and a VLBI jet, to be Lrad = 100 L. The 2$`\sigma `$ noise level of the radio maps of Irwin et al. (2000) corresponds to a radio power of 6.8 $`\times `$ 1017 W Hz-1. Integrating over 1011 Hz, this yields a radio luminosity of 180 L. Thus, even if NGC 3556 harbours a radio core of the same luminosity as NGC 4258, it would not have been detected in the observations mentioned above.

## 4 Summary and Possible Implications

In this study we have sketched a scenario whereby radio jets ejected during an active nuclear phase in the life of a large spiral galaxy could inflate out of the disk exceptionally large shells, i.e. supershells of neutral hydrogen gas. When applied to the Scd galaxy, NGC 3556, this model can account for each of the several features enumerated in Sect. 2.3.1. In the context of the HI supershells seen in this galaxy, this scheme appears to fare distinctly better than either of the two standard models for supershells which invoke either a starburst induced superwind or an infall of external gas clouds on to the galaxy disk. Thus, the jet-ISM interaction scenario for disk galaxies could effectively supplement the other two mechanisms already proposed for the supershell phenomenon, with a greater relevance whenever the shells are extraordinarily large as well as energetic and occur at large galactocentric distances.

The energy requirement is no longer a major issue in the present model. Even a modest energy input from a 107 M compact object at the galactic nucleus, such as that observed in NGC 4258, is more than adequate to supply the input energy required for the supershells of NGC 3556. This model can therefore account for HI supershells which are difficult to form via conventional processes involving massive star formation, or those occurring in galaxies lacking potentially interacting companions. Larger supershells (or parts thereof, given the probable development of instabilities) may also be predicted, since known central masses and/or accretion rates may well surpass those considered here.

Furthermore, the present model can account for the symmetric HI features in galaxies. For instance, it is more likely that bubbles will be inflated by the propagating jets at galactocentric radii where the ISM density has declined to sufficiently low values for blow-out. At the same time, symmetry is not a generic outcome our model, since local shocks, dictated by the interaction of the jet with local density perturbations, can be important in determining where the bubbles will be inflated. For example, known bubbles in the lobes of radio galaxies (Sect. 2.2) are not always seen to be equidistant on either side of the nucleus.

In our sample illustration (Sect. 3.2), the jets must exist for long enough to reach large galactocentric radii (3 $`\times `$ 107 yr) but need to inject energy over only a fraction of the lifetime of the supershells (2 $`\times `$ 107 yr, or even less, if the efficiency exceeds 1%). The jets presumably switch off thereafter, with the synchrotron emitting particles decaying in 106 \- 107 yr. It is difficult to predict whether jets and HI supershells should co-exist in spirals. Nuclear activity in galaxies is commonly believed to be a transient phenomenon, lasting for $`10^710^8`$ years (e.g. Eilek 1996, Scheuer 1995; Soltan 1982; Haehnelt & Rees 1993; Ho et al 1997). If AGN activity occurs on timescales longer than a few $`\times `$ 107 yr, this would suggest that there might be disk galaxies in which both phenomena should be observed at the same time. However, such nuclear activity timescales are likely to be more representative of radio galaxies; AGN activity in spirals is probably shorter lived. An important next step is to consider what jet parameters are required to reach large, or at least kpc-scale radii, in traversing the typical ISM of a spiral galaxy.

If jets in spirals both recur and precess, it may be possible to inflate bubbles over a variety of galactocentric radii and azimuthal angles. However, it is unlikely that the “frothy” nature of the HI seen in galaxies such as Holmberg II or the holes in M 31 could be produced by such jets directly without seeing some lingering evidence for the jets as well. Also, the jet phenomenon cannot directly account for features which correlate with the star forming disk, such as some known radio and H$`\alpha `$ halos. If there is a connection, it is more likely to be indirect. For example, once the expansion of the postulated jet-blown supershells is halted due to the gravitational pull of the host galaxy, their segments (i.e. HI clumps) would eventually shower back on to the galactic disk, giving rise to secondary shells, bubbles and cavities in the HI component of the disk. The impact could also trigger sporadic bursts of star formation, especially in the outer disk, as seen, e.g., in gas-rich low-surface-brightness galaxies (cf. O’Neil et al. 1997) which are thought to be the present-day remnants of quasars (e.g. Schombert 1998). Thus, much after its cessation, the nuclear activity in gas-rich galaxies may continue to influence the evolution of their disks.

###### Acknowledgements.

It is a pleasure to thank Paul J. Wiita for valuable comments on the manuscript. JI gratefully acknowledges a grant from the Natural Sciences and Engineering Research Council of Canada.


text:
# Universal quantum computation with two-level trapped ions

## I Introduction

Ion trap quantum computation, first introduced by Cirac and Zoller , is a potentially powerful technique for the storage and manipulation of quantum information . In this scheme, information is stored in the spin states of an array of trapped ions and manipulated using laser pulses. Reasonably long coherence times can be achieved, compared to achievable switching rates , and individual qubits can be addressed through spatial separation of the ions. Experimental implementations of this scheme have succeeded in performing simple two-qubit logic gates  and preparing entangled states .

An ion trap quantum computer may be modeled as a collection of $`N`$ particles with spin $`\frac{1}{2}`$ in a one-dimensional harmonic potential. Laser pulses incident on the ions can be tuned to simultaneously cause internal spin transitions and vibrational (phonon) excitations, thus allowing local internal states to be mapped into shared phonon states. In this manner, quantum information can be communicated between any pair of ions and logic gates can be performed.

In this paper, we consider an interesting question which arises in this scenario: what is the simplest internal spin state structure required by each ion? In the original Cirac-Zoller formulation, ions with three levels are required. However, as the Cirac-Zoller method is generalized to other physical systems, such as neutral trapped atoms or quantum dots in an electromagnetic cavity, it has become highly desirable to determine whether just two internal levels are sufficient for performing universal quantum computation.

Previously, Monroe et al. have shown how a controlled-not can be performed between an ion and a phonon state using only two-level ions . Their method depends on fine-tuning of the Lamb-Dicke parameter, which relates the laser frequncy to the scale of the ions’ wave functions, to cancel unwanted side effects. Furthermore, they do not allow the Lamb-Dicke parameter to be arbitrarily small; the lowest value quoted in  is 0.316.

Here, we provide a general and accessible technique for performing universal logic between ions with only two internal levels. This scheme operates with any sufficiently small value of the Lamb-Dicke parameter, and also introduces new ways to utilize (or avoid) particular phonon states while performing quantum logic gates. We begin in Section II by presenting the allowed operations using the usual Jaynes-Cummings model for spin-boson interactions. We then describe a multiple-pulse construction for the controlled-not gate in Section III and give further constructions in Section IV before concluding with some possible extensions of the work.

## II Theoretical ion trap model

The energy level diagram of one ion, including the motional state, is shown in Fig. 1. For the sake of definiteness, we may think of the motional levels as corresponding to the center of mass degree of freedom. However, they could just as well correspond to another mode of oscillation, such as the “breathing” mode of a pair of ions. In practice, one might wish to choose a mode other than the center of mass to achieve reduced susceptibility to decoherence .

The Hamiltonian of this system is $`H_0=\mathrm{}\omega _0\frac{\sigma _z}{2}+\mathrm{}\omega _za^{}a`$, where $`\sigma _z`$ is a Pauli spin operator for the nuclear spin and $`a`$ annihilates a phonon. Throughout this paper, we work in the frame of this Hamiltonian. Turning on the electromagnetic field of a laser gives an interaction Hamiltonian

$$H_I=\stackrel{}{\mu }\stackrel{}{B},$$
(1)
where $`\stackrel{}{\mu }=\mu \stackrel{}{\sigma }/2`$ is the magnetic moment of the ion and $`\stackrel{}{B}=B\widehat{x}\mathrm{cos}(kz\omega t+\mathrm{\Phi })`$ is the magnetic field produced by the laser. Here $`z=z_0(a+a^{})`$, where $`z_0=\sqrt{\mathrm{}/2Nm\omega _z}`$ is a characteristic length scale for the motional wave functions and $`m`$ is the mass of an ion.

We consider the regime in which $`\eta kz_01`$. In this regime, we may determine the effect of a laser pulse at a specific frequency $`\omega `$ by expanding Eq. (1) in powers of $`\eta `$ and neglecting rapidly rotating terms. Pulsing on resonance ($`\omega =\omega _0`$) allows one to perform the transformation

$$R(\theta ,\varphi )=\mathrm{exp}\left[i\frac{\theta }{2}(e^{i\varphi }\sigma ^++e^{i\varphi }\sigma ^{})\right],$$
(2)
allowing one to do arbitrary single-qubit operations on an ion’s internal state. Pulsing at the $`n`$th blue sideband frequency ($`\omega =\omega _0+n\omega _z`$) gives

$$R_n^+(\theta ,\varphi )=\mathrm{exp}\left[i\frac{\theta }{2}(e^{i\varphi }\sigma ^+a^n+e^{i\varphi }\sigma ^{}a^n)\right]$$
(3)
and pulsing at the $`n`$th red sideband frequency ($`\omega =\omega _0n\omega _z`$) gives

$$R_n^{}(\theta ,\varphi )=\mathrm{exp}\left[i\frac{\theta }{2}(e^{i\varphi }\sigma ^+a^n+e^{i\varphi }\sigma ^{}a^n)\right].$$
(4)
Here, $`\sigma ^\pm =(\sigma _x\pm i\sigma _y)/2`$ act on the internal state of the ion. In each case, the parameter $`\theta `$ depends on the strength and duration of the pulse and $`\varphi `$ depends on its phase. For an $`n`$th order transition of duration $`t`$, $`\theta `$ is given by

$$\theta =\frac{\mu Bt\eta ^n}{2\mathrm{}n!}.$$
(5)
Also,

$$\varphi =\mathrm{\Phi }+(n\mathrm{mod}\mathrm{\hspace{0.17em}4})\frac{\pi }{2}.$$
(6)
## III CNOT using the first motional sideband

We suppose that the state of the ion trap quantum computer begins in the motional $`|0`$ state. We may excite higher-order motional levels using $`R_n^\pm `$ so long as we return to the motional ground state at the end of the gate. It is simplest to consider cases where only a few higher-order levels are relevant. For a particular ion, we refer to the subspace spanned by $`\{|00,|01,|10,|11\}`$ (using the notation of Figure 1) as the computational subspace (CS).

First, we note that for the special values $`\theta =j\pi \sqrt{2}`$, where $`j`$ is an integer, $`R_1^\pm `$ preserve the CS — that is, they map states within the CS to other states within the CS. For example,

$$R_1^+(\pi \sqrt{2},\pi /2)=\left(\begin{array}{cccc}\mathrm{cos}\frac{\pi }{\sqrt{2}}& 0& 0& \mathrm{sin}\frac{\pi }{\sqrt{2}}\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ \mathrm{sin}\frac{\pi }{\sqrt{2}}& 0& 0& \mathrm{cos}\frac{\pi }{\sqrt{2}}\end{array}\right).$$
(7)
Next, note that it is easy to produce another gate which preserves the CS by conjugation. For example, consider the gate $`G=G_{\mathrm{in}}G_{\mathrm{out}}`$, where $`G_{\mathrm{in}}`$ has support only on $`\{|00,|10,|11\}`$ and $`G_{\mathrm{out}}`$ acts on the rest of the space. Since $`R_1^+`$ does not connect these subspaces, conjugating this gate by $`R_1^+`$ produces another gate which preserves the CS. In other words,

$$R_1^+(\theta ,\varphi )GR_1^+(\theta ,\varphi )$$
(8)
preserves the CS. Using this observation, we can easily diagonalize Eq. (7), giving the gate

$`P`$ $`=`$ $`R_1^+(\pi /2,0)R_1^+(\pi \sqrt{2},\pi /2)R_1^+(\pi /2,0)`$ (9)
$`=`$ $`\mathrm{diag}(e^{i\pi /\sqrt{2}},1,1,e^{i\pi /\sqrt{2}}).`$ (10)
This gate uses first-order sideband pulses for a total duration of $`(1+\sqrt{2})\pi `$, only slightly longer than the $`2\pi `$ duration required for the equivalent step in the Cirac-Zoller scheme.

With this diagonal gate and $`R_1^{}`$, it is straightforward to construct a controlled-not gate between two ions using a sequence similar to the original Cirac-Zoller construction. Note that $`R_1^{}(\pi ,\pi /2)`$ allows us to interchange $`|01`$ and $`|10`$ (up to a phase), performing

$$\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right)$$
(11)
on the subspace $`\{|00,|01,|10\}`$. Assuming the motional state is initially $`|0`$, this corresponds to a swap between the state of the ion and the motional state. In fact, we can use a version of this gate with arbitrary phase to do logic between ions. Thus we can perform a controlled-not from ion $`j`$ (the control) to ion $`k`$ (the target) using

$`CNOT_{jk}=`$ $`Z_j(\pi /(2\sqrt{2}))R_{1j}^{}(\pi ,\varphi )`$ (13)
$`H_kP_kZ_k(\pi /(2\sqrt{2}))H_kR_{1j}^{}(\pi ,\varphi ),`$
for any value of $`\varphi `$, where a subscript denotes which ion is acted on and we have introduced the single-qubit gates

$$Z(\varphi )=\left(\begin{array}{cc}e^{i\varphi }& 0\\ 0& e^{i\varphi }\end{array}\right),H=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).$$
(14)
## IV Other gate constructions

Although the construction by which we arrived at Eq. (10) is relatively straightforward, there are other ways to construct gates which preserve the CS. In fact, it is possible to perform a generalization of Eq. (10) by considering the form

$`U(\varphi )=`$ $`R_1^+(\alpha ,\pi /2)R_1^+(\beta ,\gamma )R_1^+(2\varphi ,\delta )`$ (16)
$`R_1^+(\beta ,\gamma )R_1^+(\alpha ,\pi /2).`$
Making the choices

$`\mathrm{cos}\alpha `$ $`=`$ $`\mathrm{cos}\sqrt{2}\alpha `$ (17)
$`\mathrm{cos}\beta `$ $`=`$ $`\mathrm{cos}\sqrt{2}\beta `$ (18)
$`\mathrm{sign}(\mathrm{sin}\alpha )\mathrm{sign}(\mathrm{sin}\beta )`$ $`=`$ $`\mathrm{sign}(\mathrm{sin}\sqrt{2}\alpha )\mathrm{sign}(\mathrm{sin}\sqrt{2}\beta )`$ (19)
$`\delta \gamma `$ $`=`$ $`\mathrm{sin}^1\left({\displaystyle \frac{\mathrm{cos}\alpha }{\mathrm{sin}\beta }}\right)`$ (20)
$`\gamma `$ $`=`$ $`\mathrm{tan}^1[\mathrm{cos}\beta \mathrm{tan}(\delta \gamma )]`$ (21)
results in

$$U(\varphi )=\mathrm{diag}(e^{i\varphi },e^{\pm i\sqrt{2}\varphi },1,e^{i\varphi }).$$
(22)
For example, there is such a solution with $`\alpha 298.2^{}`$, $`\beta 149.1^{}`$, $`\gamma 63.87^{}`$, and $`\delta 131.0^{}`$.

The source of the factors of $`\sqrt{2}`$ in gates derived from the first-order sideband is the matrix element $`1|a|2=\sqrt{2}`$ for harmonic oscillator states. However, note that

$$1|a^3|4=20|a^3|3=2\sqrt{6}.$$
(23)
The integer ratio of these matrix elements suggests that third-order transitions may be used to create simpler diagonal gates. Indeed, we find

$$R_3^{}(2\pi /\sqrt{6},\varphi ^{})=\mathrm{diag}(1,1,1,1)$$
(24)
for any value of $`\varphi ^{}`$. This differs from a controlled-$`Z`$ gate by only single-qubit operations, so we can simply use the Cirac-Zoller construction to produce a CNOT between ions:

$`CNOT_{jk}=`$ $`H_kZ_j(\pi /2)R_{1j}^{}(\pi ,\varphi )R_{3k}^{}(2\pi /\sqrt{6},\varphi ^{})`$ (26)
$`Z_k(\pi /2)R_{1j}^{}(\pi ,\varphi )H_k`$
Here, we require fewer pulses than in Eq. (13). However, note that the third-order sideband pulse must be longer by a factor of order $`\eta ^2`$ than the first-order sideband pulses for the same laser intensity.

Most of these results generalize to other choices of the computational subspace. For example, consider using $`\{|00,|02,|10,|12\}`$. Analogous to Eq. (13), we find

$`CNOT_{jk}=`$ $`Z_j(\pi /(2\sqrt{6}))R_{2j}^{}(\pi /\sqrt{2},\varphi )H_k`$ (30)
$`R_{2k}^+(\pi /(2\sqrt{2}),0)R_{2k}^+(\pi /\sqrt{3},\pi /2)`$
$`R_{2k}^+(\pi /(2\sqrt{2}),0)Z_k(\pi /(2\sqrt{6}))`$
$`H_kR_{2j}^{}(\pi /\sqrt{2},\varphi ).`$
Similarly, analogous to Eq. (26) (exploiting the integral relationship $`2|a^7|9=60|a^7|7`$),

$`CNOT_{jk}=`$ $`H_kZ_j(\pi /2)R_{2j}^{}(\pi /\sqrt{2},\varphi )`$ (33)
$`R_{7k}^{}(\pi /(6\sqrt{35}),\varphi ^{})Z_k(\pi /2)`$
$`R_{2j}^{}(\pi /\sqrt{2},\varphi )H_k.`$
Finally, we wish to point out that it is conceivable to treat the lowest two motional levels as an additional qubit, rather than simply an intermediary for logic, if one can perform a true swap operation over the entire computational subspace. Note that

$`R_1^{}(2l\pi \sqrt{2},\pi /2)`$ (34)
$`=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}(l\pi \sqrt{2})& \mathrm{sin}(l\pi \sqrt{2})& 0\\ 0& \mathrm{sin}(l\pi \sqrt{2})& \mathrm{cos}(l\pi \sqrt{2})& 0\\ 0& 0& 0& 1\end{array}\right).`$ (35)
We may get arbitrarily close to a swap operation by choosing large enough $`l`$ with $`\mathrm{cos}(l\pi \sqrt{2})0`$. For example, choosing $`l=2378`$ gives a swap operation to a precision of about $`10^3`$. Although this observation may not be useful in practice, it shows that universal quantum logic including a motional qubit is not forbidden in principle.

## V Conclusions

We have demonstrated the possibility of doing ion trap quantum computation with two-level ions, allowing any sufficiently small value of the Lamb-Dicke parameter. While we believe that a construction similar to Eq. (13) will be most useful in practice, we have shown that there are many other ways to realize a controlled-not gate between ions without using extra levels, some of which might be useful given appropriate experimental conditions.

Our initial motivation for gates such as Eqs. (10) and (22) came from the theory of composite pulses used in the art of NMR, in which errors in pulse length or frequency are canceled at low order using a sequence of pulses. Thus one may accomplish with a composite pulse, in the presence of errors, what a single pulse would have performed in the absence of errors . In the present case, the goal is somewhat different: we construct a sequence of gates which take the system into higher motional levels and back again to perform some logical operation. Nevertheless, the idea of stringing together several pulses to perform a simple gate has proved fruitful. It is our hope that further progress can be made in the experimental implementation of quantum information processing devices by considering how existing tools can be applied to different physical implementations.