text:
# Soft Physics, Centrality and Multiplicity at RHIC

## I Introduction

A possible conclusion to be drawn from initial observations of inclusive charged particle spectra in Au + Au collisions at RHIC, by the four detectors PHOBOS , PHENIX , BRAHMS  and STAR , which observations, of course, are only of the final state of the collision, is that soft, strongly interacting QCD plays a dominant role in the dynamics. To be sure, there do exist potentially interesting anomalies in high $`p_{}`$ charged  and neutral  pion spectra, among others. These occur for $`p_{}`$ in the range $`2`$$`4`$ GeV/c, involve only a small fraction of the integrated cross-section, and may in fact be describable by an energy loss mechanism for $`\pi `$ mesons acting in the later stages of a collision.

The hadronic simulation LUCIFER  has already given a good description  of the earliest PHOBOS  and STAR  measurements of averaged mid-rapidity cross-sections for charged particles, $`(dN/d\eta )^{ch}`$ for $`|\eta |1`$, including also a successful prediction  of the surprisingly slight ($`10`$$`14`$%) rise in this observable as the collision energy is increased from $`\sqrt{s}=130`$ GeV to 200 GeV. Also correctly anticipated was the shape of the pseudo-rapidity spectrum for charged particles at higher values of $`\eta `$ .

Here, we consider simulation results for the centrality dependence which is found by all of the detectors . These collaborations emphasize the variation of $`(dN/d\eta )^{ch}|_{\eta =0}`$ with the number of participants. This work concerns itself with comparisons of the spectral magnitudes and shapes as a function of centrality. The variation of experimentally defined centrality is directly related theoretically to changes in the impact parameter $`b`$. The energy and geometrical dependence extracted from LUCIFER speak directly to assumptions underlying the cascade simulation, and to the nature and the number density of the resonance rescatterers employed in the model.

In earlier work  we stressed that the role of our simulation was only to define the ‘ordinary,’ giving experimentalists a background calculation against which to judge potentially unusual data. There do, of course, exist a variety of other simulations , having similar goals but differing approaches: some are purely partonic models and some are hybrid partonic-hadronic.

It seems again appropriate to compare the, by now, extensive observations at RHIC with simulations which do not assume that plasma is necessarily present. The purest such comparison would employ a model involving only colour singlet, hadronic degrees of freedom. The physical picture is one in which ion-ion collisions are describable in the main by multiple interactions between excited hadrons only. In such a picture, constituent quarks, present within original and produced baryons, are perhaps raised to states differing from those present in the lowest mass baryons, but glue holding the valence quarks in place remains ‘sticky.’ So quarks and gluons continue to act as if still confined within hadrons. In addition to excited baryons, color singlet bosonic degrees of freedom are produced and play a significant role. Such a description has proven feasible in the Pb + Pb collisions examined in NA49 . In fact, excited bosons only enter at a later stage of the simulation, often after some initial baryon-baryon interactions are completed. It remains to be seen whether at the higher energies of RHIC a large number of partons, carrying a large fraction of the energy density, perhaps ‘minijets’ , are in fact free and able to roam over large spatial distances, and more importantly whether sufficient ‘free’ gluons are present to create the thermodynamic basis for material which could be described as quark-gluon plasma.

The simulation we employ has two successive stages. In phase I, no energy is removed for soft processes. Therefore, initial baryon-baryon interactions take place very nearly on the light cone. The construction of a collision history is guided by inputs from NN data, the key ones being: (1) absolute total cross-sections as a function of energy, (2) branching ratios for elastic and various inelastic processes and (3) multiplicity distributions for meson production.

Our two body model has been discussed extensively in earlier work . It attempts to emulate multi-peripheral models of nucleon-nucleon scattering which divide the total cross-section into elastic, single diffractive and non-single diffractive  processes. Our assumption is that generic $`\rho `$-like and $`K^{}`$-like resonances are the major rescatterers, only entering into phase II, of the simulation, which begins after a formation time $`\tau _f`$, a central model parameter, has elapsed for particles produced in phase I, in their various rest frames. Since these resonances themselves decay into multiple particles stable under the strong interaction on a time scale $`\tau _d`$, which is a second important model parameter, their multiplicities must be adjusted appropriately so as to reproduce free space NN data, and they are in fact reduced in number by some two to three times relative to the expected multiplicity of stable mesons. Importantly however, it then follows that any model reproducing the known nucleon-nucleon data, such as explicit string models , could have been employed up to the commencement of phase II of the simulation.

We use a resonance model of nucleon-nucleon scattering for reasons which will become clear as we proceed and because interactions which are particle-like in nature are easier to deal with than, say, the interactions of strings. We will discuss the spatial and energy densities of the ‘rescatterers’ at some length in what follows and thus make evident a simplicity in both the energy and centrality dependence of the observed pseudo-rapidity densities which arises from such considerations. It appears that over an appreciable range of both energy and impact parameter, the density of generic hadronic resonances, built up during phase I, must remain nearly constant in ion-ion collisions, if the observed centrality and energy dependence of final charged particle spectra are to be correctly described. The extracted resonance size, assuming close packing or saturation, is eminently reasonable.

We again adopt the position that the analysis for Au + Au at RHIC simply presents an extrapolation from earlier NA49 inclusive measurements  to the considerably higher energy RHIC collisions. We could equally well have normalized the absolute level of meson production to the lowest RHIC energy ($`\sqrt{s}=56`$ GeV). This normalization directly implies a determination of the average density of generic resonances, or rescatterers, in phase II of the LUCIFER simulation for Au + Au, a theoretical observable which is then a single number characterising presently measured inclusive RHIC physics.

Uniformity in resonance production is achieved by imposing an upper limit on total multiplicity in contiguous groups of interacting nucleons within an A + A event  or, equivalently, by imposing a limit directly on the two body interactions themselves. This is accomplished by introducing this limit at a single energy. The density of the saturated resonance fluid is fixed by normalising the simulated multiplicity to the experimental multiplicity at a fixed total energy in a massive (large A) ion-ion collision. No such restriction is imposed on pure NN interactions, or for example p + D collisions. Previously,  we based our approach on ideas put forward by Gottfried  in describing p + A collisions. Gottfried delayed interactions of particles produced in the successive NN interactions in a p + A collision, until the produced particles were physically sufficiently separated from each other and their progenitors to act as separate hadrons. The extra constraint which we impose, saturation of the rescatterer density in ion-ion collisions, may be occasioned when multiple resonances produced in neighbouring NN collisions interfere with one-another. It is almost self evident that individual resonances condensing from such two body collisions cannot overlap.

After a short digression on the details of LUCIFER, we develop the above ideas in two ensuing sections. First, comparisons are presented between the model and the experimental results at RHIC from the various detectors, large and small. These include charged particle pseudo-rapidity measurements at $`\eta =0`$ for varying centrality, at the first RHIC energy $`\sqrt{s}=130`$ GeV analysed in this fashion. Also treated is the energy dependence for $`(dN/d\eta )^{ch}|_{\eta =0}`$ at $`\sqrt{s}=56`$, $`130`$, and $`200`$ GeV, using for the lowest energy, data which are available only from the PHOBOS  collaboration. Also considered at $`\sqrt{s}=130`$ GeV, are the full pseudo-rapidity distributions for differing degrees of centrality.

Secondly, the distributions of underlying rescatterers, i. e.  the dependence of the number of boson-like resonances, $`N_{RES}(b,s)`$, on impact parameter and energy, are extracted from the simulations. The simple relations for the spatial and energy distributions thus obtained illuminate the physical basis for the overall simulation. In particular, the number density of meson resonances reached at the formation time $`\tau _f`$, although not large, is sufficient to explain the observed stable meson production. The baryon density from the onset of an AA event up to the earliest formation time t=$`\tau _f`$, involves far fewer particles, but could still in principle provide problems. The latter will be discussed at an appropriate point.

## II LUCIFER: Dynamics and Inputs

For completeness we present a brief outline of the dynamics in LUCIFER . Phase I has been adequately described above and sets the initial conditions for an ordinary hadronic cascade (with energy loss) in phase II . The numbers and the properties (flavour, mass, etc. ) of hadrons entering into phase II are determined in phase I by recording all baryon-baryon collisions. Providing complete initial conditions for phase II requires us to specify four momenta and positions for the produced resonances. The positions are simply randomly distributed along the paths of nucleons colliding in phase I, while the momenta must conserve overall energy and momentum and be consistent with production and phase space in NN collisions. Existing baryons acquire transverse momentum by making a random walk in the number of collisions, and softly produced meson resonances are given transverse momenta commensurate with soft production as observed in NN collisions. The total multiplicity of the resonances is restricted as described above. We reemphasise that the generic resonances first appear only in phase II. The latter is a normal cascade in which energy is lost at each binary collision and during which the resonances both interact with each other and decay into pions and other stable particles. The restricted range of masses for the resonances, 300-1200 MeV, or somewhat higher if strangeness is involved, guarantees we have only soft processes, i. e.  that only relatively low $`p_{}`$ is generated in the decays.

We have listed the important two-body inputs in phase I, but we should make more explicit a very important feature in the elementary model, viz the energy dependence of the multiplicity distributions, which obey KNO-scaling . Our description of the scaling of the multiplicity distributions comes from fitting multi-prong data mostly drawn from the UA(5) experiment . Central to this parameterization is the energy dependent average number of prongs $`\overline{n}(s)`$ entering through the definition of KNO scaling

