text:
# Mapping Deconfinement with a Compact Star Phase Diagram

## 1 Introduction

At present, the existence of exotic phases of matter at high densities is under experimental investigation in ultrarelativistic heavy-ion collisions  the most prominent being the deconfined phase of QCD . While the diagnostics of a phase transition in experiments with heavy-ion beams faces the problems of strong nonequilibrium and finite size, the dense matter in a compact star forms a macroscopic system in thermal and chemical equilibrium for which effects signalling a phase transition shall be most pronounced.

We introduce a classification of rotating compact stars in the plane of their angular frequency $`\mathrm{\Omega }`$ and mass (baryon number $`N`$) which we will call phase diagram. In this diagram, configurations with high density matter cores are separated from conventional ones by a critical phase transition line. The position and the form of these lines are sensitive to changes in the equation of state of stellar matter . Since the phase diagram of rotating compact objects seems to be a more general approach for investigations of phase transition effects in the interior of the star we assume that the deconfinement transition could be a particular case besides of other possibilities for phase transitions like pion or kaon condensation as discussed, e.g. in . Therefore, our aim is to suggest it as a heuristic tool for obtaining constraints on the EoS at high densities from the rotational behaviour of compact stars.

## 2 EoS and compact star phase diagram

Our focus is on the elucidation of qualitative features of signals from the high density phase transition in the pulsar timing, therefore we use a generic form of an equation of state (EoS) with such a transition. We use the polytropic type equation of state  for different values of the incompressibility  $`K_{L,H}(n)=9dP/dn`$ at the saturation density, see Ref. . The phase transition between the lower and higher density phases is made by the Maxwell construction  and compared to a relativistic mean field model consisting of a linear Walecka plus dynamical quark model EoS  with a Gibbs construction , see Fig.2.

In Fig. 2 we display the phase diagrams for the rotating star configurations, which correspond to the three model EoS of Fig. 2. These phase diagrams have four regions: (i) the region above the maximum frequency $`\mathrm{\Omega }>\mathrm{\Omega }_\mathrm{K}(N)`$ where no stationary rotating configurations are found, (ii) the region of black holes $`N>N_{\mathrm{max}}(\mathrm{\Omega })`$, and the region of stable compact stars which is subdivided by the critical line $`N_{\mathrm{crit}}(\mathrm{\Omega })`$ into (iii) the region of hybrid stars for $`N>N_{\mathrm{crit}}(\mathrm{\Omega })`$ where configurations contain a core with a second, high density phase and (iv) the region of mono-phase stars without such a core.

From the comparison of the regional structure of these three different phase diagrams in Fig. 2 with the corresponding EoS in Fig. 2 we arrive at the main result of this paper that there are the following correlations between the topology of the lines $`N_{\mathrm{max}}(\mathrm{\Omega })`$ and $`N_{\mathrm{crit}}(\mathrm{\Omega })`$ and the properties of two-phase EoS:

* The hardness of the high density EoS determines the maximum mass of the star, which is given by the line $`N_{\mathrm{max}}(\mathrm{\Omega })`$. Therefore $`N_{\mathrm{max}}(0)`$ is proportional to the parameter $`K_H(n_H)`$, where $`n_H`$ is the density of the transition to high density phase.
* The onset of the phase transition line $`N_{\mathrm{crit}}(0)`$ depends on the density $`n_H`$ and $`K_L(n_L)`$ where $`n_L`$ is the density of the transition to the low density phase.
* The curvature of the lines $`N_{\mathrm{max}}(\mathrm{\Omega })`$ and $`N_{\mathrm{crit}}(\mathrm{\Omega })`$ is proportional to the compressibility of the high and low density phases, respectively.

Therefore, a verification of the existence of the critical lines $`N_{\mathrm{crit}}(\mathrm{\Omega })`$ and $`N_{\mathrm{max}}(\mathrm{\Omega })`$ by observation of the rotational behavior of compact objects would constrain the parameters of the EoS for neutron star matter. Moreover, we have investigated different trajectories of rotating compact star evolution in the phase diagram in order to identify scenarios, which result in signatures of the deconfinement phase transition. To prove that the appearance or disappearance of a high density phase during the rotational evolution of the star could entail observational consequences for the angular velocity we consider three main representatives different classes of trajectories which could cross the critical line on phase diagram. These classes of tracks are: (a) spindown of isolated (non-accreting, $`\dot{N}=0`$) pulsars due to magnetic dipole radiation , (b) spin up in accreting systems with weak magnetic field  ($`N\mathrm{const}`$, vertical tracks) and (c) accretion either with strong magnetic field  or for accreting binaries emitting gravitational waves , for which $`\mathrm{\Omega }\mathrm{const}`$ (horizontal tracks). In the case of (a) the spindown $`K=K_{\mathrm{rad}}`$ or (b) spinup $`K=K_{\mathrm{acc}}`$ evolutions (in both cases $`K_{\mathrm{int}}<<K_{\mathrm{ext}}`$) the objects can undergo a phase transition if the baryon number lies within the interval $`N_{\mathrm{crit}}(\mathrm{\Omega }=0)<N<N_\mathrm{c}`$, where $`N_\mathrm{c}`$ is the end piont of the critical line $`N_{\mathrm{crit}}(\mathrm{\Omega })`$. If the core of compact star is soft enough (as in case (SH-S)) the critical line $`N_{\mathrm{crit}}(\mathrm{\Omega })`$ crosses $`N_{\mathrm{max}}(\mathrm{\Omega })`$ at some $`\mathrm{\Omega }_c(N_c)<\mathrm{\Omega }_K(N_c)`$. This means that massive mono-phase configurations with total baryon number $`N>N_c`$ rotating with angular velocities in the interval $`\mathrm{\Omega }_c<\mathrm{\Omega }<\mathrm{\Omega }_K`$ should encounter a transition to a black hole during the spin down evolution. In those cases when the core EoS is harder (H-S-H and SH-H), the region of hybrid stars is a band which separates mono-phase configurations from black holes, see upper two panels of Fig. 2. As it has been shown in  for the vertical tracks (a) and (b) in the phase diagram, the braking index $`n(\mathrm{\Omega })`$ changes its value from $`n(\mathrm{\Omega })>3`$ in the region (iii) to $`n(\mathrm{\Omega })<3`$ in (iv). This is the braking index signal for a deconfinement transition introduced in Ref. . The third evolutionary track is accretion with strong magnetic fields  and/or gravitational wave emission (horizontal tracks) .For this case the $`\dot{\mathrm{\Omega }}`$ first decreases as long as the moment of inertia monotonously increases with $`N`$. When passing the critical line $`N_{\mathrm{crit}}(\mathrm{\Omega })`$ for the phase transition, the moment of inertia starts decreasing and the internal torque term $`K_{\mathrm{int}}`$ changes sign. This leads to a narrow dip for $`\dot{\mathrm{\Omega }}(N)`$ in the vicinity of this line. As a result, the phase diagram gets overpopulated for $`N\stackrel{<}{}N_{\mathrm{crit}}(\mathrm{\Omega })`$ and depopulated for $`N\stackrel{>}{}N_{\mathrm{crit}}(\mathrm{\Omega })`$ up to the second maximum of $`I(N,\mathrm{\Omega })`$ close to the black-hole line $`N_{\mathrm{max}}(\mathrm{\Omega })`$. A population gap in the phase diagram of compact stars is appears as a detectable indicator for hybrid star configurations.

## 3 Conclusion

We have shown that population clustering of compact stars at the phase transition line could be a signal for the occurrence of stars with high density matter cores and a measure for obtaining constraints on the EoS at high densities. In the case without a phase transition, the moment of inertia could at best saturate before the transition to the black hole region and consequently $`\dot{\mathrm{\Omega }}`$ would also saturate. This would entail a smooth population of the phase diagram without a pronounced structure . The so called Z sources of QPOs in LMXBs are suggested as objects which should predominantly populate the region of the suspected phase border between hadronic stars and quark core stars in phase diagram . The existence of a population gap between the critical line and the black hole limit on the phase diagram is suggested as a signal for a phase transition at supernuclear densities in accreting compact stars . The absence of a such gap does not exclude the existence of deconfined quark matter in LMXBs as we have shown in a study of accreting strange stars .


format: text

text:
# Can the Chaplygin gas be a plausible model for dark energy?

## Abstract

In this note two cosmological models representing the flat Friedmann Universe filled with a Chaplygin fluid, with or without dust, are analyzed in terms of the recently proposed ”statefinder” parameters . Trajectories of both models in the parameter plane are shown to be significantly different w.r.t. ”quiessence” and ”tracker” models. The generalized Chaplygin gas model with an equation of state of the form $`p=A/\rho ^\alpha `$ is also analyzed in terms of the statefinder parameters.
PACS numbers: 98.80.Es, 98.80.Cq, 98.80.Hw

In the search for cosmological models describing the observed cosmic acceleration , the inspiration coming from inflation has suggested mainly models making use of scalar fields . There are of course alternatives; in particular, in  an elementary model has been presented describing a Friedmann universe filled with a perfect fluid obeying the Chaplygin equation of state

$$p=\frac{A}{\rho },$$
(1)
where $`A`$ is a positive constant (for a thorough review see Ref. ). The interesting feature of this model is that it naturally provides a universe that undergoes a transition from a decelerating phase, driven by dust-like matter, to a cosmic acceleration at later stages of its evolution (see  for details). An interesting attempt to justify this model  makes use of an effective field theory for a three-brane universe .

In the flat case, the model can be equivalently described in terms of a homogeneous minimally coupled scalar field $`\varphi `$, with potential

$$V(\varphi )=\frac{1}{2}\sqrt{A}\left(\mathrm{cosh}3\varphi +\frac{1}{\mathrm{cosh}3\varphi }\right).$$
(2)
However, since models trying to provide a description (if not an explanation) of the cosmic acceleration are proliferating, there exists the problem of discriminating between the various contenders. To this aim a new proposal introduced in  makes use of a pair of parameters $`\{r,s\}`$, called “statefinder”. The relevant definition is as follows:

$$r\frac{\stackrel{\mathrm{}}{a}}{aH^3},s\frac{r1}{3(q1/2)},$$
(3)
where $`H\frac{\dot{a}}{a}`$ is the Hubble constant and $`q\frac{\ddot{a}}{aH^2}`$ is the deceleration parameter. The new feature of the statefinder is that it involves the third derivative of the cosmological radius.

Trajectories in the $`\{s,r\}`$-plane corresponding to different cosmological models exhibit qualitatively different behaviors. $`\mathrm{\Lambda }`$CDM model diagrams correspond to the fixed point $`s=0`$, $`r=1`$. The so-called ”quiessence” models  are described by vertical segments with $`r`$ decreasing from $`r=1`$ down to some definite value. Tracker models  have typical trajectories similar to arcs of parabola lying in the positive quadrant with positive second derivative.

