text:
# New form of matter at CERN SPS: Quark Matter but not Quark Gluon Plasma

## 1 Introduction

2000 has been a very exciting time in high energy heavy ion physics, marked by two press announcements: On 10 February, CERN summarized the results of its heavy ion program with claims related to the formation of a new state of matter in fixed target 158 AGeV Pb+Pb collisions . On 13 July, the BNL reported on the observation of the first collisions of Au + Au nuclei at the recently completed RHIC accelerator as detected by the STAR, the PHENIX, the PHOBOS and BRAHMS collaborations .

The purpose of the heavy ion program at CERN SPS and at BNL RHIC accelerators is to produce a new state of matter, the quark gluon plasma (QGP), where color degrees of freedom are deconfined and the basic degrees of freedom are quarks and gluons. Due to their relatively large color degrees of freedom, the properties of QGP are dominated by that of massless gluons. In particular, the equation of state is soft, the entropy density is high due to the large number of deconfined color degrees of freedom and due to this reason it takes a long time for the system to rehadronize. Direct photons are emitted from the QGP by the fractionally charged quarks and the production of the $`J/\psi `$ mesons is expected to be suppressed due to color screening.

Based on a critical review of the experimental results of the CERN heavy ion program (summarized in ref. ), I argue that such a QGP state has not yet been reached in the CERN SPS heavy ion experiments. The circumstantial evidence seems to point towards the formation of a quark matter (QM), where the dominant degrees of freedom are massive (dressed, constituent) quarks, but the of effective role of gluons is secondary. In the followings I will argue, that in contrast to the soft QGP equation of state, the equation of state of QM has to be hard to explain the strong three-dimensional expansion observed in the hadronic final state, based on a combination of the data analysis in refs. . The entropy density of QM seems to be relatively low ( with no experimental evidence for a dominant role of gluons). This low initial entropy density and the strong three-dimensional expansion of QM may then lead to a sudden rehadronization, where the abundances of directly produced particles are determined with the help of quark combinatorics as described in refs. . Such combinatorics may also govern $`J/\psi `$ and other charmed particle production, a topic that deserves further investigations .

### 1.1 Summary of CERN Announcement

On February 10, the official CERN press release  summarized the results of its heavy ion program as follows: “Compelling evidence now exists for the formation of a new state of matter at energy densities at about 20 times larger than that in the center of atomic nuclei and temperatures about 100000 times higher than in the center of the sun. This state exhibits characteristic properties which cannot be understood with conventional hadronic dynamics (i), but are qualitatively consistent with expectations from the formation of a new state of matter (ii) in which quarks and gluons no longer feel the constraints of color confinement (iii).” Both the CERN press release  and the summary manuscript of Heinz and Jacob  clearly and consciously distinguished between the claim for evidence for a new state of matter (that they claimed) and between the “discovery of QGP” (that they did not claim). The physical picture was summarized in the following straightforward manner:

“(1.) Two colliding nuclei deposit energy in the reaction zone. The energy materializes in the form of quarks and gluons, which strongly interact with each other.

(2.) This early, very dense state has an energy density of 3-4 GeV/fm3 and the equivalent of a temperature of around 240 MeV. The conditions suppress the number of J-$`\psi `$-s (charmonia), enhance strangeness, and begin to drive the expansion of the fireball.

(3.) The “quark-gluon plasma” cools down and becomes more dilute.

(4.) At an energy density of 1 GeV/fm3 (and a temperature of 170-180 MeV), the quarks and gluons condense into hadrons, and the final abundances of the different types of particles are fixed.

(5.) At an energy density of around 50 MeV/fm3 (and a temperature of 100-120 MeV) the hadrons stop interacting completely and the fireball freezes out. At this point it is expanding at over half the speed of light.”

It had been emphasized, that the evidence for the above QGP picture is circumstantial and this picture had been put together from many little observations just like a complicated jigsaw puzzle.

## 2 Controversies in the QGP picture

I think it is important to highlight some of the controversial points in the above QGP picture, that are more or less well known but may have not been summarized before.

It is well known that the search for a new state of matter in the CERN SPS heavy ion program focused on two kind of QGP signatures: the early, penetrating probes of QGP formation (photons and lepton pairs that do not participate in the hadronic processes after their creation) and and the late, hadronic probes that are produced when the relevant kinds of strong interactions become negligible, due to the expansion and the related rarification and cooling.

It is also well known that almost any calculation of e.g. direct photon production or dilepton emission is sensitive to the time evolution and the equation of state of the system, as the penetrating probes are emitted from the whole volume integrated over the time evolution of the hot and dense strongly interacting matter. Hence the emission pattern of the penetrating probes depends sensitively on the time evolution of the temperature and density profiles, for example.

The time evolution of the temperature or the local rest energy / entropy/ baryon densities, on the other hand, depend drastically on the dimension of the expansion. The time evolution has to satisfy the boundary condition, that it ends up on the proper hadronic final state that has been determined from the analysis of the late, hadronic signals. The hadronic final state thus imposes a severe constraint on the possible time evolution scenarios.

### 2.1 The hadronic final state

Let us follow the strategy of backwards extrapolation: It has been discovered in 1994-95 that the hadronic final state can be reconstructed only from a combined analysis of single particle spectra and two-particle Bose-Einstein correlations of pions and kaons . The hadronic final state of Pb + Pb collisions has been reconstructed with this method using the Buda-Lund hydro model in ref. . For the very different kind of experiments: NA44, NA49 and the preliminary WA98 data on central Pb+Pb collisions at CERN SPS, the average value of the temperature as well as the value for of the mean transverse flow, the transverse radius for the central and the surface temperature of the fireball similar fitted values were obtained, as well as for the mean proper-time of particle freeze-out and the width of the freeze-out time distribution . Furthermore, the mean transverse flow is found to be about the same as that of an independent analysis of ref. : the mean transverse flow at the transverse radius was found to be $`u_t=0.55\pm 0.06`$, the transverse radius parameter was $`R_G=7.1`$ fm, the mean freeze-out time after the onset of the scaling longitudinal expansion was found to be $`\tau _0=5.9\pm 0.6`$ fm . The central temperature at the mean freeze-out time was found to be $`T_0139\pm 6`$ MeV, the surface temperature after the particle production was about to end was $`T_s85`$ MeV, yielding an average freeze-out temperature of about 100 -120 MeV.

The presence of a quasi-linear, Hubble -like flow in the final stage of Pb+Pb collisions is very well established, based on the observed approximately linear rise of the effective slope parameters of heavier resonances with the increasing mass:

$$T_{eff}T_{}+mu_t^2$$
(1)
The reconstructed final state of ref.  has been extrapolated backwards in time by T. S. Biró  using exact quasi-analytic solutions of three-dimensional (3d) relativistic hydrodynamics assuming the presence of mixed Quark-Gluon Plasma - hadron gas phases, as given in Figure 3. The result implies many things:

i) Due to the strong 3d expansion, the system is able to convert a large amount of latent heat into hadrons in a relatively short time. Assuming a phase mixture of QGP and hadron gas, the QGP phase should have started to hadronize as early as $`\tau =1.14`$ fm/c. This time scale is very close to the canonical guess of $`\tau _i=1`$ fm/c that is frequently used in an over-simplified Bjorken formula  to estimate the initial energy density . Bjorken’s formula is based on the assumption of a 1 dimensional self-similar flow in the beam direction:

$$ϵ_{Bj}=\frac{1}{2\pi \tau _iR_{\text{rms}}^2}\frac{dE_t}{dy}$$
(2)
which implies

$$ϵ_{Bj}=2.8\pm 0.26\text{GeV/fm}\text{3}\mathrm{at}\tau _i=1.14\mathrm{fm}/\mathrm{c},$$
(3)
in contradiction with the result of the 3d relativistic hydrodynamical solution, where

$$ϵ_{3\mathrm{d}}=ϵ_c0.6\text{GeV/fm}\text{3}\mathrm{at}\tau _i=1.14\mathrm{fm}/\mathrm{c}$$
(4)
Thus the estimation of the initial energy density with the help of Bjorken’s formula is uncertain and can become unreliable quantitatively (a factor of 4 over-estimate) and qualitatively (because it assumes a 1d instead of a 3d expansion).

ii) Even if we assume that the initial state is an equilibrated, thermalized QGP consisting of (massless) gluons and quarks, Biró’s backward extrapolation implies that the volume fraction $`x`$ of QGP decreases below $`x<1/3`$ within the first 1.3 fm/c of the 3d expansion, and the rehadronization is fully completed within a time period of 6 fm/c. This in turn suppresses the production of all the penetrating probes and implies that the signal in direct photon production, $`J/\psi `$ suppression and dilepton production must be very week, much weaker than signals calculated for an optimistic one dimensional Bjorken-type expansion with a long-lived rehadronization.

iii) Inspecting figure 3 one finds that in order to generate the strong transverse flow by the end of the expansion as required by the experimental data, the transverse flow has to be even stronger at the beginning of the rehadronization than at the end of this phase transition, if a soft equation of state is assumed. Indeed, if the pressure is constant during the rehadronization, the flow can only decrease, due to the radial expansion. However, the backward extrapolation depends not only on the boundary condition (the final state requested by the data) but also on the equation of state (EOS). So, one may assume that the hypothesis of being at the softest point is incorrect, $`c_s^2>0`$. If the equation of state is hard, then the radial flow can be generated during the transverse expansion, perhaps even a Bjorken-type initial condition can be connected to the observed final state. The sensitivity of the extrapolation on the choice of EOS to a given final state has been studied recently by Schlei and collaborators , who performed a similar but forward extrapolation using a lattice QCD inspired QGP equation of state and a hadron gas equation of state. However, the Hylander model calculations under-estimated the transverse flow in the final state even when using a hard, hadron gas equation of state, which indicates that some of the observed final transverse flow may already be generated on the primordial level in the binary nucleon-nucleon collisions .

iv). Biró’s full 3d hydrodynamical backwards solution also that the time of hand-waving arguments and order of magnitude estimates is over for the expansion stage of the “Little Bang” fireballs.

This observations can be summarized as follows: a soft equation of state together with a zero or small primodial transverse flow is in disagreement with the boundary condition imposed by the hadronic final state. If one accepts the premiss that the strong transverse flow has not been present from the very beginning in the initial state, then the strong transverse flow of the observed hadronic final state must have been generated by a hard equation of state.

### 2.2 A few caveats on single particle spectra and two-particle correlations