$$n=\overline{n}(s)f(\frac{n}{\overline{n}(s)}).$$
(1)
Although we artificially limit the number of mesons produced in ion-ion collisions by our no-overlap constraint on the generic resonances, the energy dependence for increasing production of these particles still drives the final multiplicities as a function of $`s`$. In Figure 1 we display a De facto prediction of the two body model, i. e. which arises via the energy dependence of $`\overline{n}(s)`$. This is the comparison of calculated energy variation of central charged multiplicity, $`(dN/d\eta )^{ch}(s)|_{\eta =0}`$ with data . Since one is here dealing with a collision involving two incoming nucleons (participants) this quantity is, in effect, already divided by 1/2 times the number of participants. The slight increase predicted for the rise in $`(dN/d\eta )^{ch}|_{\eta =0}`$ for Au + Au, going from $`\sqrt{s}=130`$ GeV to $`200`$ GeV, arises essentially completely from the energy dependence of $`\overline{n}`$: this dependence is weak in the two body data and equally weak for ion-ion collisions. This is a dynamic, and not merely kinematic, outcome of the simulation, but it is one which follows from the NN character of phase I.

## III Centrality Dependence of Rapidity Distributions: Lucifer vs Data.

Figure 2 reviews earlier $`(dN/d\eta )^{ch}|_{\eta =0}`$ results for the most central collisions, indicating the predicted $`12`$-$`14`$% rise between $`\sqrt{s}=130`$ and $`200`$ GeV for Au + Au at RHIC . A more complete picture of the variation in this pseudo-rapidity density $`(dN/d\eta )^{ch}|_{\eta =0}`$ for varying centrality cuts is shown in Figure 3 with LUCIFER compared to the existing detector data  for the $`130`$ GeV runs. The centrality cuts on the simulation are defined in terms of the calculated $`dN/db`$, which is integrated to give the appropriate degree of centrality between limits $`b_1`$ and $`b_2`$, $`\mathrm{\Delta }b=b_2b_1`$. Theoretical curves describing this $`b`$ dependence are shown in Figure 4. A more expeditious definition of the percentage centrality was employed earlier. This convenient geometrical definition is not necessarily in accord with actual experimental usage. The agreement between simulation and measurement persists out to rather high impact parameter, possibly deteriorating somewhat for $`b8`$ fm. For peripheral collisions the ‘constancy’ of resonance density, discussed in detail below, is not maintained in our saturated simulations. Indeed for the most peripheral collisions, with say only a few nucleons colliding, there should be no multiplicity suppression.

One should also confront the full pseudo-rapidity distributions insofar as they have been measured; this we do in Figure 5 using PHOBOS data at $`\sqrt{s}=130`$ GeV . The higher energy $`\sqrt{s}=200`$ yields similar comparisons. Again the agreement in shape is seen to be adequate. So the most inclusive meson data, taken at $`\sqrt{s}=130`$ GeV, is well described by the present two phase cascade. Finally, to demonstrate the similarity with other presentations of the data we include a plot, Figure 6, of participant nucleon number, $`N^{part}(inelastic)`$, in LUCIFER versus the ‘number of participants’ which was abstracted from two of the experiments. In the simulation we can easily estimate the total $`N^{part}(inelastic)`$ by counting the number of spectator nucleons (they are just those initial nucleons which had no inelastic collisions in phase I or II of the cascade) in a sample of events having some fixed $`b`$ and subtracting that number from the total number of initial nucleons. Of course only inelastic participants, which contribute to the multiplicity, must be selected for comparison with experiment where the distinction between ‘inelastic’ participants and ‘elastic’ participants becomes moot, so that another subtraction of purely elastic participants must be made theoretically. For the degree of centrality reached in the experiments, $`50\%`$, this additional adjustment is small, but it would become an appreciable correction, tending to $`20\%`$, in more peripheral events.

## IV Resonance Multiplicity Dependences on $`b`$ and $`\sqrt{s}`$

### A Fire Cylinder and Transverse Areal Density

Most of the simulation results can be understood in terms of simply constructed model quantities. At the commencement of the ‘soft’ cascade, phase II, we constrained the system to avoid spatial overlap of the rescatterers. The extent to which this succeeds can be quantified as a function of impact parameter. Plotted in Figure 7 is the distribution of the $`\rho `$-like resonances as a function of $`b`$, at the start of phase II. There is of course a dispersion in the numbers $`N_{RES}(b,s)`$, reflecting fluctuations from event to event, but a good estimate of the areal density of resonances can still be made. We use the term ‘areal’ density in referring to the number of resonances in the area of overlap of the two nuclei, obtained when the overlap volume from phase I is projected into a plane transverse to the beam direction, for fixed impact parameter.

The overlap volume or interaction region has, in our simulation of a gold-gold collision, for early times $`t>\tau _f`$ very nearly the form of a generalised cylinder. The directrix of this cylinder is a convex curve formed by the two circular arcs bounding the projection of the overlap volume of the spherical nuclei into the transverse plane. The generator of the cylinder is parallel to the beam direction and along the direction of longitudinal expansion of the interaction region. Perhaps this volume could be appropriately be termed a ‘fire cylinder.’ The area in question then, is the transverse cross-section of the fire cylinder. Taking the geometry of the two colliding nuclei to be sharply spherical, with radius R, the area of overlap $`A(b)`$ can be calculated:

$$A(b)=2R^2\left[cos^1\left(\frac{b}{2R}\right)\left(\frac{b}{2R}\right)\sqrt{1\left(\frac{b}{2R}\right)^2}\right].$$
(2)
For $`bR`$, the areal density of resonances is then found to be:

$$\sigma _a(b,s)=N_{RES}(b,s)/A(b).$$
(3)
Sample curves obtained for this density are exhibited in Figure 8 for varying choices of nuclear radius $`R=6.8,7.0,7.2`$ fm, the latter value being closer to the actual situation for gold nuclei, especially when the addition of a surface diffusivity converts the nuclear geometries into the more realistic forms we employ in the Monte-Carlo. This figure demonstrates, however, the near constancy of $`\sigma _a(b)`$ out to rather large separations of the colliding nuclei, $`b8`$$`10`$ fm. This perhaps surprisingly accurate, but intuitively satisfying, uniformity of the resonance population of the fire cylinder is a representation of the constraints placed on production. The result of imposing constancy strictly out to $`b=10`$ fm is included in the above Figure 3. Measurements seem to be telling us that indeed the transverse density $`\sigma _a(b,s)`$ varies very little.

A reasonable approach for still more peripheral events would be to smoothly cut off the degree of suppression of the resonance multiplicity, until for the most glancing collisions there was no multiplicity suppression at all, at which point one expects after all that only unadulterated two body reactions occur. The results for the latter eventuality are shown in Figure 1, as explained previously a prediction of the elementary two-body model used here and extracted, in fact, from p+D reactions.

### B Fire Cylinder Volume and Energy Dependency

It is more difficult to precisely define a volume for the interaction region or fire cylinder at the onset of phase II of the simulation, the soft cascade. From figure 9 it is clear that after time $`t=\tau _f`$ in the center of momentum (lab) frame, up to the time of termination of phase II, new resonances are continually appearing, at times related to their velocities $`\beta _i`$ via $`t_i=\gamma _i\tau _f`$, $`\gamma _i`$ being the Lorentz factor peculiar to the $`i`$’th resonance. In the final analysis, a choice $`\tau _f=1.21.6`$ fm/c for the formation time was made, with globally similar results. In general, it is clear there is some interplay between $`\tau _f`$ and the overall normalisation of resonance number. A more extended delay time implies a more expanded system and less cascading in phase II. Too short a formation time could in fact lead to an inability for the simulation to keep within the experimental limits for total meson numbers. No such problem arises for the range of $`\tau _f`$ selected above, and these total numbers as a function of centality are well reproduced.

For these last issues of spatial densities we choose the most central collisions with $`b4`$, although little would be changed by using a fixed impact parameter. Figure 10 demonstrates the strong interdependence between the number of resonances $`N_{RES}(s)`$ and the average longitudinal Lorentz factor $`\overline{\gamma _z}(s)`$ (the averaging here is over all resonances). $`\overline{\gamma _z}`$ is directly related to the average length of the fire cylinder in the longitudinal direction:

$$L_z(s)=2\tau _f\overline{\gamma _z}.$$
(4)
The energy dependence of $`N_{RES}(s)`$ is then very simply and naturally described in Figure 10, and is another dynamical statement of the model. The cylinder length grows linearly with $`\overline{\gamma _z}`$ for increasing $`s`$, once again a useful and easily understood, result.

It is then an easy matter to compute the volume of the fire cylinder

$$V(b)=A(b)L_z(s),$$
(5)
and thus to obtain the three dimensional number density:

$$\rho (b,s)=\sigma _a(s)/L_z(s).$$
(6)
For this selection of very central collisions $`\rho (s)`$ is found to be $`0.701.10`$ particles per fm<sup>3</sup>, depending on the specific choices for nuclear radius and $`\tau _f`$ in Figure 8. This leads to an average resonance size of $`0.70`$ fm. Of course, in the simulation we treat the generic scatterers as points and their extracted size is simply a statement of the volume each resonance occupies, more or less at the time it is created. It is, however, significant that such a size is eminently reasonable for a hadronic resonance, and that the associated number density is consequently relatively low. These averages obtain of course after the formation time $`\tau _f`$ in the rest frame of produced mesons, but considerably later in say the c. of m. Nevertheless, to use an earlier, ‘average’ time would also involve less produced particles, and the fire cylinder densities we introduce would be little influenced.

## V Transverse Number and Energy Densities