The current location of the parameters $`s`$ and $`r`$ in these diagrams can be calculated in models (given the deceleration parameter); it may also be extracted from data coming from SNAP (SuperNovae Acceleration Probe)-type experiments . Therefore, the statefinder diagnostic combined with future SNAP observations may possibly be used to discriminate between different dark energy models.

In this short paper we apply that diagnostic to Chaplygin cosmological models (the direct comparison of the already available supernova data with the Chaplygin gas model was also undertaken recently in ). We consider both the one-fluid pure Chaplygin gas model and a two-fluid model where dust is present as well. We show that these models are different from those considered in  and they are worth of further study.

To begin with, let us rewrite the formulae for the statefinder parameters  in a form convenient for our purposes. We shall need the Friedmann equation for the flat universe

$$H^2=\frac{\dot{a}^2}{a^2}=\rho $$
(4)
and the energy conservation equation

$$\dot{\rho }=3H(\rho +p).$$
(5)
Using these two equations it is easy to find that

$$q=\frac{1}{2}+\frac{3}{2}\frac{p}{\rho },$$
(6)
and then

$$r=1\frac{3\dot{p}}{2\rho \sqrt{\rho }},s=\frac{\dot{p}}{3p\sqrt{\rho }}.$$
(7)
For a one-component fluid111We confine ourselves to the case, which is fulfilled in our models, of a fluid for which the equation of state has the form $`p=p(\rho )`$. these formulae become especially simple. Since

$$\dot{p}=\frac{p}{\rho }\dot{\rho }=3\sqrt{\rho }(\rho +p)\frac{p}{\rho },$$
(8)
we easily get:

$$r=1+\frac{9}{2}\left(1+\frac{p}{\rho }\right)\frac{p}{\rho },s=\left(1+\frac{\rho }{p}\right)\frac{p}{\rho }.$$
(9)
For the Chaplygin gas one has simply that

$$v_s^2=\frac{p}{\rho }=\frac{A}{\rho ^2}=\frac{p}{\rho }=1+s$$
(10)
($`v_s^2`$ is the square of the velocity of sound) and therefore

$$r=1\frac{9}{2}s(1+s).$$
(11)
Thus, the curve of $`r(s)`$ is an arc of parabola. To find the admissible values of $`s`$, we note that Eqs. (1), (4) and (5) easily give the following dependence of the energy density on the cosmological scale factor :

$$\rho =\sqrt{A+\frac{B}{a^6}},$$
(12)
where $`B`$ is an integration constant; therefore

$$v_s^2=1+s=\frac{A}{A+\frac{B}{a^6}}.$$
(13)
When the cosmological scale factor $`a`$ changes from $`0`$ to $`\mathrm{}`$ the velocity of sound varies from $`0`$ to $`1`$ and $`s`$ varies from $`1`$ to $`0`$. Thus in our model the statefinder $`s`$ takes negative values; this feature is not shared by quiessence and tracker models considered in .

As $`s`$ varies in the interval $`[1,0]`$, $`r`$ first increases from $`r=1`$ to its maximum value and then decreases to the $`\mathrm{\Lambda }`$CDM fixed point $`s=0`$, $`r=1`$ (see Fig. 1).

If $`q0.5`$ the current values of the statefinder (within our model) are $`s0.3,r1.9`$. In  it is reported an interesting numerical experiment based on 1000 realizations of a SNAP-type experiment, probing a fiducial $`\mathrm{\Lambda }`$CDM model. Our values of the statefinder lie outside the three-sigma confidence region displayed in . Based on this fact it can be expected that future SNAP experiments should be able to discriminate between the pure Chaplygin gas model and the standard $`\mathrm{\Lambda }`$CDM model.

Let us consider now a more ”realistic” cosmological model which, besides a Chaplygin’s component, contains also a dust component. For a two-component fluid Eqs. (7) take the following form:

$$r=1+\frac{9}{2(\rho +\rho _1)}\left[\frac{p}{\rho }(\rho +p)+\frac{p_1}{\rho _1}(\rho _1+p_1)\right],$$
(14)
$$s=\frac{1}{p+p_1}\left[\frac{p}{\rho }(\rho +p)+\frac{p_1}{\rho _1}(\rho _1+p_1)\right].$$
(15)
If one of the fluids is dust, i.e. $`p_1=p_d=0`$,the above formulae become

$$r=1+\frac{9(\rho +p)}{2(\rho +\rho _d)}\frac{p}{\rho },s=\frac{\rho +p}{p}\frac{p}{\rho }.$$
(16)
If the second fluid is the Chaplygin gas, proceeding exactly as before we obtain the following relation:

$$r=1\frac{9}{2}\frac{s(s+1)}{1+\frac{\rho _d}{\rho }}.$$
(17)
To find the term $`\rho _d/\rho `$ we write down the dependence of the dust density on the cosmological scale factor:

$$\rho _d=\frac{C}{a^3},$$
(18)
where $`C`$ is a positive constant. Eq. (13) gives that $`Aa^6+B=\frac{B}{s}`$ and therefore

$$\frac{\rho _d}{\rho }=\frac{C}{\sqrt{Aa^6+B}}=\kappa \sqrt{s},$$
(19)
where the constant $`\kappa =C/\sqrt{B}`$ is the ratio between the energy densities of dust and of the Chaplygin gas at the beginning of the cosmological evolution (cf. Eq. (12)). Thus

$$r=1\frac{9}{2}\frac{s(s+1)}{1+\kappa \sqrt{s}}.$$
(20)
Graphs of the function (20) for different choices of $`\kappa `$ are plotted in Fig. 2.

In this case there are choices of the parameters so that the current values of the statefinder are close to the $`\mathrm{\Lambda }`$CDM fixed point; indeed for $`\kappa =1`$ we have $`s=0.09`$ and $`r=1.2835`$. By increasing $`\kappa `$ we get closer and closer. Already for $`\kappa =2`$ we get $`s=0.035`$, $`r=1.11`$ while for $`\kappa >10`$ the statefinder essentially coincides with the $`\mathrm{\Lambda }`$CDM fixed point (see Fig. 2).

Thus our two-fluid cosmological models (with $`\kappa `$ say bigger than 2) cannot be discriminated from the $`\mathrm{\Lambda }`$CDM model on the basis of the statefinder analysis. On the other hand this fact could also be interpreted as an argument in favor of the model because the $`\mathrm{\Lambda }`$CDM model does not contradict observations up to now.

The attractive point of our two-fluid cosmological model with a parameter $`\kappa `$ of order one is that it may suggest a solution to the cosmic coincidence conundrum: indeed here the initial values of the energies of dust and of the Chaplygin gas are of the same order of magnitude.

The comparison of the two-fluid (i.e the Chaplygin gas plus dust) cosmological model with observational data has also been studied in . An analysis of the data obtained from 26 Supernovae leads to the conclusion  that the best fitting model seems to be the pure Chaplygin gas without dust. On the other hand  by combining the results of the observations of 92 Supernovae with the matter power spectrum a wide range of values of $`\kappa `$ seems to be admissible.

Thus, for getting more precise constraints on the parameters of the Chaplygin models one needs both new observations and additional diagnostic techniques. From this point of view an application of the statefinder analysis could be useful.

Similar conclusions also hold for a one-fluid model of a generalized Chaplygin gas with a modified equation of state, as introduced in :

$$p=\frac{A}{\rho ^\alpha },$$
(21)
with $`0\alpha 1`$. This gives a cosmological evolution from an initial dust-like behavior to an asymptotic cosmological constant, with an intermediate epoch that can be seen as a mixture of a cosmological constant with a fluid obeying the state equation $`p=\alpha \rho `$ ($`\alpha =0`$ corresponds to the $`\mathrm{\Lambda }`$CDM model). This generalized model was studied in some detail in . Equation (21) gives the following dependence of the energy density $`\rho `$ on the scale factor $`a`$:

$$\rho =\left(A+\frac{B}{a^{3(\alpha +1)}}\right)^{\frac{1}{\alpha +1}}.$$
(22)
In this case the squared velocity of sound

$$v_s^2=\frac{p}{\rho }=\frac{\alpha p}{\rho }=\frac{A\alpha }{A+\frac{B}{a^{3(\alpha +1)}}}$$
(23)
varies from $`0`$ to $`\alpha `$. From Eq. (9) it follows that $`s=v_s^2\alpha `$ and therefore the admissible values of $`s`$ are now in the interval between $`\alpha `$ and $`0`$. From Eq. (9) we get

$$r=1\frac{9}{2}\frac{s(s+\alpha )}{\alpha }$$
(24)
For small values of $`\alpha `$ the generalized Chaplygin gas model becomes indistinguishable from the standard $`\mathrm{\Lambda }`$CDM cosmological model.

In conclusion. From the statefinder viewpoint both the pure Chaplygin model and the two-fluid mixture have different behaviors w.r.t. other commonly studied models. In particular, a future larger amount of data on high $`z`$ type Ia supernovae may allow to distinguish between the pure Chaplygin gas and the $`\mathrm{\Lambda }`$CDM models. Instead, the mixed two-fluid model becomes practically indistinguishable from $`\mathrm{\Lambda }`$CDM for sufficiently large values of the parameter $`\kappa `$ ($`\kappa >\kappa _05`$). On the other hand a mixed two-fluid model with $`\kappa `$ of order one is attractive from the point of view of a possible solution of the cosmic coincidence problem.

A.K. is grateful to CARIPLO Science Foundation and to University of Insubria for financial support. His work was also partially supported by the Russian Foundation for Basic Research under the grants No 02-02-16817 and 00-15-96699. A.K. is grateful to A.A. Starobinsky for useful discussions.


text:
# Results on Photon Production in Au+Au Collisions at RHIC

## 1 Introduction

The aim of this analysis is to search for direct photons, i.e. photons that are not produced in hadron decays such as $`\pi ^0\gamma \gamma `$ or $`\eta \gamma \gamma `$. Direct photons are generally subdivided into prompt photons from initial hard parton scatterings and thermal photons from a possible quark-gluon plasma and from the hadron gas.

Prompt photons are produced before a possible quark-gluon plasma has formed. Above a transverse momentum of $`p_T35`$ GeV prompt photons are expected to be the dominant source of direct photons in Au+Au collisions at RHIC. In p+p reactions direct photons are used to obtain information about the gluon distribution in the proton since their production in p+p is dominated by quark-gluon Compton scattering ($`gq\gamma q`$). By measuring direct photons in Au+Au at RHIC one is sensitive to modifications of the parton distribution functions in the nucleus (nuclear shadowing) .

The agreement between perturbative QCD calculations and experimental direct photon data in p+p and p+A can be improved if a momentum component $`k_T`$ of the initial partons transverse to the beam axis is introduced in the calculations . In nucleus-nucleus collisions one expects a stronger $`k_T`$ than in p+p due to e.g. multiple soft scatterings of the incoming nucleons. As demonstrated in  the magnitude of this intrinsic $`k_T`$ can be constrained by measurement of direct photons in nuclear collisions.