Figure 4 indicates the $`M_t`$ scaling of the effective source sizes of various particles in central Pb+Pb collisions at CERN SPS energies. Such a scaling has been predicted in an analytic calculation of the Buda-Lund hydro model (BL-H) in refs. , with a large number of caveats that cannot be discussed here. See also ref.  for some of the details and most recent results. However, I would like to highlight one aspect of the BL-H calculation: in order to be able to calculate the scaling function precisely, and in order to obtain an approximate scaling in the transverse mass variable $`M_t`$, the BL-H model had to assume a non-vanishing value for the temperature inhomogeneity inside the particle emitting source. In the center of the plane transverse to the beam axis, the temperature of the source had to be slightly hotter than at the transverse r.m.s. radius at the mean freeze-out time, so that to keep the point of maximum emittivity close to the beam axis even for large values of $`M_t`$. This in turn implied a saturation of the effective slope parameters $`T_{}`$ with the increasing values of the particle mass,

$`T_{}`$ $`=`$ $`T_0+mu_t^2{\displaystyle \frac{T_0}{T_0+m\frac{\mathrm{\Delta }T}{T}_r}}`$ (5)
$`R_{}^2`$ $`=`$ $`{\displaystyle \frac{R_G^2}{1+\frac{M_t}{T_0}\left(u_t^2+\frac{\mathrm{\Delta }T}{T}_r\right)}}`$ (6)
and the competition between the transverse flow $`u_t`$ and the transverse temperature inhomogeneity $`\frac{\mathrm{\Delta }T}{T}_r`$ controls the $`M_t`$ dependence of the effective source sizes as well as the flattening of the initial linear mass dependence of the effective slope parameters of the single-particle spectra. Very heavy particles resolve the temperature inhomogeneities of the source, as the Boltzmann factor that characterizes their abundances focuses strongly their production to the hottest central parts, where the flow effects become limited.

The central temperature is  $`T_0140`$ MeV, the flattening of the slopes sets in at about $`m=1400`$ MeV, which implies about $`10\%15\%`$ temperature inhomogeneity in the Pb+Pb source . This amount is rather small and it is in agreement with the results from a combined analysis of the single-particle spectra and the two-particle correlations in ref. , which resulted in $`u_t0.55\pm 0.06`$ and $`\frac{\mathrm{\Delta }T}{T}_r0.06\pm 0.05`$, which suggests a slope parameter of about $`T_{eff}315`$ MeV for a particle with a mass of 1400 MeV.

For more details on this specific point and on a discussion of the significance of the temperature inhomogeneity and the transverse flow in shaping the transverse density profile, I recommend ref. .

This 10% transverse inhomogeneity in the central temperature becomes rather important when discussing the hadron-chemical composition. In my opinion, this temperature inhomogeneity prohibits the association of a single temperature value to the hadron-chemical composition. However, it allows for the interpretation of the observed $`M_t`$ scaling of HBT radii and the linearly rising than flattening effective slope parameters of the single particle spectra in a self-consistent and controllable manner, within the same framework. It is also clear, that the net baryon density is inhomogeneously distributed even in the central reaction zone due to surface effects, $`\mu _B=\mu _B(r_x,r_y)`$.

Due to the importance of the inhomogeneities of the temperature, and the net baryon number distributions (baryon chemical potential), I think that it is premature to discuss the rehadronization in Pb+Pb collisions at CERN SPS in the framework of a complete hadrochemical and thermal equilibrium. The simultaneous analysis of the two-particle correlation functions and the single-particle spectra indicates  that the production of the observable hadrons happens in a relatively narrow longitudinal proper-time interval, characterized by a width of $`\mathrm{\Delta }\tau =1.6\pm 1.5`$ to be compared with the mean freeze-out proper-time of $`\tau _0=5.9\pm 0.6`$ ( as measured from the onset of the scaling longitudinal expansion). By the time the particle production is over, the surface of Pb + Pb collisions cools down from 139 MeV to about 83 MeV . It is very interesting to note, that this value is similar to the surface temperature of $`T_s=82\pm 7`$ MeV found in h+p reactions as a consequence of transverse temperature inhomogeneities, see ref. . Such “snow-balls” with relatively low values of surface temperature and a possible hotter core were reported already in $`S+Pb`$ reactions in ref. .

Other hydro parameterizations, as reviewed in ref. , frequently neglect the effects of temperature inhomogeneities during the expansion and particle production stage. Energy conservation implies that the temperature cannot be exactly constant when particles are freezing out in a non-vanishing period of time from a three-dimensionally expanding source. Fixing the temperature to a constant in this time period, one finds some approximate average values of freeze-out temperatures in the range of $`T_f=110\pm 30`$ MeV.

In the physical situation of Little Bangs, expansion competes with the drop of the pressure gradients, which in turn is induced by the drop of the temperature on the surface. If the flow is small enough, a sudden drop of the temperature on the surface may result in a sudden drop of the pressure gradients on the surface, which implies density pile-up and a formation of a “ring of fire”, frequently seen in images of planetary nebulae as well, see ref.  for further details. On the other hand, if the flow is strong enough, it blows away the material from the surface, preventing the formation of such shells of fire, and an ordinary expanding fireball is obtained. The former case seems to be realized in $`h+p`$ reactions measured by the NA22 CERN experiment: a ring of fire is formed in the transverse plane due to the low transverse flow and due to the large temperature inhomogeneities. Pb+Pb collisions at CERN SPS are characterized with large transverse flows and relatively small transverse temperature inhomogeneities, hence they correspond to a exploding fireball with a more uniform, close to Gaussian density distribution . This result indicates that non-trivial time-evolution of fireball hydrodynamics is intimately connected with the spatial inhomogeneities of the temperature and the corresponding density profile.

## 3 The new-old picture: Quark Matter

If one accepts the arguments of the previous section about the important role of temperature inhomogeneities in the source, the question arises: can we draw a phase diagram about the rehadronization process? Can we introduce the concept of a unique hadrochemical and kinetic temperature? In my opinion, these concepts are limited by the above shown variation of the local temperature during particle production (e.g. between 140 to 85 MeV at the kinetic freeze-out) so perhaps a ball-park value can be given, but the non-homogeneity is important and even in idealized cases the precision of a kinetic or chemical freeze-out “temperature” cannot be decreased below the 20 - 30 % relative error level. Thus the “data points” on the beautifully drawn phase diagrams e.g. in ref.  have large systematic uncertainties, and due to this reason I think it is premature to conclude at present about the separation of a hadrochemical freeze-out temperature of the order of 175 MeV from the and kinetic freeze-out temperature of the order of 110 MeV as well as about the precise value of the baryo-chemical potential and the temperature in various reactions: within 20 -30 % relative systematic error, these values can be the same. The validity of the method to extract these points can be questioned not only because it relies on the concept of spatially homogeneous temperature distributions, but also because the method yields a well-defined value for the hadrochemical and kinetic freeze-out temperatures even for $`e^+e^{}`$ reactions at LEP . However, we know that the $`e^+e^{}`$ reactions are characterized by jet production and non-thermal fluctuation patterns like jets within jets within jets etc. Hence the well defined values of the baryon chemical potentials and the freeze-out temperatures in these kind of analyzes seem to be more characteristic to the method than to the physical system under consideration.

If the concept of the chemical freeze-out is maintained and embedded in a model of an exploding fireball, a statistically acceptable $`\chi ^2`$ fit to the observed hadronic abundances by Rafaelski and collaborators resulted  in a chemical freeze-out temperature of the order of 140 MeV, within errors similar to the central kinetic  freeze-out temperature obtained from the statistically acceptable $`\chi ^2`$ fits to the observed single-particle spectra and two-particle correlation data of NA44, NA49 and WA98 collaborations using the Buda-Lund hydro parameterization . Hence a clear separation of the chemical and thermal freeze-out, (steps 4) and 5) of the Introduction) cannot be taken for granted at present.

## 4 Big Bang and Little Bang

Table 1 briefly summarizes the similarities between the physics of the Big Bang that resulted in our Universe and the physics of high energy heavy ion collisions or “Little Bangs” that are studied in the laboratory at CERN SPS and at Brookhaven AGS and RHIC accelerators.

## 5 Strangeness enhancement

The enhancement of strange particle production has been long thought to carry signals of QGP formation. Strangeness enhancement was intimately related to the question of chemical equilibration times in a QGP.

In the earliest discussions of high energy heavy ion collisions it was assumed, that in this reactions a long lived, thermally and chemically equilibrated QGP is formed. In such a QGP a large amount of $`s\overline{s}`$ pairs can be formed, hence the enhancement of strangeness production was proposed as a signature of QGP formation .

Later, however, it was questioned if the equilibrium value of the $`s\overline{s}`$ number can be reached or not during the lifetime of a QGP phase, and it was shown, that the rate of the $`q+\overline{q}s\overline{s}`$ reaction is too small to reach the equilibrium values  Following this observation, Rafelski and Müller showed that the inclusion of the $`g+gs\overline{s}`$ reaction a large enough strangeness production rate is achieved  which seems to be enough to reach the equilibrium value of the concentration of $`s\overline{s}`$ pairs.

The problem of the recombination of quarks and antiquarks into hadrons was studied by Biró and Zimányi in ref. . In this work the mechanism of hadronization was assumed to be a non-linear, ( quark number conserving) coalescence process. A few years later this model was extended by taking into account the effects of gluons: the gluons were assumed to fission into $`q\overline{q}`$ and $`s\overline{s}`$ pairs, increasing the number of quarks and antiquarks entering into the hadronization process described above. Further, these authors assumed a hadrochemical evolution after the hadronization in order to obtain the final hadron numbers .

If a quark-gluon plasma state is formed in a heavy ion reaction, one thus expects that the $`s`$-quark distributions are equilibrated in a few fm/c, and that the production of $`s\overline{s}`$ quarks is enhanced as compared to strange particle production in normal hadronic interactions, because the following reasons: i) In a QGP, the dominant degrees of freedom are the gluons that can easily enhance the strangeness content in the gluon fusion process $`g+gs\overline{s}`$ ; ii) at $`TT_c`$ chiral symmetry is (at least partially) restored, and the mass of strange quarks is expected to decrease to $`m_s150`$ MeV, which implies that the thermal production rate is relatively high, $`N(s)e^{m_s/T}1`$, iii) in a baryon-rich Quark Gluon Plasma, the Pauli blocking of the $`u`$ and $`d`$ quarks favours (at SPS energies also) the $`s\overline{s}`$ production over the $`u\overline{u}`$ or $`d\overline{d}`$ production. Finally, one expects that the strange quarks are converted into hyperons during the rehadronization and their abundances are thus enhanced reflecting the reduced threshold of strange quark production in the deconfined phase.