One can also estimate the transverse energy densities, $`\rho _ϵ(t)`$, associated with the above number density $`\rho (b,s)`$. Considering, for now, only the dominant bosonic contribution, one arrives at a total $`\sqrt{(1)^2+(1)^2}`$ GeV i. e.  $`1.4`$ GeV/ fm<sup>3</sup>. This since the average bosonic mass is somewhat over $`1`$ GeV and its average $`p_{}`$ in phase I is somewhat under $`1`$ GeV/c. To this must be added the energy density from baryons present at the initiation of phase II, an increase within the fire cylinder imagined above of some $`20`$%; thus overall $`\rho _ϵ(t)1.7`$ GeV/$`fm^3`$. The density $`\rho _ϵ(t)`$ is an average over many events constructed at the start of phase II, and is well below the average Bjorken value quoted by the RHIC experimental collaborations. Almost independent of collider energy, a major proportion, some $`40\%`$ of the final mesons, and hence also a large fraction of the transverse energy, is produced via the collisions taking place in phase II  of the simulation. The variation in this fraction is only a few percent from $`56`$ to $`200`$ GeV/A. Combining the contributions from I and II would yield $`2.8`$ GeV $`fm^3`$, still short of the Bjorken estimate. Part of this discrepancy in total transverse energy densities is explained by our use of a longer formation time. However, since the final phase of the cascade occurs when the system has surely expanded considerably spatially, the extra transverse energy generated in II contributes little to the early density, often used as a trip wire in theoretical models for the initiation of ‘plasma’ creation.

The question of detailed double differential cross-sections, $`d^2N/d\eta p_{}dp_{}`$, especially those at the highest transverse momenta measured, will be examined thoroughly in ensuing work. For now we present in Figure 11 a preliminary comparison between STAR data  and LUCIFER, pertaining to total negative hadron ($`h^{}`$) production below $`p_{}=4GeV/c`$. The higher ranges of this comparison, say $`p_{}2`$ GeV, are possibly beyond the limits of our current modeling which largely omitted hard processes in the underlying two body scattering. Large $`p_{}`$ events certainly require more study, but the present results are nevertheless interesting. The highest $`p_{}`$ particles are produced in the simulation by multiple scattering and are thus subject to normal energy losses expected in such scattering. Despite this reservation, the rapid drop in $`d^2N/d\eta p_{}dp_{}`$ assures one that the softer part of the spectrum dominates total transverse energy production. Thus the total transverse energy $`E_t`$ generated in the cascade is close to that seen experimentally.

Of course we cannot really say what transpires in the earliest stages of the ion collision. We have made certain assumptions about hadron tenacity and saturation and have then obtained at the close of the second stage of the cascade conditions very much like those observed for the global, inclusive, production processes. One problem, mentioned in the Introduction, concerns the earliest phase of the collision, perhaps even before the production of bosonic resonances, and involves the purely baryonic densities. If one assumes each nucleus is never less than $`1`$ fm in transverse extent, then for $`b=4`$ fm collisions the baryon density does not rise much above one per fm<sup>3</sup>.

Further, the successful use of generic resonances to describe the cascading might simply reflect a correct fixing of average cross-sections and cross-section ‘flow’ during an ion-ion event.

## VI Comments and Conclusions.

The two phases of the LUCIFER simulation have been merged into a single, idealized, dynamical model. Such a model can be only a part of the complete story, of course; partonic degrees of freedom do exist and for sufficiently high particle transverse momentum presumably will dominate the inclusive spectra. Nevertheless, in our approach centrality and energy dependence comparisons between data and simulation are understood as arising from a common physical origin: the uniform and unvarying density of resonances in the fire cylinder. The constancy of areal and volume densities of bosonic resonances and their close to linear dependence in energy with $`\overline{\gamma _z}`$ are notable results of the simulation.

The dynamics of the collision is split into two phases, sequential in time: one could have created the initial conditions for phase II in some alternate fashion, for example by attributing the initial meson distributions to a string model or to de-confinement and hadronisation of partons . Our sole and central point here is, then, that signals other than those so far observed should show up to distinguish the present cascade from such alternatives, especially from nearly complete parton de-confinement. The string approach may in fact lead to the same results in simulating a massive ion-ion collision. We reiterate that the two-body modeling and cascading in phase I leaves strong footprints in the final results, especially in the surprisingly gentle increase in predicted, and measured, $`(dN/d\eta )^{ch}|_{\eta =0}`$ seen from $`\sqrt{s}=130`$ to $`200`$ GeV.

No one can doubt that there is a saturation effect observed in meson production at RHIC. To produce this effect, we reduce cascading in two fashions: (1) through the very presence of the boson resonances which eventually decay, over both phases I and II, into many more $`\pi `$ or $`K`$ mesons ($`4`$$`5`$), and (2) through the near constant resonance density in the fire cylinder achieved in central collisions. This latter can be imposed in principle via, and is equivalent to, placing a no spatial overlap constraint on the produced resonances prior to the soft cascade. The data, Figure 3, suggest a more strict uniformity than produced in the cascade.

In any case, the relatively small observed total multiplicities, which present a problem for all models, do not strongly support a large increase in the number of degrees of freedom in the early phase of Au + Au collisions at RHIC.

One should add that no attempt was made here to divide the output of the cascade into separate dependences on numbers of participants or on the numbers of binary collisions. The cascade does this dynamically and seamlessly. One might note, however, that the participant number directly determines the overall incoming energy put into an ion collision in phase I and thus also fixes the number of resonances at reinitialisation and commencement of phase II. The division into two phases does, however, distinguish the present approach from some alternatives; principally in reducing the generation of transverse energy at very early times.

We will return in later work to more specific questions, e. g. of flavour generation and high transverse momentum spectra, but we mention here that close to the correct ratios of strange meson production at mid rapidity, are obtained in the calculations. These contribute some $`10`$% of the particle number there.

This work addresses only gross, inclusive features of existing RHIC data. Accordingly, the results only suggest that non-standard medium effects must be looked for in more exclusive, rare events, whose analysis we all await. Places to look might include still higher transverse momenta spectra, direct $`\gamma `$-rays, $`\mu ^+\mu ^{}`$ and $`e^+e^{}`$ pairs, or possibly just a sampling of events containing unusually high multiplicity fluctuations.

## VII Acknowledgements

This manuscript has been authored under the US DOE grant NO. DE-AC02-98CH10886. One of the authors (SHK) is also grateful to the Alexander von Humboldt Foundation, Bonn, Germany and the Max-Planck Institue for Nuclear Physics, Heidelberg for continued support and hospitality. Useful discussion with the BRAHMS, PHENIX, PHOBOS and STAR collaborations are gratefully acknowledged.


format: text

text:
# Solar-like oscillations in 𝛿 Scuti stars

## 1 Introduction

The $`\delta `$ Scuti stars are in general main sequence stars with masses between 1.5 M and 2.5 M. They are located inside the classical instability strip (IS hereafter) where the $`\kappa `$-mechanism drives low-order radial and nonradial modes of low degree to measurable amplitudes (opacity-driven unstable modes). Only a small number of opacity-driven modes are observed in $`\delta `$ Scuti stars (for a review see e.g. Gautschy & Saio, 1996), but their amplitudes, which are limited by nonlinear processes, are much larger than stochastically driven intrinsically stable solar-like p modes.

For main-sequence stars with surface convection zones, located outside the IS, model computations suggest all modes to be intrinsically stable but excited stochastically by turbulent convection; for models located near the red edge of the IS the predicted velocity amplitudes become as large as 15 times the solar value (Houdek et al., 1999). Moreover, these computations suggest that models located inside the IS can pulsate simultaneously with modes excited both by the $`\kappa `$-mechanism and by the turbulent velocity field.

Provided that many modes can be detected, high-frequency p modes are more easily identified than low frequency p modes. Hence there are advantages of observing simultaneously both types of modes in the same star. As a first step, high-frequency p modes can help to determine the fundamental stellar parameters (e.g., luminosity, effective temperature) more accurately, whereas low-frequency modes, which are strongly sensitive to the properties of the deep layers of the star, can then be used as a diagnostic for the inner properties of the model. Such developments are outside the scope of the present paper and we only outline briefly the underlying idea.

The nearly regular frequency spacing of solar-like modes of high order (i.e., the large frequency separation) depends predominantly on the structure of the surface layers and consequently provides further constrains on the equilibrium models. Their degree $`l`$ and azimuthal order $`m`$ can be identified with the help of the classical echelle diagram method; this method was successfully tested by the simulation results of the COROT Seismic Working Group (Appourchaux, 2002, personal communication); this severely constraints the fundamental stellar parameters (mass, age, chemical composition) of models for which the frequencies of computed oscillation modes are similar to the observed high-order modes (Berthomieu et al., 2002). Moreover, solar-type modes also provide information on the star’s mean rotation rate.

A nearly regular spacing in frequency is also observed for opacity-driven low-frequency modes (Breger et al., 1999); the large separation of these low-frequency modes has to be similar between observations and theoretical models which satisfy also the properties of the observed high-frequency solar-type p modes. However, some of the opacity-driven modes deviate from the mean value of the large frequency separation; these modes are so-called mixed modes which provide details of the stellar core and of the precise evolutionary stage of the observed star (Unno et al., 1989, and references therein). This deviation from the mean value of the large frequency separation could suggest the presence of mixed modes. The problem is further complicated by the fact that the rotational splitting frequency components are no longer equidistant for these fast rotators, i.e. these frequencies could erroneously be identified as frequencies of mixed modes. However, knowing the mean rotation rate from the high-frequency splittings of solar-type p modes, the frequency splittings of the low-frequency opacity-driven modes can be computed in the manner of Dziembowski & Goode (1992) (see also Soufi et al., 1998).

The understanding of the physics responsible for the return to stability of opacity-driven modes at the red edge of the IS is still in its infancy. As the star becomes cooler the extent of the surface convection zone increases, thereby making the effect of convection-pulsation coupling on mode stability progressively more important. Several authors have tried to model the location of the red edge, e.g., Baker & Gough (1979), Bono et al. (1995) for RR Lyrae stars and e.g., Houdek (1996) and Xiong & Deng (2001) for $`\delta `$ Scuti stars. Although the authors assumed various models for the time-dependent treatment of convection, they all concluded that convection dynamics crucially effect the location of the red edge; however, different results were reported as to whether the convective heat flux (e.g., Bono et al., 1995), the momentum flux (e.g., Houdek, 1996) or turbulent viscosity (Xiong & Deng, 2001) is the crucial agent for stabilizing the modes at the red edge. In all these investigations, the predicted position of the red edge depends crucially on the assumed convection parameters, such as the mixing-length parameter or whether acoustic emission is included or neglected in the equilibrium model (Houdek, 2000).