Thermal photons are produced over almost the entire evolution of a nucleus-nucleus collision: in a possible quark-gluon plasma phase and in the later hadron gas phase. The strongest contribution from thermal photons is expected roughly in the range $`p_T15`$ GeV. Theoretically, the entire hydrodynamical evolution of a nucleus-nucleus reaction needs to be modeled to describe thermal photon production. A comparison of these models with data constrains the initial temperature of the reaction system. In order to obtain information about the existence of a quark-gluon plasma one needs to compare models with and without phase transition to the data .

The first observation of direct photons in heavy ion collisions was achieved at the CERN SPS by the WA98 experiment . At RHIC, a substantial suppression of $`\pi ^0`$ production at high $`p_T`$ in central reactions was found that was not present at SPS energies . The reduced background from neutral pions could result in a clearer direct photon signal at RHIC. Moreover, the passage of quark jets through the excited nuclear medium might be an important new source of direct photons at RHIC, such that the subdivision of direct photons into prompt and thermal photons would have to be extended .

## 2 Data Analysis

Direct photon measurements (even in p+p) are generally considered to be very challenging. It is therefore an advantage of the PHENIX experiment that photons can be measured in different ways: directly with an electromagnetic calorimeter and also via their conversions in $`e^+e^{}`$-pairs. Here we focus on the photon measurement with the electromagnetic calorimeter (EMCal) of the PHENIX experiment. The biggest experimental uncertainties in photon measurements with calorimeters are generally background from hadrons that are misidentified as photons and errors of the energy scale of the detector. In heavy ion collisions an additional uncertainty results from the high multiplicity environment leading to overlapping showers in the detector.

In p+p reactions direct photons at sufficiently high $`p_T`$ are usually identified on an event-by-event basis with the help of certain isolation cuts. This is more difficult in heavy ion reactions and therefore the approach in this analysis is to try to find direct photons on a statistical basis. The inclusive photons (hadron decay photons and direct photons) are measured and then the background of photons from hadron decays is subtracted.

The PHENIX EMCal consists of two sub-detectors, a lead-scintillator calorimeter (PbSc, 6 sectors) and a lead-glass calorimeter (PbGl, 2 sectors). The lead-glass detector was previously used in the WA98 experiment. Each sector covers a pseudorapidity range of $`|\eta |<0.35`$ and an azimuthal angle of $`\varphi 22.5^{}`$. Both sub-detectors are highly segmented ($`\mathrm{\Delta }\varphi \times \mathrm{\Delta }\eta 0.01\times 0.01`$) such that the two decay photons of a $`\pi ^0`$ are well separated up to neutral pion momenta of $`p_T20`$ GeV. The different detection mechanisms of the two sub-detectors (measurement of scintillation light in PbSc and detection of Cherenkov photons in PbGl) result in a different response to hadrons. Thus, PbSc and PbGl provide photon measurements with different systematics.

The extraction of the inclusive photon spectra starts by identifying photon-like hits in the EMCal. Hadrons on the average produce showers with larger lateral extensions than photons. Thus, cuts on the shower shape are used to identify photons both in PbSc and PbGl. In addition, the time-of-flight of the showers is used as an identification criterion. The next analysis step is the subtraction of background from charged hadrons, neutrons, anti-neutrons, and particles not coming from the vertex. This is done statistically and not on an event-by-event basis. In order to determine the background of charged hadrons in the sample of photon-like hits, identified tracks from the PHENIX tracking system are projected to the EMCal surface. Random associations between photons and charged tracks are corrected for by an event-mixing technique. Above photon transverse momenta of $`p_T1.5`$ GeV the background of charged hits in the photon-like sample is found to be less than 10%. The background of neutrons and anti-neutrons in the sample of neutral EMCal hits, that is obtained after subtracting the charged background, has to be determined from simulation. The assumed input spectra of neutrons and anti-neutrons needed in the simulation are estimated based on measured $`p`$ and $`\overline{p}`$ spectra. The contribution of neutrons and anti-neutrons to all neutral EMCal hits is found to be less than 6%. In the last analysis step a photon efficiency correction is applied which takes the energy resolution of the detector, shower overlaps and photon losses due to identification cuts into account.

The expected photons from hadron decays are determined in a simulation that takes the measured neutral pion spectrum as input. Exactly the same data set is used for the neutral pion and the photon analysis. The neutral pion spectra obtained in this analysis agree with those presented in  within the current systematic error of 25%. An $`\eta `$-spectrum has not yet been measured. Therefore $`m_T`$-scaling is assumed for the $`\eta `$ and other hadrons. The result of the hadron decay simulation is shown in Figure 1.

## 3 Results

In order to search for non-hadronic photon sources the ratio of all measured photons to the expected background photons from hadron decays is shown as a function of $`p_T`$ for Au+Au collisions at $`\sqrt{s_{NN}}=130`$ GeV (Figure 2) and 200 GeV (Figure 3).

The results in Figure 2 are from an analysis that is similar to the one described in the previous section, but the charged background was determined from simulation and not from real data. The statistics in the 130 GeV data is small compared to the 200 GeV data such that the shape of the $`\pi ^0`$ spectrum is not as well constrained as in the 200 GeV analysis. In Figure 3 the results for Au+Au at $`\sqrt{s_{NN}}=200`$ GeV are presented in terms of the ratio of the measured $`\gamma /\pi ^0`$ ratio to the simulated $`\gamma /\pi ^0`$ ratio. In the $`\gamma /\pi ^0`$ ratio systematic errors partially cancel. The current systematic errors are dominated by uncertainties of the photon and neutral pion reconstruction efficiencies.

As can been seen from Figures 2 and 3 no direct photon signal is observed within current errors both in Au+Au collisions at $`\sqrt{s_{NN}}=130`$ GeV and $`\sqrt{s_{NN}}=200`$ GeV. The systematic uncertainties are expected to be reduced significantly in the future with further analysis.


text:
# The mass-function of primordial star clusters

## 1 Introduction

In a Cold Dark Matter (CDM) cosmology, small structures are the first to collapse and these then cluster together in a hierarchical fashion, giving rise to the bottom-up picture of galaxy formation. In this paper, we use the term ‘primordial star cluster’ for the first objects that are able to cool and form, zero-metallicity, Population iii stars. They are of interest both in their own right and because they may be responsible for reionization of the intergalactic medium.

In a primordial gas, whose main elements are hydrogen and helium and their derivatives, there are two main cooling mechanisms, dependent on the temperature: for halos with virial temperatures less than 8 600 K the cooling is dominated by roto-vibrational excitations of hydrogen molecules, while those with higher temperatures cool mainly via electronic transitions.

In a landmark paper entitled “How small were the first cosmological objects?”, Tegmark et al. (1997, hereafter T97) analytically tracked a top-hat collapse to the point of virialization, at which point the gas was cooled at constant density. They accepted an object as having cooled if it met the criterion $`T(0.75z_{\mathrm{vir}})0.75T_{\mathrm{vir}}`$, where $`T_{\mathrm{vir}}`$ is the virial temperature and $`z_{\mathrm{vir}}`$ the virialization redshift. They found that the first generation of objects that cooled in a standard CDM scenario virialized at a redshift of 27 and had a baryonic mass of about $`10^5`$$`\mathrm{M}_{}`$. In a later paper, Abel et al. (1998) redid the calculation with a different H2 cooling function and estimated a very similar virialization redshift but a smaller baryonic mass, $`7\times 10^3`$$`\mathrm{M}_{}`$.

In a paper somewhat cheekily entitled “How big were the first cosmological objects?”, (Hutchings et al., 2002, hereafter HSTC02) extended the previous work to include halos of higher mass. They found two distinct generations of halos: Generation 1 halos dominated by molecular cooling as in the previous work, and higher-mass, Generation 2 halos dominated by electronic cooling. The properties of these halos in three different versions of the CDM cosmology are listed in Table 1.

Note that the difference in cooling redshift between the two generations of halos is small (partly because the CDM fluctuation spectrum is flat on small scales and partly because cooling is more efficient in Generation 2 halos). This led HSTC02 to speculate that both generations of halos may form Population iii stars (whereas others have considered only the smaller, Generation 1 halos to be important).

There were two main deficiencies in the model of HSTC02. Firstly, they neglected substructure: the referee suggested that all Generation 2 halos will have Generation 1 halos inside them and so they will not be of primordial composition. We show below that primordial Generation 2 halos can exist. Secondly, they considered only 3-$`\sigma `$ fluctuations. In reality there will be a gaussian distribution of overdensities leading to a wide range of halo masses virialising at any given redshift.

The present paper, as a continuation of the previous work, addresses these two points by using the Block Model of Cole & Kaiser (1988) to generate a merger history of collapsed halos. This allows us to follow a wide dynamic range of halo masses very efficiently. Within the merger tree, the properties of halos are calculated using the same chemical model as in HSTC02.

Our work complements that of other authors who are investigating first object formation using numerical simulations, e.g. Abel et al. (1998), Bromm, Coppi & Larson (1999), Abel, Bryan & Norman (2000), Fuller & Couchman (2000), Nakamura & Umemura (2001), Bromm, Coppi & Larson (2002) and Nakamura & Umemura (2002). They are able to follow the dynamical evolution of single star clusters in great detail, whereas we learn instead about the properties of the cluster population.

We describe our numerical method in Section 2. The properties of individual halos for a merger tree corresponding to an overdense region of the Universe are presented in Section 3 and the corresponding mass function is described in Section 4. Finally, Section 5 explores variations of the basic model and discusses the nature of the star clusters.

## 2 Methodology

In this paper, we use the popular $`\mathrm{\Lambda }`$CDM cosmology whose parameters are tabulated in Section 2.1. The generation of a merger tree of collapsed halos is described in Section 2.2 and the criteria whereby we determine which of these form star clusters is outlined in Section 2.3.

### 2.1 Cosmology

HSTC02 investigated the cooling of halos in three different CDM cosmologies. In this paper, we restrict our attention to the currently-favoured $`\mathrm{\Lambda }`$CDM cosmology whose parameters are listed in Table 2.

We have used the transfer function calculated by CMBFAST. There have been recent suggestions that the normalisation of the power spectrum may be closer to 0.7 than 0.9 (e.g. Seljak 2002; Allen et al. 2002) this would have the effect of moving the formation epoch of the first star clusters to lower redshift and also lowering the amount of substructure.

### 2.2 Block Model

We generate a halo merger tree using the Block Model of Cole & Kaiser (1988). This starts with a ‘root’ block of mass $`M_0`$ and density fluctuation $`\delta _0`$. In this paper we fix $`M_0=10^{11}`$$`\mathrm{M}_{}`$ and choose two different values of $`\delta _0`$ corresponding to a 3-$`\sigma `$ fluctuation ($`\delta _0=10.98`$) and the mean density ($`\delta _0=0`$).