The question arises: how important is the dominance of gluons in the above picture? Actually, only the following feature matters: that the production threshold for a $`s\overline{s}`$ creation is reduced and that this process starts to compete with the production of light quark pairs. Zimányi and Biró studied the problem of how quarks recombine  into hadrons after the gluons have already disappeared from the system. For clarity, we shall refer to this situation as to the Quark Matter (QM), which name does not include the name of gluons, in contrast to the QGP acronym. The kaon and hyperon enhancement has been observed to be a common feature of both QM and QGP . The question arises, what are the most natural observables and can one generalize the results to non-equilibrium situations?

I have argued that particle production in CERN SPS heavy ion reactions happens most likely in a sudden, non-equilibrium manner. What are the consequences of this mechanism, can one quantify e.g. the effect of such a mechanism on the hadrochemical composition of the produced particles?

Fortunately, the answer to this theoretical question is positive, provided that hadron production is characterized by the sudden recombination of constituent quarks into hadrons, as described by the ALCOR model . Also, the experimental results at CERN SPS Pb+Pb measurements clearly and convincingly indicated a huge enhancement in the production of strange particles .

Bialas realized recently, that simple relations hold between various (multi-strange) antibaryon to baryon ratios  and the Budapest group has proven that these relations hold not only in a linear approximation to quark recombination, but also they are valid in general, if the hadron production happens through a sudden and complete, non-linear recombination of constituent quarks into hadrons . The simplest formulation of such a recombinative hadron production method is described by the ALCOR model (ALgebraic COalescence for Rehadronization) . The particle abundances are connected by the following simple relations :

$`{\displaystyle \frac{\overline{\mathrm{\Lambda }}+\overline{\mathrm{\Sigma }}}{\mathrm{\Lambda }+\mathrm{\Sigma }}}`$ $`=`$ $`{\displaystyle \frac{\overline{N}}{N}}\left[{\displaystyle \frac{K}{\overline{K}}}\right],`$ (7)
$`{\displaystyle \frac{\overline{\mathrm{\Xi }}}{\mathrm{\Xi }}}`$ $`=`$ $`{\displaystyle \frac{\overline{N}}{N}}\left[{\displaystyle \frac{K}{\overline{K}}}\right]^2,`$ (8)
$`{\displaystyle \frac{\overline{\mathrm{\Omega }}}{\mathrm{\Omega }}}`$ $`=`$ $`{\displaystyle \frac{\overline{N}}{N}}\left[{\displaystyle \frac{K}{\overline{K}}}\right]^3.`$ (9)
In these equations, $`N`$ is the number of directly produced nucleons, $`\overline{N}`$ is the directly produced anti-nucleons, $`\mathrm{\Lambda }`$ is the number of directly produced $`\mathrm{\Lambda }`$ baryons etc, so care must be taken when comparing with the experimental data, because of corrections from the resonance decays. Recent data from the WA97, NA44 and NA49 experiments indicate, that these relations are satisfied (after resonance decays are corrected for) in central Pb+Pb collisions at CERN SPS. This can be considered as a model-independent proof  that constituent quark degrees of freedom are liberated in these reactions and hadron production proceeds via a sudden and complete coalescence of constituent quarks to hadrons in central Pb+Pb collisions at CERN SPS.

At this point, it is natural to ask whether such constituent quark degrees of freedom reveal themselves in p+p or p+A collisions or not? This question has been investigated throughoutly by the E910 experiment at the AGS . The detailed analysis of the dependence of leading baryon production, strange particle production and pion production suggests that i) baryon stopping proceeds through a mechanism that is different from energy stopping; ii) most of the energy carried away from precursors of energetic pions can be used for strange particle (e.g. $`\mathrm{\Lambda }`$ production), iii) a break-up picture of the projectile to constituent quarks is a possible explanation of the $`\mathrm{\Lambda }`$ and $`K_S^0`$ production, which is very close in spirit to the additive quark model.

In a QGP picture, strangeness enhancement happens through the large gluon density that create $`s\overline{s}`$ pairs by gluon fusion. Chiral symmetry restoration reduces the constituent mass of strange quarks from 450 MeV to about 150 MeV, thus reducing their production threshold .

The question arises, do the gluonic degrees of freedom indeed have to play a dominant role for an enhancement of strange quark production? I think that the key requirement for strangeness enhancement is the reduction of the production threshold for $`s\overline{s}`$ pairs. If the quarks are confined, the even the minimal excitations require the presence of additional light quarks to make a hadron; if a Quark Matter is formed, constituent $`S\overline{S}`$ pairs can be created without the need for additional light quarks that also reduces the production threshold. Hence strangeness enhancement can be expected both in case of QM and in case of QGP formation, if the $`S\overline{S}`$ production treshold is reduced.

A very interesting theoretical explanation of the process of enhanced strangeness production and a hard equation of state was given by Lévai and Heinz in ref. . Using effective, massive quarks and gluons to describe the lattice QCD equation of state, they observed that the effective mass of gluons is strongly increased if one approaches from above the critical temperature. At the same time, the mass of quarks approached the constituent mass and the speed of sound remained rather large, $`c_s^2>0.15`$. The increase of the effective gluon mass with decreasing temperature may explain why the gluonic degrees of freedom seem to be less evidently required by the Pb+Pb data, than that of the constituent quarks .

The picture of the dominant role of constituent quarks and constituent antiquarks in the strangeness production, the integration of the gluonic degrees of freedom into an effective, confining equation of state was emphasized recently by the calculations of Biró, Lévai and Zimányi in refs. . The necessity of a fast rehadronization and a three-dimensional expansion has been realized already in ref. .

## 6 Charm production and $`J/\psi `$ suppression

In a QGP, color degrees of freedom are liberated, both quarks and gluons become active degrees of freedom. If a $`c\overline{c}`$ pair is created from a gluonic fusion, the color interaction between these quarks is screened and the formation of charmonium states, in particular the creation of $`J/\psi `$ mesons, is suppressed, as predicted by Matsui and Satz .

There are a number of controversies related to the presentation of the NA50 data on $`J/\psi `$ suppression. Let me mention some of them: The horizontal axes of some of the figures indicates $`L`$ a variable that is thought to be characterizing the length of a path that the $`J/\psi `$ has to travel inside a medium. In a three-dimensionally expending, rarifying and cooling fireball, such a variable is difficult to define, not mentioning the problem how to determine it experimentally. On another plot, the ratio measured/expected is plotted versus the initial energy density as calculated from Bjorken’s formula. Not only the theoretical expectations differ from model to model, but also the uncertainty in the initial energy density is at least as big as a factor of 4. Experimental data points should be determined as a function of measurable quantities, and not as a function of theoretical calculations.

A particularly interesting plot has been predicted theoretically by Kharzeev, Nardi and Satz : the well measurable mean transverse momentum $`p_t^2`$ of the $`J/\psi `$-s in a hadron gas increases monotonically with increasing transverse energy, due to increased number of rescatterings with increasing centrality. However, if a QGP phase is reached in the center of the collision zone, the production of $`J/\psi `$-s are suppressed there, which implies that the $`p_t`$ of the $`J/\psi `$-s starts to decrease with increasing transverse energy $`E_t`$ produced in the experiment, after an initial rise. In some sense, this decrease of the $`p_t^2`$ of the $`J/\psi `$-s is due to the softening of the equation of state when a QGP is produced.

This plot has been compiled from the available NA50 data by J. Nagle in ref. . The result did not show the expected decrease of the $`p_t^2`$ with increasing $`E_t`$, but followed a pattern similar to the prediction for a hadron gas equation of state , see figure 6. This result casts doubt on a $`J/\psi `$ suppression from a QGP in central Pb+Pb collisions at CERN SPS. However, the detailed study performed in ref.  indicates that the $`E_t`$ dependence of the $`J/\psi `$ yield is incompatible with a Glauber-type calculation without initial energy loss, the total $`J/\psi `$ yield can be reproduced only with the help of an approximately 15 % initial energy loss, see fig. 7. Using this energy loss, the mean $`p_t`$ transverse momenta of the $`J/\psi `$-s is under-predicted.

However, this result may be a fingerprint of the action of a hard equation of state, and a strong 3-dimensional expansion in the hot hadronic matter, that eventually enhances the $`p_t`$ of the $`J/\psi `$-s (a transverse flow effect) in agreement with the results of the combined analysis of the single particle spectra and the two-particle correlation functions of pions, kaons and protons in the hadronic final state.

If one assumes that quark degrees of freedom are liberated but gluonic degrees not, and that the resulting quark matter has a hard equation of state, the abundances of various charmed mesons and baryons can be calculated with the help of the extension of quark combinatorics to the charm flavor. The resulting ALCORc model implies simple relationships between the ratios of (multi)charmed antibaryon to baryon ratios, generalizing the results presented in the section on strangeness .

It is particulary interesting to note, that the following simple relationships are predicted by ALCORc for the multi-charmed baryon/anti-baryon ratios:

$`{\displaystyle \frac{\overline{Y}_c}{Y_c}}`$ $`=`$ $`{\displaystyle \frac{\overline{N}}{N}}\left[{\displaystyle \frac{D}{\overline{D}}}\right],`$ (10)
$`{\displaystyle \frac{\overline{\mathrm{\Xi }}_{cc}}{\mathrm{\Xi }_{cc}}}`$ $`=`$ $`{\displaystyle \frac{\overline{N}}{N}}\left[{\displaystyle \frac{D}{\overline{D}}}\right]^2,`$ (11)
$`{\displaystyle \frac{\overline{\mathrm{\Omega }}_{ccc}}{\mathrm{\Omega }_{ccc}}}`$ $`=`$ $`{\displaystyle \frac{\overline{D}}{D}}\left[{\displaystyle \frac{D}{\overline{D}}}\right]^3,`$ (12)
$`{\displaystyle \frac{\overline{\mathrm{\Omega }}_c}{\mathrm{\Omega }_c}}`$ $`=`$ $`{\displaystyle \frac{\overline{\mathrm{\Omega }}}{\mathrm{\Omega }}}\left[{\displaystyle \frac{\overline{D}_s}{D_s}}\right],`$ (13)
$`{\displaystyle \frac{\overline{\mathrm{\Omega }}_{cc}}{\mathrm{\Omega }_{cc}}}`$ $`=`$ $`{\displaystyle \frac{\overline{\mathrm{\Omega }}}{\mathrm{\Omega }}}\left[{\displaystyle \frac{\overline{D}_s}{D_s}}\right]^2,`$ (14)
$`{\displaystyle \frac{\overline{\mathrm{\Omega }}_{ccc}}{\mathrm{\Omega }_{ccc}}}`$ $`=`$ $`{\displaystyle \frac{\overline{\mathrm{\Omega }}}{\mathrm{\Omega }}}\left[{\displaystyle \frac{\overline{D}_s}{D_s}}\right]^3,`$ (15)
and the mesonic step factors are related by a simple relationship,