Although it is possible from Fig. (13) of Houdek et al. (1999) to conclude that both types of modes can be excited simultaneously in the same star, amplitudes of stochastically excited modes for stars located inside the instability strip were not explicitly carried out by Houdek et al. (1999) and their possible detection were not addressed.

The aim of this paper is to demonstrate that models of stars, located inside the IS and near the red edge, can exhibit both opacity driven modes and solar-like oscillations with sufficiently large amplitudes to be detectable with today’s ground-based instruments. Consequently the planned asteroseismology space missions, such as COROT (COnvection ROtation and planetary Transits, Baglin & The Corot Team, 1998) or Eddington (Favata et al., 2000), will detect these oscillations even more easily.

Sect. 2 describes the equilibrium models, and the linear analysis results are discussed in Sect. 3., which are obtained from solving the equations of linear nonadiabatic oscillations in which convection is treated with the time-dependent, nonlocal formalism by Gough (1976, 1977, hereafter G’MLT). Furthermore, the effect of acoustic radiation in the equilibrium model on the stability properties is taken into account in the manner of Houdek (2000, and references therein). In this paper we consider only radial p modes.

Amplitudes of solar-like oscillations result from the balance between damping and stochastic driving by turbulence. The rate at which the turbulence injects energy into the p modes is estimated in the manner of Samadi & Goupil (2001, Paper I hereafter) and is discussed in Sect. 4.

In Sect. 5 we address the possibilities and conditions for detecting solar-type oscillations in $`\delta `$ Scuti stars with ground-based instruments and propose possible candidates, some of which are listed in the catalogue by Rodríguez et al. (2000). Conclusions are given in Sect. 6.

## 2 The stellar models

Equilibrium envelope models are computed in the manner of Houdek et al. (1999) using G’MLT formulation for convection. Integration starts at an optical depth of $`\tau =10^4`$ and ends at a radius fraction 0.2. Radiation is treated in the Eddington approximation and the atmosphere is assumed to be grey and plane parallel. In G’MLT formulation two more parameters, $`a`$ and $`b`$, are introduced which control the spatial coherence of the ensemble of eddies contributing to the total heat and momentum fluxes ($`a`$), and the degree to which the turbulent fluxes are coupled to the local stratification ($`b`$). In this paper we choose $`a^2=900`$ and $`b^2=2000`$ in order to obtain stable modes in the frequency range in which the damping rates exhibit a local minimum (e.g., at about 1.1 mHz for model C; see Section 3 and Fig. 2). The mixing-length parameter $`\alpha `$ has been calibrated to a solar model to obtain the helioseismically inferred depth of the convection zone of 0.287 of the solar radius (Christensen-Dalsgaard, Gough, & Thompson, 1991).

All models assume solar chemical composition and have mass $`M=1.68`$M and luminosity $`L=11.3`$L, but differ in effective temperature $`T_{\mathrm{eff}}`$, and whether or not acoustic radiation is included in the equilibrium computations. Table 1 lists the fundamental stellar parameters of these models. The models A1, A2 , B1 and B2 are hotter than model C and are located inside the IS and close to the red edge. Models A1 and A2 differ from models B1 and B2 by the inclusion of acoustic radiation by turbulence in the envelope calculations. In this model for acoustic radiation in the equilibrium model two more parameters are introduced (Houdek & Gough, 1998): the emissivity coefficient $`\mathrm{\Lambda }`$ and the parameter $`\mathrm{\Gamma }`$ which describes the power-law dependence of the acoustic power emission on the turbulent Mach number. A Mach-number dependence of $`\mathrm{\Gamma }=5`$ assumes that acoustic emission is dominated by the energy-bearing eddies ; if acoustic emission is predominantly emitted by inertial-range eddies $`\mathrm{\Gamma }`$ has the value 7.5. Table 2 lists the values of $`\mathrm{\Lambda }`$ and $`\mathrm{\Gamma }`$ that are assumed in the models A1 and A2. The values for $`\mathrm{\Lambda }`$ provide for a solar model a similar value for the acoustic flux $`F_{\mathrm{ac}}`$ as the estimates of Stein (1968) and Musielak et al. (1994). For all the models, except for model B2, we assume for the mixing-length parameter the calibrated solar value $`\alpha =2.037`$; for model B2 the value $`\alpha =1.5`$ is assumed.

Fig. 1 displays the locations of these models in the colour-magnitude diagram. Evolutionary tracks (dashed curves) are shown for models with various masses and are obtained with the CESAM code by Morel (1997) as described by Samadi et al. (2001a). The transformation from luminosity, effective temperature and surface gravity to absolute magnitude $`M_\mathrm{v}`$ and dereddened colour indices $`(by)_\mathrm{o}`$ are obtained from the Basel Stellar Library (Lejeune et al., 1998). The blue and red edges of the fundamental radial modes (solid curves) are calculated in the manner of Houdek et al. (1999).

The positions of the observed $`\delta `$ Scuti stars (filled circles) are taken from Rodríguez et al. (2000): absolute magnitudes, derived from Hipparcos distances and dereddened colour indices were kindly supplied by E. Rodríguez (2001, personal communication; see Rodríguez & Breger, 2001 , for details).

## 3 Stability analysis

The stability computations are as described in Houdek (2000, and references therein). In particular they include the Lagrangian perturbations of the turbulent fluxes (heat and momentum) according to Gough’s (1977) time-dependent formulation. Assuming a temporal dependence, $`\mathrm{exp}(\mathrm{i}\omega t)`$, for the pulsations, the complex eigenfrequencies of the modes can be written as $`\omega =\omega _\mathrm{r}+\mathrm{i}\omega _\mathrm{i}`$, which defines the cyclic pulsation frequency $`\nu =\omega _\mathrm{r}/2\pi `$ and the damping/growth rate $`\eta =\omega _\mathrm{i}/2\pi `$. The outer boundary conditions are applied at the temperature minimum, the mechanical boundary condition being consistent with a perfectly reflecting surface; at the base of the envelope, conditions of adiabaticity and vanishing displacement are imposed.

For model C all the modes are found to be linearly stable (i.e., $`\eta >0`$) as is expected for models lying well outside the IS. This is also found for the hotter model B1. For the model A1 (resp. A2) the first four (resp. three) radial modes, $`n`$=1,…,4 (resp. $`n`$=1,2,3), are found to be overstable. With the inclusion of a model for the acoustic radiation in the equilibrium structure the efficacy with which convection transports the turbulent fluxes is decreased (see Houdek & Gough, 1998). This leads to a decrease in the turbulent Mach number and to a consequent reduction of the stabilizing influence of the perturbed momentum flux on the mode damping. The driving eventually dominates over the damping leading to overstable modes.

Reducing $`\alpha `$ has a similar effect on mode stability than the inclusion of acoustic radiation in the equilibrium model (see Houdek & Gough, 1998; Michel et al., 1998; Houdek, 2000). The model B2 was computed with the smaller mixing-length parameter $`\alpha =1.5`$, leading to overstable modes with radial orders $`n=`$1,2,3.

Table 3 displays the frequency $`\nu `$ and damping/growth rate $`\eta `$ for all overstable radial modes ($`\eta <0`$) found in the models A1, A2 and B2.

Fig. 2 displays the damping rates as function of frequency for all stable modes and for all stellar models. The coolest model C exhibits a sharp dip in $`\eta `$ at about 1.1 mHz, which we moderated by applying a median filter with a width in frequency corresponding to nine radial modes; the result is plotted by the long-dashed curve in Fig. 2. This pronounced depression in $`\eta `$ in model C is related to the structure of the outer superadiabatic boundary layer: with decreasing surface temperature the location of the superadiabatic boundary layer is shifted progressively deeper into the star. This modifies the thermodynamic properties of this boundary layer of finite thickness, in particular, the thermal relaxation time (Balmforth, 1992). The thermodynamic coupling between the pulsations and the superadiabatic boundary layer becomes more efficient in cooler models, thereby promoting the depression in the damping rates by radiative processes (see Houdek et al., 1999).

## 4 Excitation rate and amplitude spectrum

The rms value of the mode surface velocity, $`v_\mathrm{s}`$, is related to the damping rate, $`\eta `$, and to the rate at which energy is injected into the mode (excitation rate), $`P`$, by

$$v_\mathrm{s}^2=\xi _\mathrm{r}^2(r_\mathrm{s})\frac{P}{2\eta I},$$
(1)
where $`\xi _\mathrm{r}`$ is the radial displacement eigenfunction, $`r_\mathrm{s}`$ is the radius at which the surface velocities are measured and which we assume to be 200 km above $`T_{\mathrm{eff}}`$, and the mode inertia $`I`$ satisfies

$$I=_0^M\xi ^{}\xi dm.$$
(2)
The rate of energy injected into a mode is computed according to Paper I and is proportional to