Geometrically the block can be visualised as a cuboid with sides in the ratio 1:$`2^{1/3}`$:$`2^{2/3}`$, but the density fluctuations are calculated as for a spherical top-hat model of the same mass; given the uncertainty in the normalisation of the power spectrum, the distinction is of no importance. The block can be bisected by a plane perpendicular to its longest axis, creating two similar blocks of half the mass, $`M_1=M_0/2`$. To generate density fluctuations in these daughter blocks, we add power drawn at random from a gaussian distribution with variance $`\sigma ^2(M_1)\sigma ^2(M_0)`$; once again, this is an approximation, but a good one. A positive fluctuation is added to one block and an equal negative fluctuation to the other, so as to conserve the overall level of fluctuations in the root block.

The same procedure is then repeated, with each parent block being divided into two equal-mass daughters until the desired resolution is reached. In this paper, we use 21 levels, creating a total of $`2^{21}12.1\times 10^6`$ blocks with a minimum block mass of $`9.5\times 10^4`$$`\mathrm{M}_{}`$.

We use a simple model in which the collapse of blocks to form bound objects is determined only by their overdensity. To be precise, we assume them to virialise once their linear overdensity reaches $`\delta _c=1.69`$ (see Eke et al. 1996 and Lokas & Hoffman 2001 where only a very weak dependence on $`\mathrm{\Omega }`$ has been found).

The equation for the rate of growth of $`\delta `$ with redshift, $`z`$, is

$$\delta (z)=\frac{\delta (0)}{1+z}\frac{g(\mathrm{\Omega })}{g(\mathrm{\Omega }_0)},$$
(1)
where

$$g(\mathrm{\Omega })=2.5\frac{\mathrm{\Omega }}{\left(\frac{1}{70}+\frac{209}{140}\mathrm{\Omega }\frac{1}{140}\mathrm{\Omega }^2+\mathrm{\Omega }^{\frac{4}{7}}\right)}$$
(2)
is the growth suppression factor (Viana & Liddle, 1996) and

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_0}{\left(\mathrm{\Omega }_0+\left(\frac{1}{1+z}\right)^3(1\mathrm{\Omega }_0)\right)}.$$
(3)
In the binary tree generated by the Block Model, most of the blocks are contained in larger blocks of greater over-density. Often this will be the immediate parent (one of the two daughters of each parent will have lesser density) but it could also be a block further up the tree. Under these circumstances, the larger block will collapse before the smaller, and so the latter will never attain an independent existence as a virialised structure. We eliminate these under-dense blocks from the tree and call the remaining blocks ‘halos’.

We have performed one hundred realizations of the Block Model with different number seeds. We use these for both values of $`\delta _0`$ as the set of halos is the same in each case. Unless stated otherwise, the results presented below are averages over all realizations.

### 2.3 Halo evolution and the formation of star clusters

We wish to know whether a halo can form stars before it gets incorporated into some larger structure. To do this, we construct an artificial model in which the halo has no substructure and cools at constant density. The actual structure of halos will be highly complex but our model gives a reasonable estimate of what is going on, short of doing a prohibitively time-consuming simulation.

We begin with the smallest halos and work our way up the merger tree. Each halo is treated as an isolated, isothermal sphere, as described in Section 3.3 of HSTC02. The mean baryon density within the virial radius, $`r_{\mathrm{vir}}`$, is taken to be equal to

$$\rho _{\mathrm{vir}}=\left(\frac{\mathrm{\Delta }_c}{\mathrm{\Omega }}\right)\rho _{b0}(1+z_{\mathrm{vir}})^3$$
(4)
where $`z_{\mathrm{vir}}`$ is the virialization redshift, $`\rho _{b0}`$ is the current mean density of baryons in the Universe, and $`\mathrm{\Delta }_c`$ is the mean overdensity of the virialized halo in units of the critical density at that time, which we take to be $`\mathrm{\Delta }_c18\pi ^2\mathrm{\Omega }^{0.45}`$ (Eke et al., 1998).

We define the dynamical time for each halo to be

$$t_{\mathrm{dyn}}=\frac{1}{4\sqrt{2}}t_{\mathrm{vir}}$$
(5)
where $`t_{\mathrm{vir}}`$ is the age of the universe at the time of virialization.

The virial temperature of the halo is

$`T_{\mathrm{vir}}`$ $`=`$ $`{\displaystyle \frac{\mu m_H}{k_B}}{\displaystyle \frac{GM_{\mathrm{tot}}}{2r_{\mathrm{vir}}}}`$ (6)
$``$ $`40.8{\displaystyle \frac{\mu }{1.225}}(1+z_{\mathrm{vir}})\left({\displaystyle \frac{\mathrm{\Delta }_ch^2\mathrm{\Omega }_0}{18\pi ^2\mathrm{\Omega }}}\right)^{\frac{1}{3}}\left({\displaystyle \frac{M_{\mathrm{tot}}}{10^5\mathrm{M}_{}}}\right)^{\frac{2}{3}}\mathrm{K}.`$
Here $`M_{\mathrm{tot}}`$ is the total mass (dark plus baryonic), $`m_H`$ is the mass of a hydrogen atom, $`k_B`$ is the Boltzmann constant, G is the gravitational constant, and $`\mu `$ is the mean mass of particles in units of $`m_H`$. Note that equation (6) differs slightly from the equivalent expression in HSTC02 as the latter contains a typographical error.

The initial fractional abundance of molecular Hydrogen is taken to be $`1.1\times 10^6`$ as calculated by Galli & Palla (1998), and the initial ionization fraction is taken to be the maximum of the equilibrium value at $`T_{\mathrm{vir}}`$ and the residual value from the early Universe, $`1.33\times 10^4`$.

Starting from these initial conditions, we determine the time that it would take the gas to cool isochorically to $`T_{0.75}=0.75T_{\mathrm{vir}}`$, using the minimal model presented in Section 2 of HSTC02. This includes molecular hydrogen cooling, collisional excitation and ionization of hydrogen and helium, and inverse Compton cooling from cosmic microwave background photons.

If halos are unable to cool to $`T_{0.75}`$ before being swallowed by a more massive halo, then we assume that they are heated to the new virial temperature and that any substructure (which would be minimal anyway because of the long cooling time) is erased. Contrarily, halos that can cool to $`T_{0.75}`$ are assumed to be able to (instantaneously) cool further to low temperatures and to form a star cluster.

For halos whose cooling times are much shorter than their dynamical times, then the gas will likely never get heated to the virial temperature in the first place and the assumption of instantaneous star formation will be a good one. However, for halos in which the cooling time exceeds the dynamical time then it seems likely that our model will underestimate the time taken to form stars. We consider in Section 5.2 the effect of adding a time-delay before star-formation and show that it favours the second generation of halos.

Star clusters may or may not survive subsequent halo mergers but either way they are assumed to instantly contaminate their surroundings with metals. Thus the primordial star clusters are those that contain no smaller star clusters within them. We assume that the metals do not propagate into halos on other branches of the tree. Thus metals may be ejected from star clusters but are confined within the next level of the merger hierarchy. This can be justified by a self-regulated model of star formation in which star-formation is terminated once gas is expelled from the star cluster.

We also assume that there is no external radiation field, other than that provided by the cosmic microwave background. Primordial star clusters will be surrounded by neutral gas and the propagation of ionizing photons will be severely limited. Nevertheless, these first objects will be highly clustered and so they will at some stage begin to interact with each other. We hope to consider this in a future paper.

## 3 Halo properties

In this section, we present results for the overdense, 3-$`\sigma `$, root halo. A comparison with the mean-density root halo will be presented in Section 5.3.

### 3.1 Collapsed halos

We start by considering the properties of all collapsed halos, i.e. halos whose linear overdensity exceeds that of all the (more-massive) halos within which they are contained. In Figure 1 we plot the virialization redshift of such halos (drawn from all 100 realisations) against their virial temperature.

The banding comes about because the halos come in fixed masses. Thus the smallest halos of mass $`9.5\times 10^4\mathrm{M}_{}`$ correspond to the left-most band; they have a wide range of virialization redshifts and temperatures that vary between about 200 and 1200 K; the top-most point corresponds to a 6.2-$`\sigma `$ and the lowest one to a 0.9-$`\sigma `$ fluctuation. 20 other bands are then visible, one for each factor of two in mass until we reach the parent halo in the bottom-right which has a mass of $`10^{11}\mathrm{M}_{}`$, a virial temperature of $`T_{\mathrm{vir}}=8.81\times 10^5`$K, and a virialization redshift of $`z_{\mathrm{vir}}=7.1`$ (corresponding to a 3-$`\sigma `$ fluctuation on this scale). The bands are mostly parallel, except for the temperature range $`T_{\mathrm{vir}}`$10 000–20 000 K within which the ionization level is changing.

The top panel of Figure 2 shows the cooling time, $`t_{\mathrm{cool}}`$, versus the virial temperature for each of the halos shown in Figure 1. There is a sharp decline in the cooling time at $`T_{\mathrm{vir}}\mathrm{10\hspace{0.17em}000}`$ K corresponding to the ionization temperature of hydrogen. Halos with higher virial temperatures than this (Generation 2 halos) are able to cool rapidly via electronic processes and so have relatively short cooling times. Those with lower virial temperatures (Generation1) have to rely on cooling via molecular hydrogen which is formed only in very low quantities. Although our definition of $`t_{\mathrm{cool}}`$ only follows cooling down to 0.75 $`T_{\mathrm{vir}}`$, we note that Generation 2 halos have a high residual ionization at lower temperatures that acts as a catalyst to form molecular hydrogen: thus their cooling rates at low temperature are faster than for Generation 1 halos.

The lower panel of Figure 2 shows the dynamical time for the collapsed halos, as defined in equation 5. A a rule of thumb, halos with virial temperatures above about 10 000 K are able to cool in less than a dynamical time; cooler halos take longer.

### 3.2 Star clusters

The time difference, $`\mathrm{\Delta }t_{\mathrm{coll}}`$, between the collapse time of each halo and that of its parent is shown in the top panel of Figure 3. For a random location in space, one might expect that the time difference would be largest for massive, high-temperature halos. However, that is not the case for these realisations in which the top-level itself is constrained to collapse at $`t=0.69`$ Gyr (it is true for the halos considered in Section 5.3 for which the overdensity of the top-level halo is zero).

Those halos for which $`t_{\mathrm{cool}}`$ is less than $`\mathrm{\Delta }t_{\mathrm{coll}}`$ can cool to a fraction 0.75 of their virial temperature before being swallowed up by their parent halo in the merger hierarchy. To begin with, we assume that this is a sufficient criterion to allow them to form stars and we identify them with star clusters. In reality the time-delay before star formation will be larger as the gas has to cool to low temperatures and to congregate into regions of high-density. We will consider the effect of allowing a longer time-delay in Section 5.2.

The lower panel in Figure 3 plots $`\mathrm{\Delta }t_{\mathrm{coll}}t_{\mathrm{cool}}`$ against the virial temperature. Those halos that lie above the line are those that form star clusters. Just 2 percent of all collapsed halos satisfy this condition. However, these are not distributed evenly over mass. For example, in the first three levels of the merger tree (i.e. the three levels with the lowest mass) only a fraction $`1.2\times 10^7`$, $`8.4\times 10^5`$ and $`3.2\times 10^3`$ of collapsed halos are able to form stars, whereas a successively higher fraction do so at higher mass. Only for the most massive halos with virial temperatures in excess of $`10^5K`$ does the cooling time again begin to exceed the lifetime of halos.