$$\frac{\overline{D}_s}{D_s}\frac{\overline{D}}{D}=\frac{\overline{K}}{K}.$$
(16)
Quark combinatorics predicts not only the rates but also the slope parameters of charmed mesons as well, as calculated recently in refs. .

The measurements of multi-charmed anti-baryon to baryon ratios , the strange/charmed mesonic ratios of eq. (16), as well as the effective $`m_t`$ slopes $`D`$ and $`J/\psi `$ mesons can thus provide an important constraint and test of the hadronization process, and can exclude or confirm the possibility that charmed hadron production happens through quark recombination and coalescence similarly to that of the strange hadrons.

## 7 Penetrating probes

Due to lack of space, time, and expertise, I cannot discuss these results in detail. As I mentioned in the section on the particle spectra and correlations, at CERN SPS the final state indicates a strong 3-dimensional expansion with a hard equation of state, which is very efficient in reducing the large energy and entropy densities, which implies the reduction of the signal (if any) carried by the penetrating probes.

### 7.1 Direct photons

The experimental situation on direct photon production has been summarized recently in ref.  and can be briefly recapitulated as follows:

The relative increase of direct photons over photons from a conventional hadronic background is $`(N_\gamma N_{\gamma ,hadr})/N_{\gamma ,hadr}=12\%\pm 0.8\%\pm 10.9\%`$, a non-significant value. Within the errors this non-significant excess is constant as a function of multiplicity. A systematic study indicates, that the excess of directly produced photons in the most central $`Pb+Pb`$ collisions (as compared to the hadronic event generator VENUS) is less, than the same excess in Pb + Nb or in Pb+Ni collisions , see Figure 8. Indirectly, this result suggests the lack of significant amount of deconfined, undressed light quarks in the most central Pb+Pb collisions at CERN SPS.

Photon production is dominated by the $`\pi ^0\gamma +\gamma `$ decay and in the observed distributions there is no space for a larger than 10 % contribution of thermally produced photons from a Quark-Gluon Plasma.

This result of the CERES/NA45 experiment is consistent with the picture of a strong three-dimensional expansion which is able to reduce the QGP fraction (if any) to less than 1/3 within the first 1.3 fm/c of rehadronization .

One may expect that a long lived QGP is signaled with an enhanced direct photon production together with late freeze-out times and large widths of the freeze-out time distribution, as observed from two-particle correlation studies ($`R_{out}>>R_{side}`$). Thus, the lack of direct photon enhancement is consistent with a non-significant QGP production scenario and with the experimental results on a sudden particle freeze-out, $`R_{out}R_{side}`$.

### 7.2 Dilepton production

The CERES/NA45 experiment observed the inclusive $`e^+e^{}`$ invariant mass spectra and compared the result to the expectation based on dilepton production from a chemically equilibrated, thermalized hadron gas. As compared to this expectation, factor of $`2.6\pm 0.5\pm 0.6`$ enhancement of dilepton pairs in the $`0.25<m_{ee}<0.7`$ interval has been reported. The enhancement was shown to be concentrated on the low transverse momentum region, 200 MeV $`<p_t^{ee}<500`$ MeV, with almost negligible enhancement in the interval 500 MeV $`<p_t^{ee}`$. The various theoretical explanations of this effect focused on hadron modification in dense matter, as summarized recently in ref. .

However, the results can also be interpreted so that the hadron gas is not in full chemical equilibrium. An enhancement of the low $`p_t`$ ($`\eta `$,) $`\eta ^{}`$ and $`\omega `$ mesons by a factor of 5 seems to be able to describe the observed enhancements, compare the left and right panels of Fig. 6 in ref. . Such an enhancement has been proposed to signal the onset of a partial $`U_A(1)`$ symmetry restoration and should be detectable with the help of the measurement of the intercept parameter of the two-pion correlation functions as a function of the transverse mass, $`m_t`$ .

The important point of the above paragraphs is the following: the $`m_t`$ dependent production of $`\omega `$, $`\eta `$ and $`\eta ^{}`$ is constrained by the strength of the two-pion Bose-Einstein correlation function, and this constraint can be utilized to calibrate the expected number of hadronic contributions to dilepton decays in the range of the observed excess. Such consistency checks between two-pion correlation measurements and dilepton production data have not yet been performed as far as I know.

## 8 Hard Quark Matter – Soft QGP

Let me summarize the results of the CERN heavy ion program by slightly modifying the text of the original announcement (where the modifications are given in italics): “Compelling evidence now exists for the formation of a new state of matter at energy densities at about 5  times larger than that in the center of atomic nuclei… This state exhibits characteristic properties which cannot be understood with conventional hadronic dynamics (i), but are qualitatively consistent with expectations from the formation of a new state of matter (ii) in which valence quarks no longer feel the constraints of color confinement (iii) and the properties of matter are dominated by that of valence quarks, the role of gluons is secondary.”

The physical picture of Quark Matter formation can be summarized in the following straightforward manner:

(1.) Two colliding nuclei deposit energy in the reaction zone. The energy materializes predominantly in the form of constituent quarks, which strongly interact with each other.

(2.) This early, very dense state has an energy density of about 1 GeV/fm3 and the equivalent of a temperature of $`170180`$ MeV. The conditions enhance strange quark production and the hard equation of state begins to drive a strongly three-dimensional expansion of the fireball.

(3.) Due to the three dimensional expansion, the “quark matter” suddenly cools down and becomes very dilute.

(4.) The charmed, strange and light valence quarks recombine to form the hadrons. The allowed range for the “chemical” and the “kinetic” freeze-out temperatures include 140 MeV, so the chemical and the kinetic freeze-out happens almost simultaneously.

5.) At the time of the last interaction, the transverse expansion is characterized by about 10% temperature inhomogeneities , which imply large inhomogeneities (edge effects) in the transverse baryon density (and the baryon chemical potential). At the transverse rms radius the produced hadronic matterc expands at over half the speed of light.”

As mentioned, the evidence for the above Quark Matter picture is circumstantial and this picture had been put together from many little observations just like a complicated jigsaw puzzle. This picture is significantly different from the also circumstantial Quark Gluon Plasma picture reviewed in the introduction.

## 9 Summary : Quark Matter $``$ Quark Gluon Plasma

A new picture is presented here, of a formation of a Quark Matter at CERN SPS. Quark Matter is a new state of matter that can be created in high energy heavy ion collisions. Starting from the reconstructed final state of Pb + Pb collisions, certain inconsistencies are pointed out in earlier attempts that tried to put together the jigsaw picture of these heavy ion collisions on the basis of the formation of a Quark Gluon Plasma.

Constituent (valence) Quark Matter is a new state of matter, that is different from both ordinary hadronic matter and from the much expected Quark Gluon Plasma.

I think that we have to keep an open eye and look for new, unexpected phenomena. Our community has to be able to distinguish between the breakup of hadrons to a Quark Matter consisting of massive constituent quarks and between the complete melting of QCD matter to a Quark Gluon Plasma, that consists of almost massless quarks and gluons.

Quark model provides a successful description of the bulk data on hadron spectroscopy. It is not surprising that the same, constituent quark degrees of freedom are liberated in high energy nuclear collisions. Note that valence gluons are not required in hadron spectroscopy, but valence (constituent) quarks are required, as they carry conserved quantities like charge and baryon number. It is difficult to remove a constituent quark from a Quark Matter, processes like $`Q+QQ`$ are forbidden due to the conservation laws. Also, quarks are fermions so they cannot occupy the same quantum state due to Pauli blocking. I think these are the essential reasons why a Quark Matter has to be characterized by a hard equation of state.

On the other hand, if a Quark Gluon Plasma is produced, the number of degrees of freedom increases drastically and the channel $`g+gg`$ opens. As the abundant gluons do not carry any conserved charge, their number can be changed relatively easily. Also, gluons are bosons so any number of them can occupy the same quantum state. I think these are the essential reasons, why a QGP has to be characterized by a soft equation of state. For a more detailed analysis of the QCD equation of state in terms of massive quarks and gluons, let me recommed the work of Lévai and Heinz , whose phenomenological predictions (hard equation of state near the critical temperature etc) are in agreement with the qualitative analysis of the Pb+Pb data at CERN SPS.

It seems that constituent quark degrees of freedom play an important role in Pb +Pb collisions at CERN SPS energies, which is clearly demonstrated by the data on strange particle production. However, the liberation of gluonic degrees of freedom, and the softening of the equation of state due to the dominant role of gluons in a QGP is in disagreement with the available data at present, for example the mean transverse momentum of the $`J/\psi `$-s does not decrease with increasing transverse energy of the events.

A possible interpretation of this result is that a Quark Matter has been created in central Pb+Pb collisions at CERN SPS, however, a Quark Gluon Plasma phase has not yet been reached.

The discovery of QGP can be expected at higher initial energy densities, such as produced by the recently started RHIC accelerator. It is a very exciting challenge for RHIC to distinquish between various new forms of matter, for example between Quark Matter and Quark Gluon Plasma.

I would like to add, that the constituent quark dominated Quark Matter is described in earlier publications on strangeness production under various, perhaps more fancy names like CQP = constituent quark plasma, $`Q\overline{Q}P`$ = quark anti-quark plasma. This picture is somewhat similar to the valon model as well, although the valons were invented to consider soft scattering problems between quarks, as described in refs. , while the constituent quarks in the present picture emerge due to the combinatorical description of hadron production which is a bound-state formation problem.

In the light of the recent experimental data on single particle spectra, two-particle correlations, strange particle production, direct photon and dileption production as well as charmed particle production, as summarized in refs.  the picture of a formation of a Quark Matter (not QGP) emerges. This is in a qood agreement with the phenomenological analysis of the lattive QCD equation of state by Lévai and Heinz: Instead of an idealized, asymptotically free quark gluon plasma, the lattice results for the QCD equation of state can be parameterized phenomenologically as a micture of massive quarks anti-quarks, while the gluons pick up a large mass near $`T_c`$ so that they stop to play a dominant role, and as a consequence the QM has hard equation of state . The picture of Quark Matter (not Quark Gluon Plasma) formation at CERN SPS can be summarized in the following manner:

(1.) Two colliding nuclei deposit energy in the reaction zone. The bombarding energy breaks the nucleons into constituent quarks and materializes in the form of new constituent quark–anti-quark pairs. (2.) This early, very dense state has an energy density of about 1 GeV/fm3 and the equivalent of a temperature of more than 140 MeV. The conditions enhance strangeness, suppress the number of gluons, and the hard equation of state drives a strong transverse and longitudinal expansion of the fireball.