$`P(\omega )`$ $``$ $`{\displaystyle _0^M}\rho w^3\mathrm{}^4\left({\displaystyle \frac{\mathrm{d}\xi _\mathrm{r}}{\mathrm{d}r}}\right)^2𝒮(\omega ,m)dm,`$ (3)
where $`\rho `$ is the density, $`\mathrm{}`$ is the mixing length, and $`w`$ is the vertical component of the rms velocity of the convective elements. The function $`𝒮(\omega ,m)`$ describes approximately contributions from eddies with different sizes to the excitation rate $`P`$. Detailed expressions for $`𝒮(\omega ,m)`$ were given in Paper I. For estimating $`𝒮`$, assumptions for the turbulent kinetic energy spectrum $`E(k)`$, and for the turbulent spectrum of the entropy fluctuations $`E_\mathrm{s}(k)`$ have to be made, where $`k`$ is the eddy wavenumber. In this paper we assume for $`E(k)`$ the “Nesis Kolmogorov Spectrum” (NKS hereafter) as discussed in Paper I. This turbulent spectrum is obtained from observations of the solar granulation by Nesis et al. (1993), and leads to the best agreement between a solar model using our stochastic excitation theory and solar measurements (Samadi et al., 2001b).

Results for the estimated excitation rate $`P`$ are depicted in Fig. 3. For the models A1, A2 and B1 the excitation rate $`P`$ is about one magnitude larger than for model C. This is a result of the larger convective velocities in the superadiabatic boundary layers of the models A1, A2 and B1 (see Fig. 4), which are all hotter than model C. The frequency dependence of $`P`$ for the models A1 and B1 are similar as it is also the case for the models A2 and B2; the latter two models, however, are quite different from the first two models. This difference in the frequency dependence of $`P`$ is a consequence of the different profiles of the convective velocities $`w`$ (see Fig. 4); the profiles of $`w`$ are similar for A2 and B2 but differ substantially from A1 and B1. For the models A2 and B2 the efficacy of convection has been reduced severely by either including acoustic radiation in the equilibrium model (A2) or by reducing the mixing-length parameter $`\alpha `$ to a value much smaller than the calibrated value for a solar model (B2). This results in shallower superadiabatic regions and in larger superadiabatic temperature gradients; pulsation modes in A2 and B2 are therefore predominantly excited at the very top of the convection zone, whereas in the models A1 and B1 the modes are excited over a larger driving region. The two sets of values for $`\mathrm{\Lambda }`$ and $`\mathrm{\Gamma }`$, listed in Table 2, provide in a solar model approximately the same value for the acoustic flux $`F_{\mathrm{ac}}`$ (see Section 2). In model A2 (which assumes $`\mathrm{\Lambda }=2000`$ and $`\mathrm{\Gamma }=7.5`$), however, $`F_{\mathrm{ac}}`$ is about three times larger than in model A1 (which assumes $`\mathrm{\Lambda }=100`$ and $`\mathrm{\Gamma }=5`$), and the associated velocities, plotted in Fig. 4, are correspondingly greater.

It is interesting to note that in Fig. 4 the convective velocities in the superadiabatic boundary layers of models A1 and A2 are somewhat larger than the others. Although including acoustic radiation in the mean stellar structure leads to a reduction of the maximum turbulent Mach number $`M_\mathrm{t}w/c`$ (by $`1.5\%`$ for model A1 and $`5\%`$ for model A2 relative to B1), the whole structure of the outer layers changes too, thereby increasing the locally defined adiabatic sound speed $`c`$. This leads to an increase of the convective velocities $`w`$ in the very outer layers despite the decrease in $`M_\mathrm{t}`$.

The less smooth frequency dependence (wiggles) of $`P`$ for $`\nu 0.7\mu `$Hz for model B1 (and model C) is related to the location and extent of the driving region: the radial eigenfunctions $`\xi _\mathrm{r}`$ vary rapidly with depth and frequency (particularly for the high-frequency modes). As discussed above, in the models A2 and B2 the modes are predominantly excited in a shallow region beneath the surface, where the expression $`(\mathrm{d}\xi _\mathrm{r}/\mathrm{d}r)^2`$ (see Eq. 3) varies monotonically with frequency (Goldreich et al., 1994; Samadi et al., 2002), leading to the smooth frequency dependence of $`P`$ for $`\nu 0.7\mu `$Hz, as depicted in Fig. 3. The larger driving regions in the models B1 and C extend to layers where the expression $`(\mathrm{d}\xi _\mathrm{r}/\mathrm{d}r)^2`$ no longer varies monotonically with frequency leading to the frequency-dependence of $`P`$ as shown by the dot-dashed and long-dashed curves in Fig. 3.

In the top panel of Fig. 5 the surface velocity amplitudes $`v_\mathrm{s}`$ are depicted for all stellar models, computed according to Eq.(1). In the models A1, A2, B1 and B2 the amplitudes of stochastically excited p modes are larger ($`59\mathrm{ms}^1`$) than in model C ($`2\mathrm{ms}^1`$).

For estimating the luminosity amplitudes the full nonadiabatic luminosity eigenfunctions have to be used. The relative luminosity amplitudes, $`\delta L/L`$, are linearly related to the velocity amplitudes, i.e. they are proportional to the ratio of the luminosity eigenfunction over the displacement eigenfunction. This ratio is determined by the solution of the nonadiabatic pulsation equations and is independent of a stochastic excitation model (see Houdek et al., 1999). In the middle panel of Fig. 5 the amplitude ratios, $`\mathrm{\Delta }L/\mathrm{\Delta }v_\mathrm{s}`$, are plotted as a function of frequency for all stellar models. The shape of the amplitude ratios are in general similar between all the models with the smallest ratios predicted for the models A2 and A1. Only at the highest frequencies the amplitude ratios are considerably larger in A2 and B2; at high frequencies nonadiabatic effects due to radiative dissipation in the radiative zone, below the shallow surface convection zones in A2 and B2, lead to an increase in the amplitude of the luminosity eigenfunctions and consequently in the luminosity amplitudes. The velocity amplitudes in Fig. 5 are obtained 200 km above the photosphere ($`T=T_{\mathrm{eff}}`$) and do increase by a factor of about two at the outermost meshpoint of the model, i.e. at an optical depth $`\tau =10^4`$.

We predict a maximum value of the luminosity amplitude $`\delta L/L97`$ ppm for model A1, $`\delta L/L150`$ ppm for model A2, $`\delta L/L101`$ ppm for model B1 , $`\delta L/L98`$ ppm for model B2 and $`\delta L/L84`$ ppm for model C. These results are summarized in Table 4.

The dotted horizontal line in the middle panel of Fig. 5 represents an order-of-magnitude estimate of the amplitude ratio according to Kjeldsen & Bedding (1995):

$$\delta L/Lv_\mathrm{s}T_{\mathrm{eff}}^{1/2}.$$
(4)
Kjeldsen & Bedding derived this expression for a purely radiative model assuming simplified proportional relations in the adiabatic approximation. This simplified scaling law suggests smaller values for the amplitude ratios and consequently leads to smaller luminosity amplitudes $`\delta L/L`$, particularly at high frequencies, where nonadiabatic effects are important. At a frequency $`\nu 1`$ mHz, for example, the scaling law (4) predicts for model A1 a luminosity amplitude which is about three times smaller than that obtained from the nonadiabatic computation.

There is evidence that energy equipartition holds for the Sun (apparently fortuitously); an estimate of the total energy in the modes is, however, only possible for the Sun for which accurate data are available; using GONG data the total energy for modes with degrees $`l=0,\mathrm{},300`$ and with radial order $`n`$ up to 10 is found to be approximately $`2\times 10^6E_0`$, where $`E_02\times 10^{28}`$ erg is the maximum value of the kinetic energy in a particular ridge (i.e. for a particular value of $`n`$) and which is independent of $`l`$ (see Fig. 5 of Komm et al., 2000). The value $`2\times 10^6`$ is also roughly equal to the number of granules on the solar surface, a result which supports the energy equipartition principle. In other stars, however, energy equipartition does not necessarily hold, because we have a nonequilibrium dynamical (yet statistically steady) system in which the damping and excitation is balanced in a nonlinear way by the energy input and output, i.e. it is not determined by equilibrium. In such a nonequilibrium situation there is no general physical principle limiting the ratio of the energy in the oscillation mode to the energy in the convection. Another facet of such a dynamical process is provided by the reaction of convection to the acoustical radiation; although the latter contributes towards augmenting the damping of an eddy, the resultant change of the background stratification is such as to augment the driving by even more, causing the convective velocities to increase.

It is perhaps interesting to mention that the acoustic energy flux generated by the fluctuating Reynolds stress of the turbulent velocity field is relatively small compared to the total energy flux carried by the convection; the ratio between the acoustic energy flux emitted by the energy-bearing eddies and the convective energy flux is proportional to $`M_\mathrm{t}^\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }=5`$ for homogeneous, isotropic turbulence (see Lighthill, 1952). The turbulent Mach number in $`\delta `$ Scuti stars is in general much smaller than unity and consequently this ratio is small. For the Sun this ratio is of the order of $`10^3`$ (see e.g., Stein, 1968). The acoustic flux emitted predominantly by inertial-range eddies is proportional to $`M_\mathrm{t}^\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }=15/2`$ (see Goldreich & Kumar, 1990), i.e. it scales with an even higher power of the Mach number (see also Section 2). Consequently the total amount of acoustic energy injected into the p modes is small compared to the energy carried by the convection. In a fully convective envelope the total energy flux (luminosity) is carried solely by the turbulent velocity field, i.e. in that case the luminosity is a measure of the total energy in the convection. Therefore the ratio between the energy supply rate for a particular mode and the luminosity, $`P/L`$, is proportional to the ratio between the energy in that mode and the total energy in the convection. In the Sun this ratio is of the order of $`10^{11}`$ for the mode with the largest amplitude. In model A2 this ratio is $`10^8`$, which is still small.

## 5 Observational constraints for detecting solar-type oscillations