### 3.3 Primordial star clusters

We assume that metal-enrichment from star-formation is instantaneous but that it does not extend beyond the immediate environment of a star cluster and its parent halo. Then star clusters of primordial composition are simply those which do not have any smaller star clusters contained within them.

Approximately half of these clusters ($`56.8\%`$) in our 3-$`\sigma `$ realisations satisfy this condition with the bias swinging back towards low masses. The most massive primordial star cluster has a mass of $`9.76\times 10^7\mathrm{M}_{}`$ and a virial temperature of $`1.54\times 10^4K`$.

## 4 Halo mass functions

The number of star clusters as a function of virial temperature, averaged over all 100 realisations, is shown in the upper panel of Figure 4. The upper histogram shows all star clusters, whereas the lower is restricted to primordial star clusters. There is a clear minimum at about 8 600 K corresponding to the division between Generation 1 halos on the left and Generation 2 halos on the right. Note that the star clusters that make up the upper histogram are not all independent. That is to say that many of the low-temperature clusters are subcomponents of the higher-temperature ones. However, the primordial star clusters are all distinct objects. The y-scale in panel (a) could be multiplied by 3.94$`h^3`$Mpc-3 to convert to a number density but we have not done this as the 3-$`\sigma `$ region that we consider is not representative of all space.

The lower panel of Figure 4 shows the same distribution but weighted by mass. From this it is apparent that, whereas the majority of the primordial star clusters in this region are Generation 1 objects, the division of mass between the two generations is much more even, with only about twice as much mass contained in Generation 1 as compared to Generation 2 objects.

Note the sharp cut-off in the mass density of primordial star clusters at virial temperatures greater than about 15 000 K. This is because all higher temperature clusters contain a Generation 2 subcluster for which the cooling time is very short and which can itself form stars on a short timescale. This situation changes when we introduce a time-delay for star formation in Section 5.2.

Figure 5 is similar to Figure 4 except that the ordinate is now virial mass rather than virial temperature (to convert to baryonic mass, the x-scale should be multiplied by 0.12). In panel (b), the dotted line shows the contribution to the total mass of Generation 1 halos (whose virial temperature is less than $`\mathrm{8\hspace{0.17em}600}`$K), while the dot-dashed line is for Generation 2 halos. The fractional mass contained in the two generations is 0.109 and 0.049, respectively. Thus about 16 per cent of all baryons in this 3-$`\sigma `$ region of space will have at one time been part of a primordial star cluster.

Figure 6 contrasts the collapse and star-formation redshifts of both generations of halo. The spiky features visible in the distributions are due to the factor of two mass resolution of our halos and would be smoothed out in a more general merger tree. Because of the long cooling times of Generation 1 halos, the difference in the peaks of the distributions of star-formation redshifts of the two generations is not so great as for their collapse redshifts. Nevertheless, it is clear from the figure that a significant fraction of Generation 1 halos both collapse and form stars before Generation 2 halos begin to form in numbers. This highlights the need for a more sophisticated model of feedback than we attempt in this paper.

## 5 Discussion

### 5.1 Numerical considerations

It is legitimate to ask to what extent our results are limited by the factor-of-two mass resolution inherent in the Block model. If we could have subhalos with a wider range of masses would that lead to a greater probability of contamination by star formation and a reduction in the fraction of primordial halos? The dramatic reduction in the cooling time of Generation 2 halos compared to their low-mass, Generation 1 subhalos, as illustrated by the upper panel in Figure 2, suggests that this is unlikely to be the case. We expect to move to a more realistic merger tree in future work.

Meanwhile, we have tested the sensitivity of our results to the precise choice of halos masses by performing a second hundred realizations of the merger tree with the root halo mass (and hence each level of the merger hierarchy) increased by a factor of $`\sqrt{2}`$, to $`1.4\times 10^{11}`$$`\mathrm{M}_{}`$. Figure 7 shows the mass function of primordial halos for these simulations contrasted with our original simulation (dotted line). There is no significant difference between the two.

The results that we have presented so far are an average over a large number of realisations. Figure 8 shows the dispersion around the average, of five of the 100 realizations carried out in this paper. It can be seen that the scatter is significant but not enough to seriously affect the divide between the two generations of star clusters within each realisation.

### 5.2 Time-delayed star-formation

So far we have assumed that after the cooling of the gas to low temperatures (following the T97 criterion) stars form instantaneously. In reality, there will be a lapse of time until the gas reaches the high-density regime in which nuclear reactions can take place and the stars are born. In an attempt to include in our code a time-delay between initial cooling and star-formation, we consider in this section the effect of adding the dynamical time to the cooling time. The justification for this is simply that, following virialization, one would expect the gas to take at least a dynamical time to contract within the potential well of the halo (this argument has less force for high-mass halos whose cooling time is very much shorter than their dynamical time and which may therefore never attain virial equilibrium in the first place). Some evidence for this delayed star formation comes from Abel et al. (2002) who describe the formation of a primordial star using 3-D hydrodynamical simulation. Their results show that a cooled ($`200K`$) high redshift molecular cloud is formed at $`z=24`$ and then a proto-star is formed at $`z=18.2`$. The time that took to form this proto-star is of the order of the dynamical time of the cloud.

Panel (a) of Figure 9 shows a histogram of the fractional mass of primordial star clusters as a function of their virial temperatures. The dotted line corresponds to the original case in which no time delay has been added, while the full line shows the distribution when a dynamical time has been added to their cooling time. The fractional mass has changed in such a way that now we have 3 times more mass in Generation 2 halos than in Generation 1 halos, and a greater total mass fraction than before. Note also that there is no longer a sharp cut-off at virial temperatures above 15 000 K because it is possible for subhalos to have short cooling times in this model and yet not to form stars.

Panel (a) in Figure 10 shows the equivalent mass function of these clusters from which it can be seen that halos as massive as 1010$`\mathrm{M}_{}`$ can contain primordial star clusters. While this does not seen very likely, the general conclusion that we draw is that delayed, rather than instantaneous, star-formation will favour Generation 2 halos over Generation 1.

We also tried a model with instantaneous star-formation but with a time-delay before energy and metallicity feedback. The idea is that if the time difference between the collapse of a parent halo and the cooling of its child is less than the time for the formation of supernovae, $`10^7`$yr, then the parent halo will be of primordial composition and has to be added into the set of primordial objects. However, this makes only a minor difference to our results and we will not discuss it further.

### 5.3 Mean-density regions

The first star clusters will form in overdense regions of the Universe, hence our use of a 3-$`\sigma `$ root halo to this point. However, it is interesting to contrast these results with those expected for a more typical part of the Universe, with density equal to the cosmic mean.

Figure 11 contrasts the collapse redshifts of primordial star clusters for the 3-$`\sigma `$ and mean density regions.

It can be seen that the halos collapse at much lower redshifts in the mean overdensity case. Because the mean-density regions are likely to be far-removed from the regions of the first star-formation, they are unlikely to be affected by photoionizing photons at high redshift. However observations (e.g. Fan et al., 2002) and simulations (e.g. Razoumov et al., 2002) both suggest that the Universe became re-ionized at a redshift of about 6 and our model will be invalid after this time. This will mostly affect the evolution of Generation 2 halos.

The effect of the lower formation redshift on the virial temperatures and masses of the star clusters is shown in the lower panels of Figures 9 and 10. The fractional mass distribution over virial temperature is almost unchanged, but with a slight bias towards lower temperatures compared with the 3-$`\sigma `$ case. A greater effect is a shift in the mass function towards higher masses: because the halos collapse at lower redshift and hence have lower densities, they have higher masses for a given virial temperature.

### 5.4 What and where are they now?

Our model predicts the masses, virial temperatures and formation redshifts of primordial star clusters, but says nothing about their internal structure. Hydrodynamical simulations (see references in the Introduction) have made a start in this direction but are as yet still in their infancy. There has been some theoretical speculation about the masses of the first stars but no consensus has emerged. In this section, we use our results to discuss two possible fates of primordial star clusters, but note that the physics is sufficiently uncertain that we may even have got the roles of the two generations mixed up.

The baryonic mass of our root halos, $`1.2\times 10^{10}`$$`\mathrm{M}_{}`$, is similar to that of a normal galaxy of mass approximately one-tenth that of an $`L_{}`$ galaxy. The space-density for the 3-sigma fluctuations on this scale is $`3.0\times 10^3h^3`$Mpc-3, similar to that of groups of galaxies, so that we would perhaps expect one such galaxy in a typical group. The other galaxies will form slightly later and so the star clusters will be biased to higher masses, although the total mass contained in primordial star clusters will be similar (see Section 5.3). Our model therefore suggests that approximately one tenth of the baryons in a typical galaxy will have passed through a primordial star cluster. The majority of these are probably enriched with processed material without themselves forming a zero-metallicity star.

The majority of Population iii stars will be born in regions that are destined to end up in normal galaxies. However, our model does not preclude the formation of some zero-metallicity stars in low-density regions of the Universe at relatively low redshift (but before re-ionization). The resulting star clusters would be of low-density and therefore very hard to detect.

The bulk of Generation 1 star clusters at high redshift have masses in the range $`10^6`$$`10^7`$$`\mathrm{M}_{}`$; the baryonic mass is lower, $`10^5`$$`10^6`$$`\mathrm{M}_{}`$. It is natural, therefore, to identify these objects with the low-metallicity globular clusters found in the bulges and halos of normal galaxies. The relatively long cooling times of Generation 1 halos compared to their dynamical times would have aided dissipative collapse within the dark matter halo and survivability of the star cluster. One objection to this is that zero-metallicity stars have not been discovered in globulars, but of course the first stars may have been of high mass and may have burnt out long ago. A more serious objection is that we know that the amount of material in globular clusters is much less than one tenth of all the baryons in a galaxy. It is therefore probable that feedback of energy from the first supernovae (or hypernovae) will disrupt the star clusters and that the majority of zero-metallicity stars, should they still exist, will be spread throughout the bulges of normal galaxies.

The main differences with Generation 2 clusters is that they are more massive by about a factor of ten and that their cooling times are much shorter. It is interesting to speculate that electronic cooling when it did turn on would lead to catastrophic accumulation of cold gas at the centre of the collapsing halo and perhaps to the formation of a massive black hole. Observations (e.g. Ferrarese & Merritt 2000; Gebhardt et al. 2000) give a galactic black-hole to bulge mass ratio of about 0.001. To be consistent with this would require an accretion efficiency of just 2 per cent, creating seed holes of mass $`2\times 10^4`$$`2\times 10^5`$$`\mathrm{M}_{}`$. Subsequent merging of these seed holes could lead to the formation of supermassive black holes in the centres of normal galaxies today.