(3.) The Quark Matter cools down fast due to the strong three-dimensional expansion.

(4.) At the critical temperature of deconfinement, the constituent quarks suddenly recombine into hadrons, and the final abundances of the different types of directly produced particles are fixed. The number of charmed, strange and light mesons and baryons is determined from quark combinatorics and feed downs from the decay of hadronic resonances.

(5.) The chemical and thermal freeze-out temperature distributions may be rather close to each other, in the range of 140 MeV. The density distribution (baryon chemical potential) is rather inhomogeneous, and the central temperature at the mean kinetic freeze-out time decreases from about 140 MeV to about 85 MeV at the surface. At the time of the last interactions the fireball is expanding transversally at over half the speed of light.”

Table 2 summarizes briefly the main phenomenological and experimental differences between an idealized, soft Quark Gluon Plasma, and the old-new, non-ideal hard Quark Matter, whose production in central Pb + Pb reactions at CERN SPS seems to be consistent with the current, exciting experimental situation.

Figure 9 illustrates the phenomenological separation of the liberation of consituent (dressed) quark degrees of freedom in form of Quark Matter from the theoretically predicted liberation of light (undressed) quarks and gluons in very energetic heavy ion collisions. If the gluons carry an effective (temperature and density dependent) mass, the gluonic degrees of freedom appear at around the $`T=m_{}(T,\mu _B,\mathrm{})`$ line, which effectively separates Quark Matter from Quark Gluon Plasma on this illustrative example of fig. 9. Note however, that the “line” separating the QM from QGP may correspond to a crossover instead of a phase transition in a strict sense, similarly to that of the transition of a heated mono-atomic gas to a plasma state. At this line, the dominant degrees of freedom change, but the functional form of the equation of state may not necessarily change.

Finally, Table 3 summarizes the 5 simple but different steps of a hadronizing QGP and a hadronizing Quark Matter.

## Acknowledgments

I would like to thank B. Cole, R. Hwa, P. Lévai and J. Zimányi for inspiring discussions, the Organizers of Torino 2000 for creating a pleasant atmosphere and an inspiring working environment.

This research was supported by a Bolyai Fellowship of the Hungarian Academy of Sciences and by the grants OTKA T024094, T026435, T029158, the US-Hungarian Joint Fund MAKA grant 652/1998, NWO - OTKA N025186, OMFB - Ukraine S& T grant 45014 and FAPESP 98/2249-4 and 99/09113-3.


format: text

text:
# Exploding and Non-exploding Stars: Coupling Nuclear Reaction Networks to Multidimensional Hydrodynamics

## Introduction

Thermonuclear reactive flows are ubiquituous in astrophysics and occur in non-explosive environments as, e.g., in most (hydrostatic) stars as well as in explosive events, for which novae and supernovae are examples. Often they provide the energy which powers stellar outbreaks (as in the case of novae, X-ray flashes, and thermonuclear, i.e. Type Ia, supernovae) and even for stellar explosions where this is not the case (as e.g. in core collapse supernovae, which are driven by neutrino heating), the strong coupling of hydrodynamic advection and thermonuclear reactions is of utmost importance for the nucleosynthesis which accompanies these events. It is by a proper numerical modelling of this coupling through which a more detailed insight into the origin of the nuclear abundances in the solar system can be gained, which are themselves the result of a superposition of material which has been processed in explosive and non-explosive thermonuclear environments. By comparing the results of numerical models with the observed solar abundance pattern, on the other hand, one might also hope to learn more about the thermodynamic conditions in the otherwise unaccessible nucleosynthetic sites and events themselves.

The high precision with which nuclear abundances can be measured nowadays poses great demands on the accuracy of the numerical models, especially since it was convincingly demonstrated in recent years that due to the importance of hydrodynamic instabilities, rotation, and other effects, most of the nucleosynthetic sites do not possess spherical symmetry. Thus a reliable computation of the highly non-linear interaction of hydrodynamic advection and nuclear burning requires multidimensional numerical models. In the following sections we give a general overview of the methods which are currently employed for modelling thermonuclear flows and discuss some of the problems which are hereby encountered. Further reviews on reactive flow modelling can be found in LeVeque98 , Mueller98  and Oran\_Boris87 .

## The Governing Equations

A rather wide range of astrophysical reactive flows, in which relativistic effects, viscosity and magnetic fields can be neglected, is described by the well-known (reactive) Euler equations. This system of non-linear partial differential equations which expresses the conservation of the total mass, momentum, total (i.e. kinetic + internal) energy and baryons of the fluid reads

$`{\displaystyle \frac{\rho }{t}}+\left(\rho 𝐯\right)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{\rho 𝐯}{t}}+\left(\rho \mathrm{𝐯𝐯}\right)+P`$ $`=`$ $`\rho 𝐠+\rho 𝐟_{\mathrm{add}}`$ (2)
$`{\displaystyle \frac{\rho E}{t}}+\left(\left[\rho E+P\right]𝐯\right)`$ $`=`$ $`\rho 𝐯𝐠+\rho \dot{Q}_{\mathrm{add}}+\rho \dot{Q}_{\mathrm{nuc}}`$ (3)
$`{\displaystyle \frac{\rho X_i}{t}}+\left(\rho X_i𝐯\right)`$ $`=`$ $`\rho \dot{X}_i`$ (4)
$`{\displaystyle \underset{i}{}}X_i`$ $`=`$ $`1,`$ (5)
where $`\rho `$, $`𝐯`$, $`E=v^2/2+e`$, and $`P`$ have their usual meanings, $`X_i`$ is the mass fraction of nucleus $`i`$, and $`\rho \dot{X}_i`$ as well as $`\rho \dot{Q}_{\mathrm{nuc}}`$ are source terms due to nuclear transmutations. If self-gravity is important, the gravitational acceleration

$$𝐠=\mathrm{\Phi }$$
(6)
which appears in the source terms of Eqs. (2) and (3) and which depends on the gravitational potential, $`\mathrm{\Phi }`$, has to be obtained from a solution of the Poisson equation

$$\mathrm{\Delta }\mathrm{\Phi }=4\pi G\rho .$$
(7)
In mathematical terms Eqs. (17) describe a mixed initial/boundary value problem due to the hyperbolic and elliptic nature of the Euler and Poisson equations, respectively. Given appropriate initial and boundary conditions, an equation of state relating $`\rho `$, $`P`$ and $`e`$, and the additional source terms $`𝐟_{\mathrm{add}}`$ and $`\dot{Q}_{\mathrm{add}}`$, which in general will be problem-dependent, Eqs. (17) can be solved after an appropriate flow representation as well as a suitable numerical scheme have been adopted.

## Flow Representations and Numerical Schemes

There are two primary approaches to solve the homogeneous part of the system of equations (15). In the Eulerian framework the system of conservation laws is solved on a grid which is fixed in space and the evolution of the flow is followed by advecting the fluid through the computational cells. The principal assets of this method are its straightforward extension from one to two or three spatial dimensions and the simplicity of its implementation on serial and parallel computer architectures. If an appropriate shock-capturing, finite-volume numerical scheme is used, it is equally straightforward to obtain strict numerical conservation of all physically conserved quantities and a sharp resolution of shocks. The major drawback is numerical diffusion. Consider the continuity equation (1) in its Eulerian form, which can be written as

$$\frac{\rho }{t}+𝐯\rho +\rho 𝐯=0,$$
(8)
where the second term describes advection and the third term compression. Numerical diffusion is introduced into a numerical solution of this equation as a result of discretization errors of the $`𝐯`$ operator. There appears to be a simple remedy to this problem: using the comoving derivative $`\mathrm{d}/\mathrm{d}t=/t+𝐯`$ we can rewrite Eq. (8) in the frame comoving with the matter to obtain its Lagrangian form

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}+\rho 𝐯=0.$$
(9)
Note that in this frame the advection term $`𝐯\rho `$ has vanished. Therefore the Lagrangian approach is (in principle) not prone to numerical diffusion of mass (or composition). In Lagrangian methods each cell of the numerical grid represents a discretized fluid element which evolves subject to forces which are due to interactions with its neighbors and the time rate of change of the density of such a fluid element is solely determined by the compression (or expansion) that it experiences. Density interfaces (contact discontinuities) as well as composition discontinuities can be easily aligned with the boundaries of grid cells and do not have to be advected through the grid in the course of the calculation.

While this very desirable property of the Lagrangian approach has made it the method of choice for one-dimensional nucleosynthesis calculations, considerable difficulties are experienced when Lagrangian schemes are applied to multidimensional flows. Shear and vortices can severely distort a Lagrangian grid. The discrete approximation of differential operators over such a grid results in large errors in the numerical derivatives, and in the extreme case that the grid lines cross (grid tangling) the calculations have to be stopped. Some remapping procedure to a new, more regular grid must then be applied which unavoidably introduces numerical diffusion to the solution. The distortion problem can be overcome if triangular instead of quadrilateral grids are used Oran\_Boris87  or if (as in the Smoothed Particle Hydrodynamics, or SPH approach) no grid at all is adopted and instead the flow is sampled by a finite number of particles. In the former method considerable logic overhead is added in restructuring the deformed triangular grid, while in the latter case, due to the Monte-Carlo nature of the sampling, Poisson noise is introduced.

Due to the aforementioned drawbacks and due to significant progress in the development of accurate Eulerian schemes in the early 1980’s, Lagrangian methods employing quadrilateral or triangular grids have not been used extensively in multidimensional calculations of astrophysical flows (see Woodward76 , Pen98  as well as Oran\_Boris87  and the references therein for examples). On the other hand, the simplicity of SPH has made this method very popular for astrophysical (especially cosmological) simulations. Without attempting to escalate the very vigorous discussion, whether SPH or grid-based Eulerian schemes are to be prefered in astrophysical calculations (see e.g. Mueller98 ), we will argue below that, due to its Monte-Carlo nature, the SPH scheme appears to be rather unsuited for multidimensional nucleosynthesis calculations, especially in cases where hydrodynamic instabilities are known to be important.

Among Eulerian schemes, the so-called shock-capturing schemes have proven to be the most accurate ones for problems which involve discontinuities in the flow as shock waves (see WC84  for details). The latter are very frequently encountered in explosive events, since in these cases the flows can attain supersonic speeds. Shock-capturing schemes derive their accuracy from a discretization of the hydrodynamic equations which closely mimics the physics of compressible flows by making use of the Riemann problem, i.e. the dissolution of an arbitrary flow discontinuity into a set of simple waves (shocks, contact discontinuities and rarefaction waves). Suitably constructed Riemann problems at the interfaces between adjacent computational cells are solved within each time step, from which the complete solution of the system of conservation laws is constructed. This allows one to avoid the use of large amounts of artificial viscosity in order to obtain a well-behaved numerical scheme in the vicinity of shocks. One of the most accurate shock-capturing schemes, which has been widely applied in astrophysics, is the (direct Eulerian) PPM scheme of CW84 , a second order extension of Godunov’s original (and rather diffusive) first-order shock-capturing scheme Godunov . In addition to its accurate treatment of shocks PPM includes a special detection and steepening algorithm to minimize numerical diffusion across contact discontinuities.