There have been recent reports on the possible detection of solar-type oscillations in $`\alpha `$ Cen (HD 128620) by Bouchy & Carrier (2001), in $`\beta `$ Hydri (HD 2151) by Bedding et al. (2001) and in Procyon A (HD 61421) by Martic et al. (1999, see also , ), who obtained spectroscopic surface velocity measurements of these bright stars (the apparent magnitude $`V=2.80`$ for $`\beta `$ Hydri, $`V=0.34`$ for Procyon and $`V=0.1`$ for $`\alpha `$ Cen) from the ground. The maximum values of the observed peak-velocity amplitudes are of the order $`35\mathrm{cms}^1`$ for $`\alpha `$ Cen, $`50\mathrm{cms}^1`$ for $`\beta `$ Hydri and $`50\mathrm{cms}^1`$ for Procyon. Current ground-based instruments are able to detect oscillations with velocity amplitudes of the order predicted for our models A1, A2 and B1, B2, but only for stars with an apparent magnitude $`V`$ of less than $`34`$ (Bouchy, personal communication). The HARPS (High-Accuracy Radial-velocity Planetary Search) project (Bouchy & Carrier, 2001), for example, will be able to detect oscillations with our predicted velocity amplitudes for stars with an apparent magnitude smaller than $`45`$. This detection threshold is still too small for detecting solar-type oscillations in $`\delta `$ Scuti stars located near the red edge of the IS, particularly in view of the fact that most of the currently known $`\delta `$ Scuti stars are even fainter. For example, the apparent magnitudes of known $`\delta `$ Scuti stars located nearest to the red edge (see Fig. 1) are between $`V=5.7`$ and $`V=9.3`$.

Future space missions with instruments dedicated to asteroseismology, however, will be able to detect solar-like oscillations in $`\delta `$ Scuti stars: the forthcoming space project COROT (Baglin & The Corot Team, 1998), for example, will reach a noise level of $`0.7`$ ppm (Auvergne & the COROT Team, 2000) for a star with an apparent magnitude of $`V=6`$, using photometric measurements. Therefore, in stars with similar magnitudes, COROT will be able to detect oscillation amplitudes as small as $`3`$ ppm, a value which is similar to that measured in the Sun. The instrument on COROT will be limited by the photon noise only for stars with magnitudes larger than $`V9`$: i.e., for a star with magnitude $`V8`$ the detection threshold will be $`5`$ ppm. This threshold is small enough to detect and measure many solar-like oscillations in $`\delta `$ Scuti stars which are similar to the $`\delta `$ Scuti models considered in this paper.

## 6 Conclusion

We studied oscillation properties in $`\delta `$ Scuti stars located near the observed red edge of the classical instability strip. Such stars can pulsate with both opacity-driven modes and intrinsically stable stochastically driven (solar-like) modes. The estimated velocity amplitudes of the stochastically driven modes in our $`\delta `$ Scuti models are found to be larger than in cooler and pulsationally stable models lying outside the IS. This result supports the idea that solar-like oscillations in $`\delta `$ Scuti stars may be detected.

Including a model for the acoustic radiation in the equilibrium model results in a cooler red edge and does effect the properties of the excitation rate of p modes (see also Houdek & Gough, 1998, Houdek, 2000); in particular the pulsation amplitudes do become larger and are predicted to be largest for a model with the largest acoustic flux $`F_{\mathrm{ac}}`$ (i.e., model A2). Moreover, for the $`\delta `$ Scuti models considered in this paper, overstable modes were predicted only if either acoustic emission in the mean stratification was included or if the mixing-length parameter was reduced to a value smaller than suggested by a calibrated solar model.

A potential target star should neither be too cool (i.e., no opacity-driven modes) nor too hot (i.e., stochastically excited modes with amplitudes too small to be detectable). We quantify this with the illustrative case of our $`\delta `$ Scuti models with a mass $`M=1.68`$M and we identify the following $`\delta `$ Scuti stars from the Rodríguez et al. (2000) catalogue, located near the red edge, as potential candidates for the target selection of upcoming observing campaigns: HD57167, HD14147, HD208999 and HD105513.

Although the amplitudes of the solar-type oscillations, predicted in our $`\delta `$ Scuti models, are large enough to be detected from ground, today’s ground-based instruments will detect such oscillations only in brighter $`\delta `$ Scuti stars with an apparent magnitude of up to $`V34`$ (Bouchy, 2001, personal communication). However, new ground-based observing campaigns, such as the HARPS project (Bouchy & Carrier, 2001) will be able to detect stochastically excited oscillations in $`\delta `$ Scuti stars with an apparent magnitude of up to $`V45`$. Unfortunately, there are no such bright stars in the Rodríguez et al. (2000) catalogue which are located near the red edge, although some bright stars near the red edge may have opacity-driven modes with amplitudes too small to be detectable with today’s ground-based instruments and are therefore not classified as $`\delta `$ Scuti stars.

The forthcoming space missions for asteroseismology, such as COROT (Baglin & The Corot Team, 1998) and Eddington (Favata et al., 2000) will be able to detect solar-like oscillations in faint $`\delta `$ Scuti stars. The large instrument on the Eddington spacecraft will measure stellar oscillations with amplitudes as small as $`1.5`$ ppm in stars with an apparent magnitude of $`V11`$ assuming an observing period of 30 days. Moreover, Eddington’s large field of view will allow it to monitor a large number of stars simultaneously. This will be helpful for detecting and classifying new $`\delta `$ Scuti stars and for measuring the location of the red edge of the IS with greater precision than it was possible before.

###### Acknowledgements.

We thank E. Rodríguez for providing the $`\delta `$ Scuti data set in a convenient and immediate usable form, T. Lejeune for allowing us to use the Basel library and D. Cordier for providing it on the Internet. We thank A. Baglin for useful discussions on the COROT specifications, F. Bouchy for providing valuable information on the HARPS project and related experiments, and C. Catala and E. Michel for useful discussions on the possibilities of detecting new $`\delta `$ Scuti stars. We are grateful to Douglas Gough for very helpful discussions on stochastic mode excitation and to Mike Montgomery for improving the English. GH and RS acknowledge support by the Particle Physics and Astronomy Research Council of the UK. RS’s work has been supported under the grant PPA/G/O/1998/00576.


text:
# Limits on the Short Term Variability of Sagittarius A* in the Near-Infrared

## 1 Introduction

The variability of Sagittarius (Sgr) A\* at X-ray wavelengths (Baganoff et al., 2001a) has bolstered the case for associating this source with the suspected 2.6 x 106 M black hole at the center of our Galaxy (Eckart & Genzel, 1997; Genzel et al., 1997, 2000; Ghez et al., 1998, 2000). In two Chandra observations separated by almost a year and having a total of 76 ksec of exposure time, Sgr A\* was detected at X-ray wavelengths for the first time and was also seen to flare in intensity over a time scale of 3 hours (Baganoff et al., 2001a, b). While the flare’s short duration implied a small region of origin, $``$400 Rs (where $`R_s`$ is the Schwarzschild radius = 2GM/c2), its large amplitude, a factor of 50, has raised the possibility of detecting corresponding intensity enhancements at wavelengths outside the X-ray regime.

Existing models for Sgr A\*’s flared state make very disparate predictions for the emission at wavelengths between the X-ray and radio regimes (Markoff et al., 2001; Liu & Melia, 2002). The wide differences between these models are a result of assuming different geometries (disk vs. jet) and emission mechanisms for the flaring process (e.g., enhanced accretion rates vs. magnetic reconnection). In some models, the predicted emission in the flared state, at infrared (IR) wavelengths, dramatically exceeds that of existing detection limits (Genzel & Eckart, 1999; Stolovy et al., 1999; Morris et al., 2001). For example, the preferred model of Markoff et al. (2001) predicts an observed 2µm flux density of $``$13 mJy, or a dereddened flux density of $``$300 mJy, during the flared state. Unlike the situation at radio wavelengths, where long-term monitoring campaigns have been used to constrain the flared state of Sgr A\* (Zhao, Bower, & Goss, 2001), the limited time coverage and spatial resolution of published IR experiments prevent meaningful constraints on the flared state’s IR emission from being inferred111While several papers have reported the possible detection of a variable near-infrared source coincident with Sgr A\* (Herbst, Beckwith, & Shure, 1993; Close, McCarthy, & Melia, 1995; Genzel et al., 1997), subsequent high resolution observations have identified this emission to be from high proper motion sources (Eckart et al., 1995; Eckart & Genzel, 1997; Ghez et al., 1998, 2002). and, thus, the reported limits are assumed to be associated with Sgr A\*’s quiescent state.

The W. M. Keck Observatory dynamical study of stars in the central stellar cluster (Ghez et al., 1998, 2000; Gezari et al., 2002) provides a rich source of high angular resolution 2µm data between 1995 and 2001. In this paper, we present 2µm flux density limits from maps that were each composed of data from a single night. The elapsed time of 3-4 hours in each map is approximately the same as the time scale of the observed X-ray flare, making this data set ideally suited for possibly detecting a flare of this type. Given the quantity of these observations, our non-detections establish a robust upper limit for the flared state’s 2µm emission intensity.

## 2 Observations

High resolution, near-infrared observations of the Galactic Center were conducted from 1995 June to 2001 July using both speckle and adaptive optics (AO) imaging techniques on the Keck 10-m telescopes. The speckle observations were obtained in the K-band ($`\lambda _o`$ = 2.2µm, $`\mathrm{\Delta }\lambda `$=0.4µm) using the Keck I facility near-infrared camera (NIRC; Matthews & Soifer 1994) with external re-imaging optics. This resulted in a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$0203 and a field of view (FOV) of 5$`\stackrel{}{\mathrm{.}}`$12$`\times `$5$`\stackrel{}{\mathrm{.}}`$12 (Matthews et al., 1996). During each night of observations several thousand short exposures (texp= 0.137 sec) were taken in sets of $``$200. A more limited set of data was collected using two different science cameras behind the Keck II AO system (Wizinowich, 2000). The first AO data set was collected in the K’-band ($`\lambda _o`$ = 2.1µm, $`\mathrm{\Delta }\lambda `$=0.35µm) in early 1999 with the near-infrared engineering camera (KCAM; Wizinowich et al. 2000), which had a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$0175, a FOV of 4$`\stackrel{}{\mathrm{.}}`$4$`\times `$4$`\stackrel{}{\mathrm{.}}`$4. Each image had an exposure time of 5 sec. The slit-viewing camera of NIRSPEC (SCAM; McLean et al. 1998) provided a second set of AO images for this study. These images, like the speckle images, were made in the K-band and have a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$0170 and a FOV of 4$`\stackrel{}{\mathrm{.}}`$4$`\times `$4$`\stackrel{}{\mathrm{.}}`$4, and exposure time of 10 sec. USNO 0600-28579500 served as the natural guide star for all of these AO observations. Since this guide star is both faint (R = 13.2) and distant from the target (r$``$30″), the AO performance was non-optimal. Table 1 provides a summary of all observations.