## Acknowledgements

FS thanks his parents for their support; PAT is a PPARC Lecturer Fellow.


text:
# Simulated X-ray Cluster Temperature Maps

## 1 INTRODUCTION

Observations of X-ray clusters of galaxies over the last few years have shown that clusters are not the smooth, featureless systems they were expected to be. X-ray surface brightness and temperature observations indicate the presence of substructure and support the view that cluster formation occurs through the infall and merger of subclusters.

Merger shocks would be expected to occur in such a scenario (e.g. Markevitch & Vikhlinin, 2001; Markevitch et al., 2002) but additionally Chandra observations have revealed a new phenomenon of ‘cold fronts’ (Markevitch et al, 2000; Forman et al., 2001; Mazzotta, Fusco-Femiano & Vikhlinin, 2002). At a cold front, the entropy jump across the sharp gas density discontinuity is in the opposite sense to that expected for a shock, with the high surface brightness side of the dense edge corresponding to low temperature.

The variety of features that have been observed in the clusters mapped so far are of great interest since they may shed light on the physical processes which are occurring as well as containing information on the stage, geometry, scale and velocity of the mergers. A number of authors have carried out simulations of controlled, single mergers (see Ritchie & Thomas 2002 for a review). For example, Ritchie & Thomas (2002) carried out high-resolution simulations of the merger of idealized clusters containing both dark matter and gas. They studied the effect on the observable properties of clusters of single head-on and off-centre mergers between both equal and unequal mass objects. A sequence of maps of emission-weighted temperature with superimposed X-ray surface brightness and velocity fields for mergers between equal mass systems showed the compression and shocking of the gas as the merger progressed. Recently Nagai & Kravtsov (2002) carried out a detailed study of cold fronts in high resolution simulations of two clusters forming in different Cold Dark Matter models (standard CDM and $`\mathrm{\Lambda }`$CDM). Their results indicate that cold fronts are probably fairly common but are non-equilibrium transient phenomena.

In this paper, we present preliminary results from temperature maps of an ensemble of clusters that form within a cosmological simulation, already shown to reproduce the observed X-ray scaling relations at low redshift (Muanwong et al. 2002). This way, we are able to directly assess the range of substructure present in the cluster population. In common with observations, we see significant temperature fluctuations in the hot gas, even when there is little information present in the surface brightness distribution. In particular, we discuss the presence of a cold front in one of our bimodal clusters, which has properties consistent with the cold front in Abell 2142 (Markevitch et al. 2000).

## 2 The simulations

Simulation data were generated using a parallel version of the hydra $`N`$-body/hydrodynamics code (Couchman, Thomas & Pearce 1995; Pearce & Couchman 1997) on the Cray T3E supercomputer at the Edinburgh Parallel Computing Centre. Full details of the simulations have already been presented elsewhere (Thomas et al. 2002; Muanwong et al. 2002, hereafter M2002), so we summarize pertinent details only.

Results are presented for a flat, low-density cosmology ($`\mathrm{\Omega }_0=0.35,\mathrm{\Omega }_\mathrm{\Lambda }=0.65,h=0.71,\sigma _8=0.9`$), using $`160^3`$ each of gas and dark matter particles within a box of comoving length, $`100h^1`$Mpc. We consider 2 models: the Radiative model includes radiative cooling of the gas, and the Preheating model for which we additionally preheat the gas by raising its specific thermal energy by 1.5 keV per particle at $`z=4`$. In both models, we adopted a time-dependent global metallicity, $`Z=0.3(t/t_0)Z_{}`$, where $`t_013`$Gyr is the current age of the universe. As discussed in M2002, both models reproduce the X-ray cluster scaling relations, although the Radiative model contains a significantly higher cooled fraction than the Preheating model (15 per cent, as opposed to $`0.5`$ per cent). In the former model, cooling is limited by numerical resolution. The cooled gas fraction is larger than that determined by Balogh et al. (2001), using results from the 2MASS and 2dF galaxy surveys (Cole et al. 2001), but in close agreement with the observed value from the SDSS (Blanton et al. 2001).

Simulated cluster catalogues were produced using the same method as described in M2002. In summary, clusters were identified as clumps of particles within spheres of average overdensities compared to the comoving critical density, centred on the position of the densest dark matter particle. For this paper, we adopt an overdensity of 1000; the associated radius, $`R_{1000}`$ is approximately half the size of the virial radius, and is comparable to the extent currently probed by Chandra.

## 3 The temperature and surface brightness maps

Maps of emission-weighted temperature were made for each of the 9 most massive clusters ($`M_{1000}=14\times 10^{14}h^1\mathrm{M}_{}`$) in the Preheating catalogue, and for the corresponding clusters in the Radiative catalogue. Each map was produced by first locating all hot ($`T>10^5`$K) gas particles within a cube of length, $`2R_{1000}`$ along each side, centred on the cluster of interest. The desired quantity was then smoothed onto a 3D array using the SPH interpolation method (e.g. Monaghan 1992)

$$A(𝐫_j)=\frac{\underset{i}{}A_iw_i𝒲(\mathrm{\Delta }r_{ij},h_i)}{_iw_i𝒲(\mathrm{\Delta }r_{ij},h_i)},$$
(1)
where

$$𝒲(\mathrm{\Delta }r_{ij},h_i)=\frac{W(\mathrm{\Delta }r_{ij},h_i)}{_kW(\mathrm{\Delta }r_{ik},h_i)}.$$
(2)
The sums $`i`$ extend over all particles and $`k`$ over all voxels. $`A_i`$ is the value of the quantity for particle $`i`$ at position $`𝐫_i`$, $`𝐫_j`$ is the centroid position of voxel $`j`$, $`\mathrm{\Delta }r_{ij}=|𝐫_j𝐫_i|`$, $`w_i`$ is a weight factor, $`W(\mathrm{\Delta }r,h)`$ is the same SPH smoothing kernel used by hydra (Thomas & Couchman 1992) and $`h_i`$ is the smoothing length of particle $`i`$ ($`W0`$ as $`\mathrm{\Delta }r2h`$). Equation 2 normalizes the contribution of each particle to be $`A_i`$ when summed over all voxels. For all mass-weighted quantities, $`w_i=m_i`$ and for emission-weighted quantities, $`w_i=m_in_i\mathrm{\Lambda }(T_i,Z)`$, for particle $`i`$ with mass $`m_i`$, density $`n_i`$, temperature $`T_i`$ and X-ray emissivity, $`ϵ_i=n_i^2\mathrm{\Lambda }(T_i,Z)`$; the cooling function, $`\mathrm{\Lambda }(T,Z)`$ is the same function adopted for our simulations, using tables published in Sutherland & Doptia (1993). The smoothed distribution was then projected onto a $`64\times 64`$ grid, with the width of each pixel being approximately 7 to 10 $`h^1\mathrm{kpc}`$, smaller than the gravitational softening length used in the simulations (25 $`h^1\mathrm{kpc}`$).

Figs. 1 & 2 show the emission-weighted temperature maps (with superimposed X-ray surface brightness contours and projected velocity vectors) for the Preheating and Radiative runs respectively. Surface brightness contours are normalized to the maximum value, with each contour representing variation by an order of magnitude. Velocities (in km/s) are mass-weighted and the length of each velocity vector represents the tangential speed (i.e. across the map plane) of the gas at that point.

Positive temperature fluctuations are evident of up to 3 times the average, associated with compression of the gas. There also exists bright clumps of cool gas, with temperatures less than half the mean value. These features are far more abundant in the Radiative maps than in the Preheating maps because the energy injection in the Preheating run was large enough to erase a significant amount of substructure.

Among the 9 clusters shown in Fig. 1, several show clear bimodal structure (clusters 1, 2 and 8). This is most clear for the 2nd and 8th maps, in each of which 2 cold clumps appear to be moving away from one another. In the Radiative maps (Fig. 2) the bimodal structure of map 1 is not so clear. The 8th cluster shows evidence of a sharp boundary at the leading edge of the cold clump in the upper right quadrant; when the maps for this cluster are re-centred to show the cold clump in the lower right of the map, a similar sharp edge can be seen. These features are indicative of observed cold fronts and this is investigated in the next section.

The velocity vectors superimposed on the cluster maps show clear evidence for rotation in some cases (for example maps 3, 5 and 9). Ritchie & Thomas (2001) found similar rotation of the velocity field in their simulated off-centre merger of two equal mass clusters. Except for map 9. where the surface brightness contours are essentially round, the contours are generally elliptical, again as found for merging clusters. The spiral patterns seen in some of the temperature maps (e.g. 3 & 5) may be due to heating by shocks produced during the merger as was seen in the simulations of Ritchie & Thomas during the interaction of cluster cores. With its rounder surface brightness contours, cluster 9 may be in a later, more relaxed state.

The centres of all clusters show a decline in temperature. This is in contrast to the simulated temperature maps of Loken et al. (2002), produced using an Adaptive Mesh-Refinement code, where the radial temperature profiles continue to rise at small radii. The Santa Barbara cluster comparison project (Frenk et al., 1999) found the general result that for non-radiative gas, SPH codes produce a central flat or slightly declining temperature profile, while all of the grid codes produce a temperature profile which continues to rise to the resolution limit. This conflict between results using the different codes still needs to be resolved. However, Fabian (2002) noted that all of the cores mapped so far by Chandra with radiative cooling times of a few Gyr show significant central temperature drops. This temperature drop of a factor of about 3 or more in the central region of the clusters was explained by Fabian as possibly due to a combination of radiative cooling and gas introduced from dense cooling subclusters, as we see here.

We present circularly-averaged emission-weighted temperature profiles for the Radiative maps in Fig. 3. The more regular of our clusters (4, 7, 9) exhibit the type of gradually declining temperature profile found by Loken et al. (2002). The others exhibit a variety of profiles with the irregular clusters (2, 6) even appearing approximately ishothermal when circularly averaged. It therefore seems doubtful that a universal temperature profile can be applied to galaxy clusters, at least in the inner regions (out to about 500 $`h^1`$ kpc, approximately half the virial radius).

## 4 Identification of a possible cold front in map 8

Of the nine cluster maps produced here, map 8 appeared to have the strongest visual evidence for the presence of cold fronts. To investigate whether the physical properties are consistent with a cold front, profiles were produced of emission-weighted temperature, X-ray surface brightness, pressure and density across the surface brightness and temperature discontinuity in the upper right quadarant of the Radiative map. The pixels with position angles between 30 and 60 with respect to the image centre were binned as a function of radius for each of the four parameters. The profiles are shown in Fig. 4.

A sharp rise in temperature by a factor of about 3 over a radial distance of $`<0.2h^1`$ Mpc was found, corresponding to a sudden change in the slope of the surface brightness profile. Comparison with the observed temperature profile across the northwestern edge of the X-ray image for Abell 2142 in Markevitch et al. (2000) shows that the radial scale is very similar and the temperature for both increases by approximately the same factor. In the region of the cold front, the pressure profile of the simulated cluster was smooth and fairly flat, falling off at larger radius. For the pressure profile to flatten in this region, where the temperature is rising sharply, implies that the gas density must be dropping sharply, again similar to the profile in Markevitch et al. (2000). The radial velocity of the gas in the cold clump is highly subsonic with respect to the gas to the right of the cold front.