The superiority of shock-capturing schemes in computing compressible flows has been demonstrated e.g. in WC84 , and their performance for computing reactive astrophysical flows was studied in FMA89  and Plewa\_Mueller99 . Fig. 1 shows a representative result from FMA89  in which PPM was compared to a number of older Eulerian schemes which were in wide-spread use until the mid 1980’s. The figure shows the width of a contact discontinuity as a function of the number of zones that it has travelled across a numerical grid. Most Eulerian schemes tend to smear such fluid (and also composition) interfaces without limit, i.e. the width of the “discontinuity” tends to grow with time. Of all the schemes investigated, only PPM maintained a sharp resolution of the interface within two zones. Still however, numerical diffusion cannot be completely avoided in Eulerian calculations and its minimization necessitates an adequate spatial resolution in addition to an excellent advection scheme. This has led to the development of adaptive mesh refinement methods BC89 , which concentrate the computational effort in critical regions of the flow and thereby often allow for substantial savings in computer time.

## Additional Physics

While the numerical problems encountered in solving the homogeneous part of the Euler equations are difficult to overcome, they represent only a part of the computational difficulties for a realistic simulation. The source terms, which are usually taken into account using the operator splitting technique LeVeque98 , often require much more computer time than the solution of the hydrodynamic equations themselves. This holds, e.g. if large nuclear networks need to be evolved with the hydrodynamics or transport processes need to be taken into account (as e.g. neutrino transport in core collapse supernovae, see Rampp\_Janka00  and the references therein and A. Burrows, this volume). In some cases even phenomenological (sub-grid) models might have to be introduced. This is e.g. the case for turbulent combustion in thermonuclear supernovae, where a white dwarf is incinerated by a deflagration front whose propagation speed is impossible to compute in a direct simulation since this would require a resolution of the turbulent energy cascade down to the dissipation length scale Reinecke+99 . Exacerbating the situation is the fact that stellar models, which serve as initial data for supernova simulations, might be affected by considerable uncertainties. In the absence of computational schemes and resources which allow for a consistent multidimensional treatment of stellar convection and rotation over stellar evolutionary time scales, one is forced to describe these phenomena by one-dimensional appoximations (see the contribution of N. Langer, this volume). Finally, uncertainties in nuclear reaction rates enter the calculations.

It is apparent that progress in only a single of the involved fields is not going to improve the accuracy of the desired nucleosynthetic yields considerably. In fact, a concerted effort in all areas appears to be required, since, as we will show in our example below, the different effects can conspire in falsifying the nucleosynthetic yields.

## Core Collapse Supernovae: A case study

Nucleosynthesis in core collapse supernovae is a good example for illustrating the aforementioned problems and we will start with a discussion of numerical diffusion using results of simple one-dimensional calculations. We subsequently address the complications introduced by convection in multidimensional calculations as well as by “additional physics”, i.e. neutronization due to neutrino matter interactions. Finally we show how a multidimensional numerical failure, the so-called “odd-even-decoupling” phenomenon, an instability which appears to plague most shock-capturing schemes and whose effects have not yet been discussed extensively in the numerical astrophysics literature, can enhance neutronization by strengthening hydrodynamic convection and affect the nucleosythetic yields in multidimensional simulations.

### Nucleosynthesis in a 15 $`\mathrm{M}_{}`$ star (1D)

In core collapse supernovae nucleosynthesis is triggered by a shock wave which forms after the collapse of the iron core of a massive star has proceeded to supranuclear densities. The shock, while initially powerful, stalls after a few milliseconds due to the energy losses from which it suffers while propagating through the outer iron core, but eventually ejects the outer stellar layers if heating by neutrinos from the collapsed core is able to overcompensate for the energy losses. In most nucleosynthesis calculations, however, these processes are not modelled in detail and instead a shock is initiated by simply depositing the typical observed supernova energy of $`10^{51}`$ erg near the center of a presupernova model.

In the left panels of Fig. 2 we show snapshots of the mass fractions from the first 500 ms of such a calculation, focusing on the silicon-rich layers just outside the iron core in which explosive nucleosynthesis takes place. Only 100 ms after the start of the calculations explosive silicon burning has frozen out and has left behind a significant abundance of $`{}_{}{}^{56}\mathrm{Ni}`$ as well as $`{}_{}{}^{40}\mathrm{Ca}`$, $`{}_{}{}^{36}\mathrm{Ar}`$ and $`{}_{}{}^{32}\mathrm{S}`$. Particularly noteworthy for the following discussion is the nucleus $`{}_{}{}^{44}\mathrm{Ti}`$. From one-dimensional Lagrangian nucleosynthesis calculations TNH96 , Woosley\_Weaver95  it is known that this isotope should be primarily synthesized in the innermost stellar layers which experience an $`\alpha `$-rich freezeout. However, in the present Eulerian calculation a significant abundance of $`{}_{}{}^{44}\mathrm{Ti}`$ has formed in zones with mass coordinates around 1.36 $`\mathrm{M}_{}`$ (marked with arrows in the left panels of Fig. 2), i.e. at the interface of the regions enriched in $`{}_{}{}^{4}\mathrm{He}`$ and $`{}_{}{}^{40}\mathrm{Ca}`$. This $`{}_{}{}^{44}\mathrm{Ti}`$ “bump” results from the reaction $`{}_{}{}^{40}\mathrm{Ca}(\alpha ,\gamma ){}_{}{}^{44}\mathrm{Ti}`$ and the amount of $`{}_{}{}^{44}\mathrm{Ti}`$ thereby produced is very sensitive to numerical diffusion in this region. Consequently, the strength of the $`{}_{}{}^{44}\mathrm{Ti}`$ “bump” varies with the diffusivity of the numerical scheme which is used to advect the nuclear species. This is illustrated in the right panels of Fig. 2. The top right panel depicts results from a calculation with Godunov’s first-order scheme. In this case, all mass fraction profiles are heavily smeared due to strong diffusion, as can be seen by a comparison to the middle right and bottom right panels which show results that were obtained with the FMA and CMA advection schemes of FMA89  and Plewa\_Mueller99 , respectively. Of all these schemes, CMA is the least diffusive since it is the only method which includes a detection and steepening algorithm for composition interfaces which was derived from PPM’s original detection and steepening algorithm for contact discontinuities. Note the size of the $`{}_{}{}^{44}\mathrm{Ti}`$ bump for the three different runs. The more diffusive schemes produce much more $`{}_{}{}^{44}\mathrm{Ti}`$. This is also illustrated in Fig. 3 which summarizes how the $`{}_{}{}^{44}\mathrm{Ti}`$ yield depends on the adopted advection scheme and the spatial resolution. While the CMA results (bottom curve) are already converged for a resolution of $`\mathrm{\Delta }r=40`$ km, FMA (middle curve) needs a resolution of about $`10`$ km. The first-order scheme (top curve) would need much finer zoning than $`\mathrm{\Delta }r=5`$ km to yield results of comparable quality. Note also that, if a diffusive advection scheme and coarse resolution are used, the errors might be as large as a factor of four!

It should be pointed out, however, that $`{}_{}{}^{44}\mathrm{Ti}`$ is a somewhat extreme (though very important) example. The (relative) errors due to numerical diffusion are usually smaller for the more abundant nuclei. This is illustrated in Fig. 4 which shows the dependence of the yields of different $`\alpha `$-nuclei on resolution in a 1D calculation from Kifonidis00 , in which no ad hoc energy deposition was adopted, but where the shock revival phase was followed in detail by including the effects of neutrino heating from a central light bulb neutrino source JM96 . It can be seen from this figure that individual elemental yields are more accurate than the $`{}_{}{}^{44}\mathrm{Ti}`$ yield by more than an order of magnitude. However, if one is aiming at a numerical accuracy of about 1% for all yields, about 3000 radial zones are required for this calculation. This amount of spatial resolution makes accurate multidimensional simulations very expensive.

### Convection and $`{}_{}{}^{56}\mathrm{Ni}`$ synthesis in core collapse supernovae

Neutrino matter interactions play a crucial role in the explosion of core collapse supernovae. They heat the matter behind the stalling shock and thereby trigger the explosion. On the other hand they determine the electron fraction per baryon, $`Y_\mathrm{e}`$, (or equivalently the ratio of protons to neutrons) in the ejecta and thus influence the nucleosynthetic yields. If the $`Y_\mathrm{e}`$ value of material that has been photodissociated by the shock is significantly reduced below 0.5 by neutrino/matter interactions, this matter will recombine mainly to neutron-rich nuclei after expanding and cooling. In that case nuclei with $`Z=N`$ like $`{}_{}{}^{56}\mathrm{Ni}`$ will not form in the ejecta. Fig. 5 which shows results of a two-dimensional simulation of shock revival illustrates this effect. In this calculation the luminosities of $`\nu _e`$ and $`\overline{\nu }_e`$ were such that $`\overline{\nu }_e`$ absorption on protons was favored against $`\nu _e`$ absorption on neutrons and thus the heated gas was neutronized JM96 . The sharp division between this material which is visible in the bubbles behind the shock in Fig. 5 and the lepton-rich material farther out which has formed $`{}_{}{}^{56}\mathrm{Ni}`$ (the region enclosed by the white contour line in Fig. 5) is clearly visible. The negative entropy gradient that the heating has imprinted on the layers between the radius of maximum neutrino heating (close to the center) and the shock farther out, has also led to strong convective motions: bubbles of heated deleptonized gas rise toward the shock while lower entropy flux tubes transport lepton-rich matter to deeper layers where it interacts with the neutrino fluxes much more efficiently. This interplay of convection and deleptonization is crucial for the $`{}_{}{}^{56}\mathrm{Ni}`$ yield. The shock is only able to heat a certain amount of lepton-rich material to temperatures in excess of the $`5\times 10^9`$ K which are required for $`{}_{}{}^{56}\mathrm{Ni}`$ synthesis. If convection is strong enough to advect significant amounts of this matter close to the neutron star, where the gas will experience deleptonization, the $`{}_{}{}^{56}\mathrm{Ni}`$ yield will be lower than in a model with no or weak convection.