## 3 Data Analysis & Results

Three basic steps constitute the data analysis process in this program. First, high angular resolution maps are generated from the individual short exposure frames (§3.1). Second, all point sources in the FOV are identified and a direct detection of Sgr A\* is ruled out (§3.2). Third, limits for Sgr A\* are derived from the residual maps, in which all identified point sources have been removed (§3.3).

### 3.1 Construction of Images

Image processing proceeds similarly to that carried out for the dynamical experiment with one exception. Rather than combining all the data from the duration of an observing run, typically 2-3 nights, we synthesize the data over each night to produce 27 maps, each of which is limited to an elapsed time of 3-4 hours. Since the details of this method are described elsewhere (Ghez et al., 1998), only a brief summary is provided here. Standard image reduction techniques are applied to all the individual speckle and AO frames. For the speckle data, a two stage shift-and-add (SAA; Christou 1991; Ghez et al. 1998) analysis then produces the final high resolution maps. In the first stage, the 200 frames in each set are combined to form an intermediate SAA image. Then, these multiple intermediate SAA images (from throughout the night) are combined to form one final SAA map. This allows each intermediate image to be examined for seeing quality. In combining the intermediate images, a seeing cut is established so as to exclude those images with the worst seeing from the final map. For the AO data, this cut is also carried out on the individual AO images before they are registered and averaged together. Figure Limits on the Short Term Variability of Sagittarius A\* in the Near-Infrared displays representative final nightly speckle and AO maps.

### 3.2 Point Source Identification & Search for Sgr A\*’s near-infrared emission

In all maps, stars are identified using StarFinder, an IDL package developed for astrometry and photometry in crowded stellar fields (Diolaiti et al., 2000). This package iteratively generates estimates of the point spread function (PSF) from a few selected bright stars and then identifies point sources over the entire FOV through cross-correlation of the map with the PSF model. For the PSF extraction, we found that the most reliable PSF models are obtained with a support size of $``$2″, which represents a compromise between needing to accommodate the large PSF halos and yet having a limited FOV. This choice limits our analysis to images with PSF halo sizes of 0$`\stackrel{}{\mathrm{.}}`$4 or less, as the PSFs of the remaining two images are poorly characterized by this process. The PSF model is based on four of the five brightest stars in the FOV (IRS 16NE, 16C, 16NW, and 29N; see Figure Limits on the Short Term Variability of Sagittarius A\* in the Near-Infrared); IRS 16SW is avoided as a PSF model star, since it is surrounded by relatively bright stars in its immediate vicinity. For point source identification, a correlation coefficient greater than 0.8 between the PSF model and the actual stellar image is required to avoid spurious detections. This process results in the identification of $``$100 point sources in each map.

The speckle and AO images have significantly different PSFs. Nonetheless, both PSFs are composed of a compact core on top of a broader halo. Table 1 provides the characterization of the PSF in each map based on the radial profile of the PSF model. While the speckle images have PSF core FWHM that are nearly diffraction limited ($``$0$`\stackrel{}{\mathrm{.}}`$05) and $``$40% smaller than that of the AO images ($``$0$`\stackrel{}{\mathrm{.}}`$08), the AO PSF core contains $``$30% of the total energy, $``$12 times more than the typical speckle PSF.

At this stage in the analysis, it is possible to look for direct detections of Sgr A\*. We use proper motion acceleration vectors (Ghez et al., 2000) to pinpoint the central black hole’s position relative to the nominal radio position of Sgr A\* (Menten et al., 1997) to within 0$`\stackrel{}{\mathrm{.}}`$04 (1$`\sigma `$). Within 0$`\stackrel{}{\mathrm{.}}`$08 of this location, 4 sources (13.9 $`<`$ K $`<`$ 16.5) are identified, all of which were previously detected in ”monthly” maps made from all data in a single observing run (Ghez et al., 1998, 2002) and, furthermore, have significant proper motions. This high stellar density emphasizes the need for improved accuracy in Sgr A\*’s position in the IR reference frame in order to measure or constrain its emission. With no stationary source identified in this region, we conclude that Sgr A\* has not been detected.

### 3.3 Flux Density Limits for Sgr A\*

In order to determine an accurate detection limit at the position of Sgr A\*, it is necessary to remove the contaminating seeing halos from nearby sources. A ‘stars only’ map is created using the PSF model and list of stars generated by StarFinder. This is then subtracted from the original map, producing a residual map. With the residual map, a 3$`\sigma `$ point source detection limit for Sgr A\* is established based on three times the RMS of 25 aperture photometry values, which are calculated using $``$60 mas radius apertures and sky annuli extending from $``$60 mas to $``$90 mas. The 5$`\times `$5 grid of apertures in the residual map corresponds to an area of $``$0$`\stackrel{}{\mathrm{.}}`$6$`\times `$0$`\stackrel{}{\mathrm{.}}`$6, more than two orders of magnitude larger than the uncertainty in the location of Sgr A\*.

Zero points are obtained through carrying out the same aperture photometry in the original maps (prior to the ‘stars only’ subtraction) on all known non-variable sources brighter than K=10.5,using the flux densities reported in Blum, Sellgren, & Depoy (1996), and that occur in more than 30% of the frames for each night. The photometric calibration sources used are IRS 16NW, 16C, 16CC, and 16NE, when the FOV allows its inclusion; IRS 29N is omitted as it is found to be marginally variable at the 2$`\sigma `$ level. Typical photometric zero point 1$`\sigma `$ uncertainties of $``$0.04 mag result from this procedure.

Table 1 and Figure Limits on the Short Term Variability of Sagittarius A\* in the Near-Infrared contain the resulting 3$`\sigma `$ point source detection limits for Sgr A\*. The lowest of these upper limits gives an observed limit for the quiescent state of Sgr A\* of 0.09$`\pm `$0.005 mJy, or, equivalently, a dereddened value of 2.0$`\pm `$0.1 mJy, while the highest limit constrains our analysis of the variable state of Sgr A\* to $``$0.8 mJy (observed) or 19 mJy (dereddened). Although these upper limits vary significantly from epoch to epoch (due to variable observing conditions and/or a variation in the number of contaminating sources detected and removed) the lowest of them is lower than the best previously reported limits at 2µm for Sgr A\*’s quiescent state (dereddened 4 mJy; Genzel & Eckart 1999) by a factor of 2.

## 4 Discussion

Using the coverage of the X-ray and near-infrared experiments and assuming random 3-hr flaring events that produce both enhanced X-ray and near-infrared emission, we consider the likelihood that a near-infrared flare in excess of our weakest limit occurred over the course of the IR experiment. The probability of seeing such a flare is given by the following binomial distribution:

$$P(\nu \text{ detections in }n\text{ trials})=\frac{n!}{\nu !(n\nu )!}p^\nu q^{n\nu }$$
(1)
where $`p`$ is the probability of detecting a flare in one trial and $`q`$ is the probability of no detection (1-$`p`$). Here we will refer to $`p`$ as the duty cycle and note that in our context, this quantity describes what fraction of a set of 3-hr images should contain a flare. Figure Limits on the Short Term Variability of Sagittarius A\* in the Near-Infrared shows the probability distributions as a function of underlying duty cycle for both the X-ray experiment, in which $`\nu `$=1 and $`n`$=7.6 (obtained by dividing the total time, 76 ksec, by the length of the observed X-ray flare, 10 ksec) and the near-infrared experiment, in which $`\nu `$=0 and n=25. The X-ray and near-infrared probability distributions have only a small overlap, which suggests that if the flaring activity at X-ray and near-infrared wavelengths are coupled, a flare most likely occurred during the near-infrared experiment. The joint probability distribution quantifies this and suggests there is, at most, a probability of 9% of a random 3-hr near-infrared flare in excess of 19 mJy (dereddened). We therefore assume that at the 2$`\sigma `$ confidence level, a flare occurred during our experiment and use our limits to constrain the variability models.

While only a limited amount of modeling of the recent 3-hr X-ray flare detected at Sgr A\* has been carried out, existing models predict 2µm dereddened emission as high as 300 mJy in the model preferred by Markoff et al. (2001) but as small as 0.4 mJy in Liu & Melia (2002). The former model explains the elevated X-ray emission, produced by synchrotron self-Compton, by an enhanced temperature for the relativistic electron population, as might arise in a magnetic reconnection event. On the other hand, Liu & Melia present a flare model in which the flare arises due to an enhanced accretion rate and bremsstrahlung emission is dominant at both near-infrared and X-ray wavelengths. The lack of a near-infrared detection of Sgr A\* makes the Markoff et al. (2001) model and any other mechanism that produces flared 2µm emission in excess of 19 mJy (dereddened) unlikely.

## 5 Conclusions