Similar profiles, consistent with cold fronts, were found when the binning was repeated across the surface brightness edge in the lower left quadrant, and also for the Preheating run for the same cluster. The overall picture for this cluster, taking into account also the direction of the velocity vectors, is of two cold subclumps moving away from each other which fits in with the merger scenario suggested by Markevitch et al. (2000) of subclumps falling in along filaments. The cold cores are undisrupted by the merger and continue on through the cluster.

## 5 Discussion

Temperature maps of the 9 largest clusters in the simulated cluster catalogues of M2002 show clearly that the clusters do contain substructure. Features that are seen include cold clumps, some of which may be consistent with the infall of cold subclusters, hot spiral features consistent with shocks produced in mergers, and at least one map containing cold fronts. Elliptical surface brightness contours are the norm, and rotation as expected in off-centre mergers is common. The Preheating maps appear ’cleaner’, with fewer cold clumps present than for the Radiative maps, since the gas has been heated at z=4, wiping out many of the colder blobs. The remaining cold clumps seen in these maps did seem to be associated with motion into or through the cluster.

In this preliminary paper our main aim has been to demonstrate that simulated maps contain a variety of substructure. The scale of the structure is similar to that seen observationally and so gives some indication of the features which should be observable. Future work will focus on the evolution of structure in clusters and on a larger range of cluster masses extending down to the size of groups.

## 6 Acknowledgements

The simulations used in this paper were carried out on the Cray T3E at the EPCC as part of the Virgo Consortium programme of investigations into the formation of structure in the Universe. LIO was supported in part by the Leverhulme Trust and STK by PPARC; PAT is a PPARC Lecturer Fellow.


text:
# Integrability in 1D Quantum Chaos

## I Quantum Networks

As a model of a quantum chaotic system, let us consider a quantum particle that moves on a quasi one-dimensional network. On every bond $`B_{ij}`$ of the network one may define a potential $`U_{ij}`$, which turns the graph into a dressed quantum graph . Below we shall consider the case of constant scaled bond potentials $`U_{ij}=E\lambda _{ij}`$, where $`E`$ is the energy and $`0<\lambda _{ij}<1`$ are real scaling constants. Such potentials will change the action lengths of the bonds and the values of the reflection and transmission amplitudes of the quantum particle at the vertices of the graph. The wave function of the particle is given by

$$\mathrm{\Psi }_{ij}(x)=A_{ij}e^{ik_{ij}x_{ij}}+B_{ij}e^{ik_{ij}x_{ij}}.$$
(1)
The scaling assumption allows to avoid unnecessary complications, such as the geometrical restructuring of the periodic orbits for different energies.

The boundary conditions consist of continuity and flux conservation requirements for the wave function $`\mathrm{\Psi }_{ij}`$ across the vertices $`V_i`$ of the graph . Such boundary conditions produce a spectral determinant $`\mathrm{\Delta }(k)`$, whose zeroes define the spectrum of the problem. Usually, the quantum spectrum determines the energy levels for a given potential. However, for the scaling problems considered here, the scaling constants are given and the potential is determined by the spectrum.

The spectral determinant for the quantum networks is a finite exponential sum,

$`\mathrm{\Delta }(k)=1+e^{i2(S_0k\pi \gamma _0)}{\displaystyle \underset{i=1}{\overset{2N_\mathrm{\Gamma }}{}}}a_ie^{i2(S_ik\pi \gamma _i)},`$ (2)
where $`S_i`$ are constants that can be expressed via the action lengths of the graph bonds and the $`\gamma _i`$ are constants. The total length $`S_0`$ (the 1D volume) of the graph is the largest frequency that appears in the exponential sum (2), $`S_0>S_i`$. The number of terms $`N_\mathrm{\Gamma }`$ in the sum (2) is finite.

We are interested in studying the roots of $`\mathrm{\Delta }(k)`$. Such functions (for a more general case of a complex argument) were studied from the mathematical perspective in , where it was shown that the roots $`z_n`$ of (2) form an almost periodic set, which is uniformly bounded in the complex plane, i.e. is contained in a stripe of finite width,

$$\left|Imz_n\right|<M<\mathrm{},\mathrm{\Delta }(z_n)=0.$$
(3)
Moreover, it was shown that the real parts of the roots have the form

$$Rez_n=\frac{\pi }{S_0}n+\phi (n),$$
(4)
where $`\phi (n)`$ is an almost periodic function, whose frequency set contains linear combinations of the frequencies $`\pi S_i/S_0`$, with integer coefficients. In  these mathematical results were put into the perspective of quantum chaos theory that describes the spectral properties of quantum networks.

## II Quantum Chaos

The spectral problem of quantum networks has a long history (see  and the references therein). In fact, the problem of obtaining the energy levels in the step potential