### Multidimensional numerical failures: Odd-even decoupling

Quirk Quirk94  has reported on a subtle flaw in a number of shock capturing schemes which becomes evident when calculating multidimensional flows with strong, grid-aligned shocks. He has dubbed this failure the “odd-even decoupling” phenomenon. The problem shows up only if a sufficiently strong shock is either fully or nearly aligned with one of the coordinate directions of the grid, and if, in addition, the flow is slightly perturbed. This can be due to either perturbations intentionally introduced in order to study *physical* instabilities, as it is done in all studies of convection in supernovae, or due to perturbations caused by other flow features. Many Riemann solvers show the tendency to allow these perturbations to grow without limit along the shock surface, thus triggering a strong rippling of the shock front as well as the post shock state. In supernova simulations these perturbations, whose amplitudes can exceed those of the seed perturbations by several orders of magnitude, enhance the growth of hydrodynamic instabilities. In case of neutrino driven convection they lead to large-scale overturn and angular wavelengths of convective bubbles which are significantly larger than in a “clean” calculation. This artificial enhancement of convection is demonstrated in Fig. 6 where the entropy distribution of a simulation exhibiting odd-even decoupling (top panel) is compared to one in which the numerical failure has been cured (bottom panel). A modification of PPM’s original dissipation algorithms CW84  was necessary for this purpose. Alternatively the hybrid Riemann solver method of Quirk94  might be used. Note that the calculations have been carried out in spherical coordinates $`(r,\theta )`$ so that the (initially spherical) shock wave was fully aligned with the grid. However, cylindrical coordinates which have been used for plotting are indicated in the figures. The difference in the final $`{}_{}{}^{56}\mathrm{Ni}`$ yield for these two simulations was about 40%, the calculation not exhibiting odd-even decoupling showing the larger yield, as expected from our discussion in the previous section. This demonstrates that due to the strong coupling of neutrino physics and convection, the $`{}_{}{}^{56}\mathrm{Ni}`$ yield in multidimensional calculations is much more difficult to calculate correctly than in one spatial dimension. Therefore the error of at most a few percent which can be deduced for the latter case from Fig. 4 can be deceptive. The results shown in Fig. 6 also suggest that numerical noise, whatever its origin is, leads to a grossly overestimated efficiency of convection. Thus, schemes which are known to suffer from this problem (like SPH) do not appear to be suited for nucleosynthesis calculations in core collapse supernovae.

## Conclusions

Realistic nucleosynthesis calculations in astrophysical contexts represent a challenge in many respects. The difficulties involve the numerical treatment of multidimensional hydrodynamic advection, complex physics in addition to hydrodynamics and burning, disparate length and time scales, realistic initial conditions and uncertainties in reaction rates. It is our conviction, that substantial efforts are required in each of these fields in order to obtain reliable yields in multidimensional nucleosynthesis calculations.


text:
# Stellar Kinematics of the Spiral Structure in the Solar Neighborhood

## 1 Introduction

Spiral structures of galaxies have been studied for a long time in order to understand how these structures are formed (e.g., Roberts, Roberts & Shu 1975; Rohlfs 1977; Binney & Tremaine 1987). One of the proposed models to explain spiral arms is that they are just material arms, where the stars originally making up a spiral arm remains in the arm even at the later time. However this simple model holds a wellknown problem which is called “winding problem”: the differential rotation in galactic disks winds up the arm in a short time compared with the age of galaxies, so that the spiral pattern would be too tightly wound compared with the observed spiral structures. In contrast, the currently most popular model, which is free from the winding problem, is the density wave theory (Lin & Shu 1964), where a spiral arm is regarded as a wave and wavelike oscillation of stellar motions propagates through galactic disks. In this picture, the global spiral pattern is sustained independently of individual stars moving at different angular velocities. For a comprehensive review of the density wave theory, see, e.g., Rohlfs (1977) and Binney & Tremaine (1987).

The density wave theory has been suggested by various observational aspects in spiral galaxies, including the relative distributions of dust lane, interstellar gas, and H II regions across the arms (Fujimoto 1968; Roberts 1969; Rohlfs 1977), intensity distribution of radio continuum radiation (Mathewson, van der Kruit, & Brouw 1972), and systematic variation of gaseous velocity fields near the arms (e.g., Visser 1980). In particular, recent high-resolution observations using CO emission have revealed detailed streaming motions of molecular gas, which are generally in agreement with predictions of the density wave models (e.g., Kuno & Nakai 1997; Aalto et al. 1999). However, we note that these observational results provide only an outcome of nonlinear interaction between interstellar matter and background stellar arms, and it is yet unknown whether the motions of stars themselves, which make up spiral pattern, actually follow those predicted by the density wave theory. In this regard, the direct access to detailed stellar kinematics in disks is possible only in our Galaxy.

Here, we present a method to clarify this issue, based on the analysis of local kinematics of disk stars. We then apply the method to 78 Cepheids in the solar neighborhood, for which the precise data of proper motions are available from the Hipparcos Catalogue (ESA 1997). Also, the distances to these sample stars can be accurately estimated from the period-luminosity relation, so that combined with the radial velocity data, the full three-dimensional velocities are available. We note here that although we focus on the local kinematics of spiral arms in this work, the method we develop here can be applied to the motions of more remote stars distributed over a whole disk, for which precise astrometric data will be provided by the next-generation satellites such as FAME and GAIA.

Our paper is organized as follows. In §2, we describe the method to determine whether or not the motions of stars agree with those expected in the density wave theory. In §3, we show the detail of the sample stars and the fundamental parameters of our Galaxy adopted in this work. The application of our method to the sample stars is shown in §4. Finally, §5 is devoted to discussion and conclusions.

## 2 METHOD

All orbits of the disk stars in our Galaxy are not perfectly circular. The discrepancy between the motion of a star and its reference circular orbit, called orbit of the guiding center, can be represented by an epicyclic motion of the star (Binney & Tremaine 1987). We describe here a method to compare the epicyclic motions of the stars observed near the Sun with those expected in the density wave theory.

First, we define the “position phase”, $`\chi `$, as a function of the position of a star in the Galactic plane, by assuming that the shape of the spiral arm is logarithmic as follows:

$$\chi =\mathrm{ln}\frac{(r/r_g)}{\mathrm{tan}i}\theta +\mathrm{\Omega }_pt,$$
(1)
where $`r`$ is the distance between the star and the Galactic Center (GC), $`r_g`$ is the distance between the Sun and GC, $`i`$ is the pitch angle of the spiral arm, $`\theta `$ is the angle between $`r`$ and $`r_g`$, and $`\mathrm{\Omega }_p`$ is the angular velocity of the spiral pattern in the azimuthal direction. We note that all the stars on the same arm have the same value of $`\chi `$ (Figure 1).

Second, we define the phase of the epicyclic motion of a star, $`\varphi `$, which is the angle between the velocity vector of the epicyclic motion and the direction perpendicular to a spiral arm (Figure 1). Then, if the stellar motions in the arm obey the density wave theory, all the stars on the same arm have the same value of $`\varphi `$ as well as the same value of $`\chi `$. Also, when the number of arms is $`m`$, the epicyclic frequency $`\kappa `$ of the stellar motions is equal to $`m`$ times as large as the angular frequency $`\mathrm{\Omega }\mathrm{\Omega }_p`$ \[$`\kappa =m(\mathrm{\Omega }\mathrm{\Omega }_p)`$\]. If so, since $`d\varphi /dt=\kappa `$, and $`d\chi /dt=\mathrm{\Omega }_pd\theta /dt=\mathrm{\Omega }_p\mathrm{\Omega }`$ from equation (1), we obtain the following linear relation between $`\varphi `$ and $`\chi `$,

$$\varphi m\chi (mod2\pi ),$$
(2)
if the stellar motions obey the density wave theory.

Thus, by assessing this linear relation between observed $`\varphi `$ and $`\chi `$, it is possible to investigate whether or not the observed stellar motions follow those predicted by the density way theory. We can also derive the number of the arms, $`m`$, from the slope in the $`\varphi \chi `$ relation. In Figure 3, we show the case $`m=4`$ (solid lines) as an example.

## 3 DATA

We adopt the Galactic Cepheids as the tracers of stellar motions in spiral arms, because these bright stars can be seen from a long distance and so we can investigate a wide region around the Sun. Also, these young populations show only a small deviation from circular rotation, thereby allowing us to analyze the internal kinematics of spiral arms alone; relatively old stars, such as dwarf stars, have too large velocity dispersions, possibly due to repeated gravitational interactions with massive clouds (Spitzer & Schwarzschild 1953) in addition to the effect of spiral arms. Furthermore we can obtain accurate distances for Cepheids, based on the relation between pulsation period and absolute magnitude. We do not use the Hipparcos parallaxes, which have generally large errors for many Cepheids beyond $`100`$ pc from the Sun. We adopt the Cepheid catalog compiled by Mishurov et al. (1997) for distances and radial velocities, and the Hipparcos Catalogue for proper motions, to calculate the individual motions of the Cepheids in the Galactic plane.

In the Mishurov et al. catalog, Cepheids in the region of $`r>4`$ kpc and those in a binary system are excluded. Also, nearby Cepheids in the region of $`r<0.5`$ kpc are excluded in order to reduce the local effects like Gould’s Belt. Furthermore, Cepheids whose pulsation periods exceed 9 days are also excluded, because they are supposed to be extremely young objects. In addition, we further exclude the Cepheids within 1 kpc for the following reason. When we compare the observed distribution of Cepheids in the $`\varphi \chi `$ plane with the theoretical prediction in a quantitative manner, as will be described later, we assign larger statistical weights to the Cepheids having a smaller observational error \[eq. (3)\]. This leads to larger weights to the Cepheids located close to the Sun, say $`r<1`$ kpc, so that the result of the analysis will be largely determined by only small number of Cepheids in a small region near the Sun. We also exclude two Cepheids whose peculiar velocity is exceptionally large, over 50 km s-1, compared to other ones. As a consequence, we adopt 78 Cepheids in this work, and their spatial distribution is shown in Figure 2.

We adopt $`r_g=8.3`$ kpc in our analysis, which is approximately an average of observed values ranging from 8.1 kpc to 8.5 kpc (Kerr & Lynden-Bell 1986; Hanson 1987; Pont, Mayor & Burk 1994; Feast & Whitelock 1997). As for the pitch angle of spiral arms, $`i`$, in our Galaxy, several authors investigated one of the conspicuous spiral arms near the Sun, the Sagittarius arm, based on the spatial distributions of open clusters, CO emissions, or O-B2 clusters, and arrived at several to about $`20^{}`$ (Pavlovskaya & Suchkov 1984; Dame et al. 1986; Grabelsky et al. 1988; Alfaro, Cabrera-Cano & Delgado 1992). In this work, we adopt $`i=8.0^{}`$ in the solar neighborhood, where the negative value for $`i`$ denotes a trailing arm. We have found that even if we change the values of these parameters over a likely range, the result shown below remains essentially unchanged.