This paper summarizes a search for a near-infrared counterpart to Sgr A\* in the flared state. From the length of our observations, this search was sensitive to variability on time scales of 3 hours. No such counterpart was detected. However, by identifying and removing all the stars in the crowded inner $``$0$`\stackrel{}{\mathrm{.}}`$6$`\times `$0$`\stackrel{}{\mathrm{.}}`$6 of the Galactic Center, an upper limit for the emission from Sgr A\* has been inferred for each observation epoch. These limits constrain the quiescent emission from Sgr A\* to $``$0.09 mJy (2.0 mJy, dereddened) and the variable component to $``$0.8 mJy (19 mJy, dereddened) at the 2$`\sigma `$ confidence level. More X-ray data will improve the estimated flaring duty cycle and is likely to increase the confidence level of our limits for the near-infrared component of the flared emission.

Valuable input on the manuscript was contributed by Gaspard Duchêne, Caer McCabe, and Andreas Eckart. Data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This work was supported by the NSF both through grant No. 9988397 and the Science and Technology Center for Adaptive Optics, managed by the University of California at Santa Cruz under cooperative agreement No. AST - 9876783. AMG also thanks the Packard Foundation for its support.


text:
# Conditional phase shifts using trapped atoms

## Abstract

We describe a scheme for producing conditional nonlinear phase shifts on two-photon optical fields using an interaction with one or more ancilla two-level atomic systems. The conditional field state transformations are induced by using high efficiency fluorescence shelving measurements on the atomic ancilla. The scheme can be nearly deterministic and is of obvious benefit for quantum information applications.

It has recently been shown KLM  that nonlinear phase shifts on two photon states can be produced by coupling the mode of interest to ancilla modes via a beam splitter and making photon counting measurements on the ancilla modes. Such conditional nonlinear phase shifts can be used to perform two qubit operations for logical states encoded in photon number states. If such conditional gates are used to prepare entangled states for teleportation, efficient quantum computation can be performed which with suitable error correcting codes can be made fault tolerant KLM . In this paper we show that if the ancilla modes are replaced with a two level atom similar conditional nonlinear phase shifts can be achieved by near deterministic postslection on atomic measurements. The atomic measurements can made with fluorescence shelving techniques which are very much more efficient than single photon counting measurements, thus reducing the need for new photon counting technologies inherent in the KLM scheme. Recently, a different scheme for conditional quantum gates based on atomic systems was presneted by Protsenko et al.protsenko

Consider a single optical mode prepared in an arbitrary two photon state

$$|\psi =c_0|0+c_1|1+c_2|2$$
(1)
Our objective is to find a way to produce the nonlinear phase shift transformation defined by

$$|\psi c_0|0+c_1|1c_2|2$$
(2)
To achieve this result we assume that at some fixed time an interaction between the field mode and a single two level atom is switched on. After some interaction time $`t`$ the interaction is turned off and the atomic state is measured by fluorescence shelving. The resulting conditional state of the field will then depend on the initial state of the two-level atom and the interaction time. We will show that these can be so arranged as to effect the nonlinear phase shift required.

We have in mind a quantum computing communication protocol in which the optical field mode is derived from a transform limited pulsed field which is rapidly switched into the cavity mode containing the atomic systems at fixed times determined by the pulse repetition rate. Similar systems have been proposed as a quantum memory for optical information processing Pittman2002 . When the atomic measurement yields the required result the field may be switched out again for further analysis or subsequent processing through linear and conditional elements.

Once the cavity field is prepared, we need to switch on the interaction with the atomic system. In order that we can switch this interaction at predetermined times we propose that an effective two level transition connected by a Raman process with one classical field and the quantised signal field, be used. A similar scheme has recently been proposed as the basis of a high efficiency photon counting measurement Imamoglu ; James\_Kwiat . The process is also used in the EIT schemes for storing photonic information Lukin2000  and for quantum state transfer between distant cavities Cirac1997 . The level diagram is shown in figure 1. The nearly degenerate levels $`|1`$ and $`|2`$ are connected by a stimulated Raman transition to level $`|3`$. The detuning of the Raman pulse from the excited state $`|3`$ is $`\mathrm{\Delta }`$, which is approximately the same as the detuning of the signal mode form the same transition. An advantage of using a stimulated Raman process of this kind is that the excited state $`|2`$ can be a metastable, long lived level. We thus do not need to consider spontaneous emission from this level back to the ground state. The readout of the atomic system may be achieved by using a cycling transition between the excited state $`|2`$ and another probe level $`|4`$. Such measurements are routinely performed in ion trap studies Rowe2000  and can have efficiencies greater than 99%.

The interaction between the single mode field and a two-level atom is described by the effective Hamiltonian

$$\widehat{𝑯}=\kappa (\widehat{𝒂}^{}\widehat{𝝈}^{}+\widehat{𝒂}\widehat{𝝈}^+)$$
(3)
The interaction strength is given by $`\kappa =\mathrm{\Omega }g/2\mathrm{\Delta }`$ where $`\mathrm{\Omega }`$ is the Rabi frequency for the Raman pulse and $`g`$ is the one photon Rabi frequency for the signal field. The unitary transformation that acts when this interaction is applied for a time $`t`$ is

$$\widehat{𝑼}(\tau )=\mathrm{exp}[i\tau (\widehat{𝒂}^{}\widehat{𝝈}^{}+\widehat{𝒂}\widehat{𝝈}^+)]$$
(4)
where $`\tau =\kappa t`$.

If the atom is prepared in the ground state and found in the ground state after the interaction, the conditional state of the field is given by

$$\widehat{𝚼}_{gg}(\tau )|\psi =\mathrm{cos}(\tau \sqrt{\widehat{𝒂}^{}\widehat{𝒂}})|\psi $$
(5)
On the other hand if the atom is prepared in the exited state and found in the excited state after the interaction, the conditional state is given by

$$\widehat{𝚼}_{ee}(\tau )|\psi =\mathrm{cos}(\tau \sqrt{\widehat{𝒂}\widehat{𝒂}^{}})|\psi $$
(6)
There is considerable practical advantage to using the excited state preparation rather than the ground state as there is always a signal for correct operation, however the analysis is the same.

Now suppose we send in a generic two photon state

$$|\psi =c_0|0+c_1|1+c_2|2$$
(7)
which interacts with a single two level atom, prepared in the ground state, and after the interaction the atom is found still to be in the ground state. In this case we need to apply the measurement operator $`\widehat{𝚼}_{gg}`$, and the resulting state of the field is

$$|\varphi _g=A_0c_0|0+A_1c_1|1+A_2c_2|2$$
(8)
where $`A_0=1`$, $`A_1=\mathrm{cos}(\tau )`$ and $`A_2=\mathrm{cos}(\sqrt{2}\tau )`$. Note that the frequencies of the $`A_1`$ and $`A_2`$ terms are irrational multiples of each other. As we sweep $`\tau `$, $`A_1`$ and $`A_2`$ should explore their entire phase space. It should be possible to find values of $`\tau `$ for which $`(A_1,A_2)`$ approaches arbitrarily close to $`(1,1)`$. These solutions can be found trivially for small $`\tau `$ by plotting $`A_1`$ and $`A_2`$ and visually inspecting for intersections near 1. Some high-probability solutions are summarised in table 1.

With two atoms, one initially prepared in the ground state and the second prepared in the excited state, the conditional state given that both atoms are found in their initial state after the interaction is

$`|\varphi _{ge}`$ $`=`$ $`\widehat{𝚼}_{gg}(\tau _1)\widehat{𝚼}_{ee}(\tau _2)|\psi `$ (9)
$`=`$ $`A_0c_0|0+A_1c_1|1+A_2c_2|2`$ (10)
where $`A_0=\mathrm{cos}(\tau _2)`$, $`A_1=\mathrm{cos}(\tau _1)\mathrm{cos}(\sqrt{2}\tau _2)`$ , and $`A_2=\mathrm{cos}(\sqrt{2}\tau _1)\mathrm{cos}(\sqrt{3}\tau _2)`$.

By using two atoms, each with a different interaction time, we have more freedom in locating solutions for which the magnitude of the $`A_i`$’s are closer together. This occurs at the expense of a more complex experimental scheme. To find solutions we employed a simulated annealing algorithm on an initial ensemble of randomly chosen points. After the points where suitably ‘cooled’ relative to a penalty function, we applied a simplex minimisation method on the top contenders to find the local minima. Some interaction times that result in implementing nearly ideal nonlinear-sign gates with high probability are given in table 2.

In an experiment it would be necessary to find a way to calibrate the interaction time until the desired phase shift had been reached. One way to do this is depicted in figure 2. Two single modes, each prepared in a one photon state are incident on a 50/50 beam splitter. The two photon interference then results in the state $`|2|0|0|2`$. The conditional phase shift can then be inserted on one arm so that the state is transformed to $`|2|0e^{i\theta }|0|2`$. This state is then run through an identical 50/50 beam splitter to the first and the probability for coincidence counts is sampled. This probability is given by $`P_c=\mathrm{cos}^2\theta /2`$. The coincidence detection rate then drops to zero at the required conditional phase shift of $`\theta =\pi `$.

We now estimate some typical values for the parameters. In a recent experiment a similar stimulated Raman process was observed using single rubidium atoms falling through a high finesse optical cavity Henrich2000 . The following parameters are typical of that experiment: $`g=2\pi \times 4.5`$ MHz, $`\mathrm{\Omega }=2\pi \times 30`$ MHz and $`\mathrm{\Delta }=2\pi \times 6`$ MHz. This gives a coupling constant of the order of $`70`$ Mhz. To achieve effective interaction constants of the order of those in the table 1 requires interaction times of the order of $`0.15`$ $`\mu `$s. In this paper we have neglected cavity decay which obviously needs to be kept small over similar time scale, which while difficult is not impossible.

We would like to thank the Computer Science department of the University of Waikato for making available their computing resources. AG was supported by the New Zealand Foundation for Research, Science and Technology under grant UQSL0001. GJM was supported by the Cambridge-MIT Institute while a visitor at University of Cambridge.