$$V(x)=\{\begin{array}{cc}0,\hfill & \text{for }0<xb\text{,}\hfill \\ V_0=\lambda E,\hfill & \text{for }b<x<1\text{,}\hfill \\ \mathrm{},\hfill & \text{for }x<0\text{ or }x>L\text{,}\hfill \end{array}$$
(5)
is the one that is usually presented following the elementary “infinite square well” problem in many standard quantum mechanics textbooks . However, unlike in the square well problem, the spectral equation for the step potential (5) is not easy to solve analytically. In fact, the only known way to obtain the roots of the spectral equation are numerical or graphical solutions.

The complexity of the spectral equation for the system (5) as well as for a generic quantum network is in fact due to the dynamical complexity of the underlying classical system, which turns out to be chaotic .

In the classical limit the quantum graph systems reduce to a classical particle that scatters stochastically at the vertices of the graph (Fig. 3). The scattering probabilities are obtained in the $`\mathrm{}=0`$ limit from the quantum-mechanical vertex transition amplitudes $`t_{ji,ij}`$ . Every time the particle encounters a vertex $`V_i`$ traveling along the bond $`B_{ji}`$, it can reflect back from it with the probability $`|t_{ji,ij}|^2`$ or pass through to another bond $`B_{ik}`$ with the transmission probability $`|t_{ji,ik}|^2`$. After every reflection or transmission (scattering) event, the particle completely loses its memory about the previous stage of its motion.

Using the coefficients $`t_{ji,ik}`$, it is easy to obtain the probability amplitudes with which a particular periodic orbit is realized. Indeed, if a certain orbit $`\gamma `$ reflects $`\sigma _i(\gamma )`$ times off the vertex $`V_i`$ and passes through it $`\tau _i(\gamma )`$ times, then the probability amplitude for such an orbit is

$$A_\gamma =t_{ji,ij}^{\sigma _i(\gamma )}t_{ji,ik}^{\tau _i(\gamma )},$$
(6)
where the product is taken over all the vertices visited by $`\gamma `$. It is usually convenient to distinguish the so-called “primitive” orbits, i.e. the ones that traverse a certain sequence of bonds only once, and the “multiple traversal” orbits. If a certain primitive orbit $`\gamma _p`$ with the weight $`A_p`$ is traversed $`\nu `$ times, then the corresponding amplitude is $`A_p^\nu `$.

Due to the stochastic scattering, the trajectories of the particle are geometrically very complex and the number of possible orbits (and in particular the periodic orbits), increases exponentially . The number $`N`$ of the $`n`$-bond periodic orbits grows as

$$N\frac{e^{\tau n}}{n},$$
(7)
where the exponential proliferation rate $`\tau `$ depends only on the topology of the graph.

Such behavior mimics closely an important phase space feature of deterministically chaotic systems. Hence quantum networks provide very convenient and simple models for studying various features of quantum systems that correspond to classically chaotic systems in the context of quantum chaos theory.

## III Periodic orbit theory

Let us outline the main attributes of quantum chaos theory , applied to quantum networks. In  it was shown that for dressed quantum graphs there exists an exact periodic orbit expression for the density of states,

$`\rho (k)={\displaystyle \frac{S_0}{\pi }}+{\displaystyle \frac{1}{\pi }}Re{\displaystyle \underset{p}{}}S_p^0{\displaystyle \underset{\nu =1}{\overset{\mathrm{}}{}}}A_p^\nu e^{i\pi \nu S_p^0k},`$ (8)
which generalizes the earlier results obtained in . Here $`S_p^0`$ and $`A_p(E)`$ are correspondingly the action length and the weight factor (6) of the prime periodic orbit labeled by $`p`$, $`\nu `$ is the repetition index and $`S_0`$ is the total action length of the graph.

Another important characteristic of the spectrum is the spectral staircase function

$$N(k)=\underset{j}{}\mathrm{\Theta }(kk_j),$$
(9)
for which there also exists an exact periodic orbit expansion. In accordance with Weyl’s law, the average spectral staircase depends linearly on the momentum,

$$\overline{N}(k)=\frac{S_0k}{\pi }+\overline{N}(0),$$
(10)
where $`S_0`$ is the total length (1D volume) of the network. This statement implies that the average dependence of the momentum eigenvalues $`k_n`$ on their index  is

$$\overline{k}_n=\frac{\pi }{S_0}\left(n\overline{N}(0)\frac{1}{2}\right),$$
(11)
which complies with the relationship (4).

Weyl’s average (10) of the spectral staircase is the same for all networks with the same action length $`S_0`$. However, different networks vary significantly in how well the average approximates the actual staircase function $`N(k)`$. It is interesting that there exist certain graphs with a particularly regular behavior of their spectra. For these quantum networks (called regular in ), the average spectral staircase (10) pierces every stair-step of $`N(k)`$. It is clear that for the regular graphs the intersection points $`\widehat{k}_n`$ between the spectral staircase and its average,

$$N(\widehat{k}_n)=\overline{N}(\widehat{k}_n),$$
(12)
occur periodically on the momentum axis . In other words, for such systems there exists exactly one quantum eigenvalue $`k_n`$ between every two such intersections.

For such graphs one can obtain immediately the exact periodic orbit series representation for the individual eigenvalues of the momentum. Indeed, since every interval $`I_n=[\widehat{k}_{n1},\widehat{k}_n]`$ contains one delta-peak of $`\rho (k)`$, the corresponding root $`k_n`$ can be obtained via

$$k_n=_{\widehat{k}_{n1}}^{\widehat{k}_n}k\rho (k)𝑑k.$$
(13)
Using the periodic orbit expansion for the density of states (8), the integration (13) can be performed explicitly and yields

$`k_n={\displaystyle \frac{\pi }{S_0}}n{\displaystyle \frac{2}{\pi }}{\displaystyle \underset{p}{}}{\displaystyle \frac{1}{S_p^0}}{\displaystyle \underset{\nu =1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{A_p^\nu }{\nu ^2}}\mathrm{sin}\left({\displaystyle \frac{1}{2}}\nu \omega _p\right)\mathrm{sin}(\nu \omega _pn),`$ (14)
where $`\omega _p=\pi S_p^0/S_0`$, and the $`A_p`$’s are assumed to be real (no mixed boundary conditions).

The spectral formula (14) provides an explicit harmonic expansion of the function $`\phi (n)`$ in (4) for the roots of the spectral determinant of the regular graphs.

It can be shown  that the series (14) is convergent. An analytical criterion for the regularity (12) of a given network comes from analyzing the spectral determinant (2). After removing the complex phase of $`\mathrm{\Delta }(k)`$, the spectral equation for quantum graphs can be written as $`\mathrm{cos}\left(S_0k\pi \gamma _0\right)=\mathrm{\Phi }(k)`$. It was shown in , that if the characteristic function $`\mathrm{\Phi }(k)`$ of the graph is bounded by $`1`$ for all $`kR^1`$, then the piercing average condition (12) is satisfied  and the expansion (14) exists. In terms of the coefficients of (2), $`|\mathrm{\Phi }(k)|<1`$ is certainly satisfied if

$$\underset{i=1}{\overset{N_\mathrm{\Gamma }}{}}|a_i|<1.$$
(15)
In  we called graphs satisfying (15) regular quantum graphs. In general, the regularity condition is satisfied only for a special choice of the potentials (hence the importance of using various “graph dressings”) .

## IV Linear chain graphs

For many network topologies there exist no regular regimes at all . As an example of graphs which can be made regular by an appropriate choice of parameters, let us consider a dressed chain graph (a sequence of $`1D`$ square wells), with the vertices positioned at $`x=b_0,x=b_1,\mathrm{},x=b_n`$ starting at $`b_0=0`$. For every region $`b_{i1}<x<b_i`$, the wave function is given by

$$\psi _i(x)=\frac{A_i}{\sqrt{k_i}}e^{ik_ix}+\frac{B_i}{\sqrt{k_i}}e^{ik_ix},$$
(16)
where $`k_i=\sqrt{k^2\lambda _ik^2}=\beta _ik`$, and where the physically meaningful flux amplitudes $`1/\sqrt{k_i}`$ are factored out of the arbitrary coefficients $`A_i`$ and $`B_i`$. The boundary conditions at every vertex,

$`{\displaystyle \frac{A_i}{\sqrt{\beta _i}}}e^{ik\beta _ib_i}+{\displaystyle \frac{B_i}{\sqrt{\beta _i}}}e^{ik\beta _ib_i}`$ $`=`$ $`{\displaystyle \frac{A_{i+1}}{\sqrt{\beta _{i+1}}}}e^{ik\beta _{i+1}b_i}+{\displaystyle \frac{B_{i+1}}{\sqrt{\beta _{i+1}}}}e^{ik\beta _{i+1}b_i},`$
$`A_i\sqrt{\beta _i}e^{ik\beta _ib_i}B_i\sqrt{\beta _i}e^{ik\beta _ib_i}`$ $`=`$ $`A_{i+1}\sqrt{\beta _{i+1}}e^{ik\beta _{i+1}b_i}B_{i+1}\sqrt{\beta _{i+1}}e^{ik\beta _{i+1}b_i}`$
can be written in the matrix form

$$M_{i+1,i}\stackrel{}{C}_{i+1}=M_{i,i}\stackrel{}{C}_i,$$
(17)
where

$$M_{ij}\left(\begin{array}{cc}\frac{1}{\sqrt{\beta _i}}e^{ik\beta _ib_j}& \frac{1}{\sqrt{\beta _i}}e^{ik\beta _ib_j}\\ \sqrt{\beta _i}e^{ik\beta _ib_j}& \sqrt{\beta _i}e^{ik\beta _ib_j}\end{array}\right),\stackrel{}{C}_i\left(\begin{array}{c}A_i\\ B_i\end{array}\right).$$
(18)
The amplitudes $`\stackrel{}{C}_2`$ at a vertex $`V_i`$ can be expressed via the amplitudes at the previous vertex by means of the transfer matrix $`T_i`$,

$$\stackrel{}{C}_{i+1}=M_{i+1,i}^1M_{i,i}\stackrel{}{C}_iT_i\stackrel{}{C}_i,$$
(19)
where

$$T_i=\frac{1}{t_i}\left(\begin{array}{cc}e^{ik\left(\beta _i\beta _{i+1}\right)b_i}& r_ie^{ik\left(\beta _{i+1}+\beta _i\right)b_i}\\ r_ie^{ik\left(\beta _{i+1}+\beta _i\right)b_i}& e^{ik\left(\beta _{i+1}\beta _i\right)b_i}\end{array}\right)$$
(20)
is unitary, $`detT_i=1`$. Here $`r_i`$ is the reflection and $`t_i=\sqrt{1r_i^2}`$ is the transmission coefficient,

$`r_i`$ $`=`$ $`{\displaystyle \frac{\beta _{i+1}\beta _i}{\beta _{i+1}+\beta _i}},`$
$`t_i`$ $`=`$ $`{\displaystyle \frac{2\sqrt{\beta _i\beta _{i+1}}}{\beta _i+\beta _{i+1}}}.`$
Applying the “boundary matching” equation (19) consecutively to all vertices of the chain yields the equation

$$\stackrel{}{C}_n=T\stackrel{}{C}_1,$$
(21)
where

$$T=T_{n1}T_{n2}\mathrm{}T_1.$$
(22)
Due to the Dirichlet boundary conditions at the end vertices $`V_0`$ and $`V_n`$ of the chain,

$$\stackrel{}{C}_1=\left(\begin{array}{c}A_1\\ A_1\end{array}\right),\stackrel{}{C}_n=\left(\begin{array}{c}A{}_{n}{}^{}e_{}^{ik_nb_n}\\ A{}_{n}{}^{}e_{}^{ik_nb_n}\end{array}\right),$$
(23)
equation (21) can be written in terms of the matrix elements of $`T`$ only,

$$e^{ik_nb_n}\left[t_{11}(k)t_{12}(k)\right]=e^{ik_nb_n}\left[t_{21}(k)t_{22}(k)\right].$$
(24)
Since the spectral equation (24) is linear, the common factor

$$\tau ^1=\underset{i=1}{\overset{n}{}}t_i^1$$
(25)
that multiplies the matrix $`T`$ (22) cancels out, so that (24) contains only the reflection coefficients $`r_1`$, $`r_2`$,…, $`r_n`$, and no transmission coefficients $`t_i`$. From the explicit form of $`T_i`$ it follows that the matrix elements $`t_{ij}`$ of $`T`$,

$`\left({\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{1}{t_i}}\right)\left(\begin{array}{cc}e^{ik_nb_n}e^{ikS_0}+e^{ik_nb_n}P(r_i,k)& e^{ik_nb_n}Q(r_i,k)\\ e^{ik_nb_n}Q^{}(r_i,k)& e^{ik_nb_n}e^{ikS_0}+e^{ik_nb_n}P^{}(r_i,k)\end{array}\right),`$

are certain multivariable polynomials of the variables $`r_i`$ that may include monomials of all the degrees between $`1`$ and $`n`$. Hence the coefficients $`a_i`$ in the spectral determinant

$`\mathrm{\Delta }(k)=e^{2ik_nb_n}\left(t_{11}(k)t_{12}(k)\right)+t_{21}(k)t_{22}(k)`$ (26)
$`=1+e^{i2(S_0k\pi \gamma _0)}{\displaystyle \underset{i=1}{\overset{2N_\mathrm{\Gamma }}{}}}a_ie^{i2(S_ik\pi \gamma _i)},`$ (27)
are also polynomials in $`r_i`$ of degree larger or equal to $`1`$.

Hence, by choosing the graph dressing parameters $`\beta _i`$ sufficiently close to one another, one can make the coefficients $`r_i`$ as small as necessary in order to satisfy the regularity condition (15). Hence the regular regime (12) exists for all the linear chain graphs. Physically this corresponds to the case where the heights of the bond potentials do not differ too much from one another.

As an example, consider the case of a two-bond linear graph with Dirichlet boundary conditions at the end vertices positioned at $`x=0`$ and $`x=b_2`$ . In this case the transfer matrix is

$$T=\frac{1}{t}\left(\begin{array}{cc}e^{ik\left(\beta _1\beta _2\right)b_1}& re^{ik\left(\beta _1+\beta _2\right)b_1}\\ re^{ik\left(\beta _1+\beta _2\right)b_1}& e^{ik\left(\beta _2\beta _1\right)b_1}\end{array}\right),$$
(28)
where $`r`$ and $`t`$ are correspondingly the reflection and the transmission coefficients at the middle vertex $`V_2`$. The spectral equation obtained from (28) is

$`\mathrm{sin}\left[k\left(\beta _1b_1+\beta _2\left(b_2b_1\right)\right)\right]+r\mathrm{sin}\left[k\left(\beta _1b_1\beta _2\left(b_2b_1\right)\right)\right]=0.`$ (29)
Since $`r1`$, the condition (15) is satisfied and therefore a two-bond chain graph is always regular .

## V Irregular regimes

As a network moves from the regular to the irregular regime, the principle of “one root per interval $`I_n`$” is violated. In the irregular regime certain intervals contain no roots of the spectral equation at all, while others may host two or more roots. The piercing property of the Weyl average  also disappears, so in order to separate the intervals that contain exactly one root, one cannot use Weyl’s average and the intersection points (12).

As outlined in , in order to repeat the root expansion procedure discussed above, one would need to find some smooth monotonic function $`\overline{N}_1(k)`$, which pierces every stair of $`N(k)`$, and find the intersection points $`\overline{N}_1(\widehat{k}_n^{})=N(\widehat{k}_n^{})`$. Solving the equation $`\overline{N}_1(\widehat{k}_n^{})=n`$ would lead to the “generalized separators” $`\widehat{k}_n^{}`$, that isolate the roots $`k_n`$ from each other. Unlike in the regular regime, where the separators $`\widehat{k}_n^{}`$ are periodic, the separators $`\widehat{k}_n^{}`$ form an almost periodic set .

## VI Discussion

Since quantum chaos theory deals with non-integrability in the semiclassical regime, the systems it usually addresses are fairly complicated. One of the drawbacks of such complexity is that the various methods used in the theory (such as periodic orbit theory, random matrix theory, etc.) describe the system from quite different perspectives, and the direct connections between these methods are often obscured.

Quantum networks provide a convenient exception to this. For these systems, many essential attributes of quantum chaos theory, and in particular the periodic orbit theory, can be directly established, and the results turn out to be exact. As a result, quantum networks provide a convenient model for testing various phenomenological and approximate methods.

Due to their simplicity, analytical solvability and overall transparency, quantum graphs play the role of the “harmonic oscillator” of the quantum chaos theory.

## Acknowledgments

Y.D. and R.B. gratefully acknowledge financial support by NSF grant PHY-9984075; Y.D. and R.V.J by NSF grant PHY-9900746.