In order to analyze the epicyclic motion of the Cepheids, we require to subtract both the local solar motion with respect to the local standard of rest and the effect of the Galactic differential rotation from the observed motions. As the local solar motion, we adopt 15.5 km s-1 in the direction to $`l_{\mathrm{}}`$=45 and $`b_{\mathrm{}}=23.6^{}`$ (Kulikovskij 1985). To estimate the effect of the Galactic differential rotation, we assume that in the concerned region within about $`4`$ kpc from the Sun, the Galactic rotation is monotonously changing with distance from the Sun, where the Oort constant $`A`$ is assumed to be $`13.0`$ km s-1kpc-1. We again note that the result shown in §4 is essentially independent of the value of $`A`$ or the assumption for differential rotation.

## 4 RESULTS

We present the relation between the phases of the epicyclic motions of the Cepheids $`\varphi `$ and their position phases $`\chi `$ in Figure 3. We also show the case for the 4-armed galaxy which obeys the density wave theory (solid lines). If the motions of the Cepheids in the solar neighborhood are in accordance with those expected in the density wave theory, we expect the similar linear relation between $`\varphi `$ and $`\chi `$. However, such a linear relation is not apparent for the Cepheids, which may imply that these objects in the region we investigate do not obey the density wave theory.

To be more quantitative, we analyze the motions of the Cepheids by using the sum of squares of the deviations for the observed $`\varphi `$’s from the theoretically expected ones. We define this deviation, $`\mathrm{\Delta }^2`$, as follows:

$$\mathrm{\Delta }^2\underset{i}{}\frac{\delta \varphi _i^2}{\sigma _i^2},$$
(3)
where $`\delta \varphi _i`$ is defined as,

$$\delta \varphi _imin(|\varphi _i\varphi _{DW}(\chi _i)|,2\pi |\varphi _i\varphi _{DW}(\chi _i)|),$$
(4)
where $`\sigma _i`$ is the error in the phase $`\varphi _i`$ of the star $`i`$, and $`\varphi _{DW}(\chi _i)`$ is the phase when the motion of the star $`i`$ obeys the density wave theory. Here we note that the error $`\sigma _i`$ includes both the observational error and the velocity dispersion of the Cepheids. We adopt 13 km s-1 for the velocity dispersion. The quantity $`\delta \varphi _i`$ denotes the discrepancy between $`\varphi _i`$ and $`\varphi _{DW}(\chi _i)`$, for which we take the smaller value of $`|\varphi _i\varphi _{DW}(\chi _i)|`$ and $`2\pi |\varphi _i\varphi _{DW}(\chi _i)|`$, because $`\varphi _i`$ has a period of $`2\pi `$. The phase $`\varphi _{DW}(\chi _i)`$ is determined by minimizing the deviation $`\mathrm{\Delta }^2`$ for each value of the arm number $`m`$.

If the motions of stars are in accordance with those expected in the density wave theory, the expected value of the deviation $`\mathrm{\Delta }^2`$ is about the same as the number of the sample stars (78 in this work), because the distribution of $`\delta \varphi _i^2`$ is normalized by its standard dispersion $`\sigma _i^2`$. Thus, in this case, we will obtain a minimum of $`\mathrm{\Delta }^2`$ with $`\mathrm{\Delta }_{min}^280`$ at a specific arm number $`m=m_g`$, whereas at other $`m`$’s, the value of $`\mathrm{\Delta }^2`$ will be systematically larger than 80; the larger ratio of $`\mathrm{\Delta }^2/\mathrm{\Delta }_{min}^2`$ implies more likely the observed stellar motions match the density wave theory. On the other hand, if $`\mathrm{\Delta }^2`$ is always larger than 80 at all $`m`$’s, we may conclude that the stellar motions are totally inconsistent with those expected in the theory.

We plot the deviation $`\mathrm{\Delta }^2`$ as a function of $`m`$ in Figure 4. Solid line shows the result for our Cepheid sample. It is apparent that $`\mathrm{\Delta }^2`$ is around $`250`$ for all values of the parameter $`m`$, without showing a noticeable minimum with $`\mathrm{\Delta }^280`$.

In order to examine the significance of this result, we investigate three hypothetical models for comparison. As the first model, we randomly select $`\varphi `$, independently of $`\chi `$, so that there is no correlation between $`\varphi `$ and $`\chi `$. We then assign the same errors $`\sigma _i`$ as those for the Cepheids to each $`\varphi `$. Thick dashed line in Figure 4 shows $`\mathrm{\Delta }^2`$ vs. $`m`$ derived from this “random model”. It follows that the properties of $`\mathrm{\Delta }^2`$ as a function of $`m`$ are basically the same as for the Cepheids, thereby implying that the phase $`\varphi `$ of the Cepheids is random in the solar neighborhood. The second model, as shown by thin dashed line, follows the density wave theory with $`m=4`$. We assign the same errors $`\sigma _i`$ as those for the Cepheids. The value of $`\mathrm{\Delta }^2`$ is around $`250`$, except for the case of $`m=4`$ at the minimum of $`\mathrm{\Delta }^2`$ ($`\mathrm{\Delta }_{min}^280`$). This behavior of $`\mathrm{\Delta }^2`$ is in sharp contrast to the case of the Cepheids. The last model (dotted line) is similar to the second one, except for the assignment of larger errors $`\sigma _i`$: 1.4 times as much as those of the Cepheids, in order to see the effect of $`\sigma _i`$ on the result. In this case, the value $`\mathrm{\Delta }^2`$ is smaller than other cases: $`\mathrm{\Delta }150`$ except for the case of $`m=4`$ at the minimum of $`\mathrm{\Delta }^2`$ ($`\mathrm{\Delta }_{min}^280`$). Thus, even if the errors $`\sigma _i`$ are large, which gives rise to small $`\mathrm{\Delta }^2/\mathrm{\Delta }_{min}^2`$, $`\mathrm{\Delta }^2`$ has a noticeable minimum at a specific value of $`m`$, in sharp contrast to the case of the Cepheids.

Therefore, these results suggest us that the motions of the Cepheids in the solar neighborhood are inconsistent with those expected in the density wave theory, whereas the random model, in which the phases $`\varphi `$ and $`\chi `$ are randomly selected, reproduces the observation in a reasonable manner.

## 5 DISCUSSION AND CONCLUSIONS

We have presented a method, based on the analysis of local kinematics of disk stars, to clarify whether the motions of the stars, which make up the spiral arms, follow those expected in the density wave theory. The method utilizes the comparison with the expected linear relation between the “position phases” of the stars $`\chi `$ and those of their epicyclic motions $`\varphi `$, as given in equation (2). The application of the method to the 78 Galactic Cepheids within 4 kpc from the Sun, for which accurate proper motions are available from Hipparcos, has revealed that the relation between $`\chi `$ and $`\varphi `$ for the Cepheids does not show the expected linear relation. Based on the quantitative analysis using the deviation $`\mathrm{\Delta }^2`$, we conclude that the observed motions of the Cepheids are well reproduced by the random model having no correlation between $`\chi `$ and $`\varphi `$.

There are a couple of possibilities to explain the current results, even if the spiral arms follow the density wave theory. First, the spiral structure around the Sun may not be simple as given in equation (1). In fact, many of our sample Cepheids belong to the local spiral structure called the Orion arm, for which the definite conclusion on its spatial structure is yet to be reached. It is frequently expressed as “Orion spur”, having a rather irregular pattern compared to other large-scale arms, Sagittarius and Perseus arm (see, e.g., Gilmore, King, & van der Kruit 1989). The existence of the Orion arm may give disturbances on the density wave motions of stars induced by these large-scale arms. Second, the Cepheids we have adopted here may still convey systematic velocities of dens gas clouds from which these stars were formed, in the form of the streaming motions. If there still exist some individual streaming motions among the sample stars, such motions may violate the ideal linear relation between $`\chi `$ and $`\varphi `$ expected for the density wave motions. Third, the Cepheids we have adopted here may have already experienced some scattering by dens gas clouds, thus having large velocity dispersions (Spitzer & Schwarzschild 1953). However, our experiment in §4 (dotted line in Figure 4) implies that the effect of velocity dispersions of our sample on the result appears to be minor.

In order to settle the last issue described above more clearly, we have repeated our analysis using younger populations with smaller velocity dispersions than the Cepheids. As such young stars, we have adopted the O-B5 stars, although due to their fainter luminosities than Cepheids, the sample with available proper motions is confined to the narrower region near the Sun. These sample stars are taken from the NASA SKY2000 Master Star Catalog Ver. 2 (Sande et al. 1998) which provides almost 300,000 stars brighter than 8 mag. The catalog contains many basic quantities, such as MK classification, luminosity class, apparent magnitude, color, radial velocity, and so on. We have then calibrated distances using Hipparcos parallaxes or spectroscopic distances using the program kindly supplied by Drs. M. Sôma and M. Yoshizawa, and also obtained accurate proper motions by the cross-identification with the Hipparcos and ACT Reference Catalogs (Urban, Corbin, & Wycoff 1998). After removing binaries and multiples, we have selected 773 O-B5 stars for which full three-dimensional velocities are available. Then, the application of the method we have developed here has revealed that the deviation $`\mathrm{\Delta }^2`$ as a function of $`m`$ remains essentially constant of the order of 3000, without showing any noticeable minimum at a specific value of $`m`$. Thus, even the motions of such young populations with small velocity dispersions are in contradiction with the density wave theory. We note here that most of these O-B5 stars are located within $`1`$ kpc from the Sun, so the effect of the local irregular spiral on the result cannot be negligible.

More definite conclusions on the issue we have addressed here require the assembly and analysis of much larger numbers of stars with accurate distances and proper motions, so that the statistical fluctuation in the result can be significantly reduced. Also, it is necessary to assemble the data of more remote stars over a large fraction of the disk, thereby diminishing effects of local irregular spiral structures on the kinematic analysis. Indeed, next-generation satellites such as FAME and GAIA will provide very precise astrometric data for huge numbers of the Galactic stars, and will thus offer us an opportunity to assess detailed motions of disk stars in conjunction with the density wave theory.

We are grateful to M. Sôma and M. Yoshizawa for providing us the program for calculating spectroscopic distances. T.Y. would like to thank K. Okoshi and S. Bouquillon for useful discussion. This work was supported in part by Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists (No.00074).