text:
# String order in spin liquid phases of spin ladders

## I Introduction

The properties of spin ladder systems attracted considerable interest in the past decade.review  It was shownferro  that even-leg ladder models develop a gap in their excitation spectrum if the rung coupling is ferromagnetic. This result could be considered as a generalization of Haldane’s conjucturehaldane  about the existence of a finite gap in integer spin Heisenberg chains, since for strong enough ferromagnetic rung coupling the spin ladder with spin-$`S`$ on the legs is equivalent to a spin chain of spin-$`2S`$. In this case the ground state is an AKLT like stateAKLT  in which short range valence bonds couple spins on neighboring rungs. Later it was demonstrated numerically,antiferro  that a gap is generated for antiferromagnetic rung coupling as well, in which case the singlet valence bonds are formed predominantly on the rungs. It may be called a rung singlet state. Both are Haldane-type spin liquids in which there are no broken local symmetries, and the excitations are coherent magnons with a finite gap.

Similar situation occurs if the two legs of the ladder are coupled not by rung couplings, but diagonally or by plaquette couplings.legeza  When these couplings compete with each other, non-Haldane spin liquid phases, characterized by the absence of coherent magnon excitations and spontaneous dimerization,legeza ; nersesyan may occur due to frustration effects. Although the thermodynamic properties of Haldane and non-Haldane spin liquids are identical, the correlation functions differ due to the difference in the excitations spectrum. The two phases could also be distinguished by a local dimer order parameter.

Beside these spin liquid phases several other phases may occur, like ferromagnetic or incommensurate phases. When drawing the phase diagram, it was thought for some time that the concept of Haldane-like spin liquids is unambiguous for spin ladders, and therefore the AKLT-like and rung singlet phases are identical not only thermodynamically, but also in the sense, that in a large parameter space, where rung, diagonal and plaquette couplings are all taken into account, one could go from one to the other without any discontinuity. In fact, it was arguedwhite1  that there is a controlled way of changing the parameters gradually and interpolate continuously between two limiting cases belonging trivially to one or to the other phase. The procedure contains a shift of one of the legs of the ladder by one unit, thereby transforming a rung singlet into a diagonally situated valence bond, characteristic of the the AKLT-like state, or vice versa. However, as was pointed out by Kim et al.kim  this shift changes the topology of the system drastically. When we start from an open ladder, the above transformation deforms the usual straight-end ladder into a “ladder” whose ends are cut diagonally. Such a deformation has serious consequences on the degeneracy of the ground state for open boundary condition. In a usual ladder the AKLT state has a fourfold degeneracy due to the unpaired end spins, while the rung singlet state is unique. On the other hand, in the diagonally cut ladder the ground state is unique for the AKLT state, and fourfold degenerate in the rung singlet state. This is best seen in Fig. 1(a) and (b), where the typical valence-bond configurations of the AKLT-like state and the rung singlet state are shown. In Fig. 1(c) we also show a typical valence bond configuration of the dimer phase. In this phase the invariance of translation by one rung is broken spontaneously.

While in the AKLT state the number of valence bonds crossing an arbitrary vertical line between rungs is always odd, in the rung singlet state it is always even. This is true not only in this special configuration but remains valid generally, as long as the range of valence bonds is short compared to the size of the system.bonesteel  The two states can be distinguished by a topological quantum number. Thus for topological reasons in a spin ladder the two Haldane-type spin liquid phases are different and in the space of couplings they must be separated by a phase transition hyperplane,kim ; nedelcu  although extended gapless (critical) regions between them cannot be excluded a priori either.

Since the topological quantum number is not accessible to direct computation, a generalization of the hidden order parameter of the spin-1 chaindennijs  was proposed,kim  that could distinguish between the even and odd topological sectors of the Haldane phase. Although the relationship between the two new string order parameters and the two types of Haldane phases was made plausible by looking at special configurations, it was not possible to calculate their value, since the bosonization procedure applied to obtain the phase transition lines did not distinguish between the two parameters.

The aim of this paper is to do a numerical study of the relationship between the possible spin liquid phases of the spin ladders and the string order parameter, using the density matrix renormalization group (DMRG) procedure.dmrg  The setup of the paper is as follows. In Sec. II we give a short description of the ladder models and a brief account of what is known about the possible phases. The results of our numerical calculations are presented in Sec. III. Finally Sec. IV contains a brief summary.

## II Gapped phases in ladder models

Following the notations of Ref. IV we will write the most general form of the Hamiltonian for isotropic two-leg spin ladder models if only spins on the same and neighboring rungs interact, as

$$=\underset{i=1}{\overset{N}{}}h_{i,i+1},$$
(1)
where

$`h_{i,i+1}`$ $`=`$ $`J_\mathrm{l}\stackrel{}{\sigma }_i\stackrel{}{\sigma }_{i+1}+J_\mathrm{l}^{}\stackrel{}{\tau }_i\stackrel{}{\tau }_{i+1}`$ (2)
$`+{\displaystyle \frac{1}{2}}J_\mathrm{r}\stackrel{}{\sigma }_i\stackrel{}{\tau }_i+{\displaystyle \frac{1}{2}}J_\mathrm{r}^{}\stackrel{}{\sigma }_{i+1}\stackrel{}{\tau }_{i+1}`$
$`+J_\mathrm{d}\stackrel{}{\sigma }_i\stackrel{}{\tau }_{i+1}+J_\mathrm{d}^{}\stackrel{}{\tau }_i\stackrel{}{\sigma }_{i+1}`$
$`+V_{\mathrm{ll}}(\stackrel{}{\sigma }_i\stackrel{}{\sigma }_{i+1})(\stackrel{}{\tau }_i\stackrel{}{\tau }_{i+1})`$
$`+V_{\mathrm{rr}}(\stackrel{}{\sigma }_i\stackrel{}{\tau }_i)(\stackrel{}{\sigma }_{i+1}\stackrel{}{\tau }_{i+1})`$
$`+V_{\mathrm{dd}}(\stackrel{}{\sigma }_i\stackrel{}{\tau }_{i+1})(\stackrel{}{\tau }_i\stackrel{}{\sigma }_{i+1}).`$
In Eq. (2) $`\sigma _i`$ and $`\tau _i`$ are spin-1/2 operators corresponding to spins sitting on the two legs of the ladder. The schematic plot of the spin couplings in $`h_{i,i+1}`$ between the spins on the two legs of the ladder are shown in Fig. 2.

The terms with $`J_\mathrm{l}`$ and $`J_\mathrm{l}^{}`$ couple neighbouring spins on the same legs, the terms with $`J_\mathrm{r}`$ and $`J_\mathrm{r}^{}`$ spins on the same rung, while the terms with $`J_\mathrm{d}`$ and $`J_\mathrm{d}^{}`$ spins that are situated diagonally. In addition to these bilinear terms there could be three biquadratic terms with couplings $`V_{\mathrm{ll}}`$, $`V_{\mathrm{rr}}`$ and $`V_{\mathrm{dd}}`$.

In order to decrease the number of parameters from now on we only consider the case when $`J_\mathrm{R}J_\mathrm{r}=J_\mathrm{r}^{}`$, $`J_\mathrm{L}J_\mathrm{l}=J_\mathrm{l}^{}`$ and $`J_\mathrm{D}J_\mathrm{d}=J_\mathrm{d}^{}`$. The energy scale will be set by fixing $`J_\mathrm{L}=1`$.

It is known from previous studieslegeza  that all the interleg couplings $`J_\mathrm{R},J_\mathrm{D},\mathrm{},V_{\mathrm{dd}}`$, when acting alone, generate a Haldane-like phase. This is an AKLT-like phase with odd topology of the number of valence bonds if $`J_\mathrm{R}<0`$, $`J_\mathrm{D}>0`$, $`V_{\mathrm{ll}}>0`$, $`V_{\mathrm{rr}}<0`$ or $`V_{\mathrm{dd}}>0`$, while for opposite sign of the couplings a rung singlet like state with even topology is generated.

As mentioned before, the goal of the present paper is to identify an appropriate order parameter that can characterize the spin liquid phases. Starting from the string order parameter of the $`S=1`$ chain proposed by den Nijs and Rommelsedennijs , it is natural to generalize it to spin ladders by defining the order parameter of AKLT-like Haldane state as

$$𝒪_{\mathrm{odd}}^\alpha =\underset{|ij|\mathrm{}}{lim}S_i^\alpha \mathrm{exp}\left(i\pi \underset{l=i+1}{\overset{j1}{}}S_l^\alpha \right)S_j^\alpha $$
(3)
where $`S_i^\alpha `$ is the total spin on rung $`i`$, $`S_i^\alpha =\sigma _i^\alpha +\tau _i^\alpha `$, $`\alpha =x,y,z`$.

Since in the extreme limits the pure rung singlet state can be related to the AKLT state by a shift of one of the legs, the string order parameter of the Haldane state with even topological quantum number is expected to be

$$𝒪_{\mathrm{even}}^\alpha =\underset{|ij|\mathrm{}}{lim}S_i^\alpha \mathrm{exp}\left(i\pi \underset{l=i+1}{\overset{j1}{}}S_l^\alpha \right)S_j^\alpha $$
(4)
where $`S_i^\alpha `$ is now the sum of two diagonally situated spins, $`S_i^\alpha =\sigma _i^\alpha +\tau _{i+1}^\alpha `$.kim ; nishi

In what follows we will calculate the order parameters for the $`z`$ component only, they will be denoted as $`𝒪_{\mathrm{odd}}`$, and $`𝒪_{\mathrm{even}}`$, respectively. We will show, that these two order parameters are mutually exclusive. $`𝒪_{\mathrm{odd}}`$ is finite for the AKLT-like Haldane state, while $`𝒪_{\mathrm{even}}`$ is finite for the rung singlet like Haldane phase. Moreover the order parameter goes to zero continuously at the transition point when the transition is of second order, while it vanishes with a jump in case of first order transitions.

## III Numerical results

To investigate this problem we performed numerical calculations applying the DMRG methoddmrg  on the model defined by the Hamiltonian in Eq. (1). Since the DMRG procedure is more accurate for systems with free ends, we consider our ladder model with open boundary condition. In order to eliminate the fourfold degeneracy of the ground state in the conventional Haldane phase due to the free $`S=1/2`$ spins at the ends of the ladder, we exploited the left-right reflection symmetry of the system. Out of the four degenerate ground state we have selected the one with $`S_{\mathrm{tot}}^z=0`$ and positive parity. The calculations were done for chains with up to 150 sites keeping 100 to 200 states per block. The truncation error varied in the range $`10^5`$ to $`10^7`$. In most of the calculations we used the finite-lattice method with two or three iteration cycles.

In order to reduce the undesirable boundary effects we fixed site $`i`$ in Eqs. (3) and (4) at the center of the chain and let site $`j`$ run through half of the chain. The finite size scaling procedure was used to determine the thermodynamic limit of the quantitities measured.

In what follows we present our numerical results for the string order parameters for various choices of the coupling constants. Since our aim was not to explore the full phase diagram, but rather to check how the string order parameters $`𝒪_{\mathrm{even}}`$ and $`𝒪_{\mathrm{odd}}`$ vary as we go from one Haldane-like phase to another, topologically different gapped phase, we calculated the string order parameters along specially chosen paths connecting points in the parameter space, where the model is known to belong to different phases.

It is worth pointing out that since the ground state wave function and therefore the ground state expectation values of operators can be calculated more accurately than those of excited states, the order parameters provide a better way to obtain phase diagrams than investigating the energy spectrum.

### III.1 Haldane phases due to rung, diagonal or plaquette couplings alone

The effect of $`J_\mathrm{R}`$ alone was analysed numerically for short systems in Ref. IV. It was found there that for $`J_\mathrm{R}<0`$ the order parameter $`𝒪_{\mathrm{odd}}`$ is finite, while for $`J_\mathrm{R}>0`$ the other order parameter $`𝒪_{\mathrm{even}}`$ is finite. Our results with the DMRG method confirmed these results.

From earlier workslegeza  it is known, that not only the rung but also the diagonal and plaquette inter-leg couplings are relevant perturbations which generate a gap and drive the system into phases with topological order. Depending on the sign of the perturbation the emerging phases belong to one of the two possible topological sectors.

We have therefore investigated the effect of these couplings by doing calculations at two specialy chosen values. The finite size scaling analysis of the results obtained for couplings $`J_\mathrm{d}=\pm 1`$ or $`J_\mathrm{D}=\pm 1`$ indicated that irrespective whether only one or both diagonal couplings are present, $`𝒪_{\mathrm{odd}}0`$ and $`𝒪_{\mathrm{even}}=0`$ for positive $`J_\mathrm{d}`$ or $`J_\mathrm{D}`$, while we have $`𝒪_{\mathrm{even}}0`$ and $`𝒪_{\mathrm{odd}}=0`$ for negative $`J_\mathrm{d}`$ or $`J_\mathrm{D}`$.

Similarly, calculations for the values $`\pm 1`$ of the couplings $`V_{\mathrm{ll}}`$, $`V_{\mathrm{rr}}`$ and $`V_{\mathrm{dd}}`$ (the other couplings set to zero) convinced us that the even order parameter is finite, $`𝒪_{\mathrm{even}}>0`$, for $`V_{\mathrm{ll}}<0`$, or $`V_{\mathrm{rr}}>0`$ or $`V_{\mathrm{dd}}<0`$, while the odd order parameter is finite, $`𝒪_{\mathrm{odd}}>0`$, for opposite signs of the couplings, i.e., for $`V_{\mathrm{ll}}>0`$, or $`V_{\mathrm{rr}}<0`$ or $`V_{\mathrm{dd}}>0`$.

Our results obtained for the string order parameters allows us to conclude that the topology of the valence bond structure can be identified directly by calculating the order parameters $`𝒪_{\mathrm{odd}}`$ or $`𝒪_{\mathrm{even}}`$, since they mutually exclude each other.

### III.2 The $`J_\mathrm{R}`$ \- $`J_\mathrm{D}`$ phase diagram

Next we consider the subspace spanned by $`J_\mathrm{R}`$ and $`J_\mathrm{D}`$. According to the recent bosonization studykim  the phase transition line separating the two topological sectors is of first order in the quadrant $`J_\mathrm{R}>0`$, $`J_\mathrm{D}>0`$, while of second order in the quadrant $`J_\mathrm{R}<0`$, $`J_\mathrm{D}<0`$. In order to check the behavior of the string order parameter, we have calculated them along two paths in the $`(J_\mathrm{R},J_\mathrm{D})`$ parameter space.

Path 1 is defined by $`J_\mathrm{R}=1J_\mathrm{D}`$ with $`0<J_\mathrm{D}<1`$, and path 2 is parametrized by $`J_\mathrm{R}=1J_\mathrm{D}`$ with $`1<J_\mathrm{D}<0`$. They are indicated by solid lines on Fig. 3.

Our numerical results for $`𝒪_{\mathrm{even}}`$ and $`𝒪_{\mathrm{odd}}`$ obtained in the thermodynamic limit along path 1 is plotted on Fig. 4. The lines connecting the points are only guides to the eye. It is apparent from the figure that $`𝒪_{\mathrm{even}}`$ is finite in one particular region while $`𝒪_{\mathrm{odd}}`$ in the other.

As one sees from the figure the string order parameters, extrapolated to the infinite ladder limit, have a jump discontinuity at the transition point. In order to further support the first order nature of the phase transition we have also investigated the low lying energy spectrum. We found that there is a sudden change in the degeneracy of the ground state at the transition point. In the case of open ladders the fourfold degeneracy of the AKLT phase associated with spin-1/2 end-spins disappear, and the ground state becomes unique as a rung-dimer-like structure emerges. Moreover, by examining higher lying energy levels we found that the energy gap remains finite throughout the transition. The ground state energy shows a cusp-like singularity, a characteristic feature of first order transitions, as is shown on the inset of Fig. 4. Note, however, that for finite chains there is no observable level crossing in the ground state due to hybridization.

Next we consider the behavior of the order parameters calculated along path 2 in the negative $`J_\mathrm{R}`$ and $`J_\mathrm{D}`$ quadrant. Our results are plotted in Fig. 5. In this case again, only one of the order parameters is finite in a particular phase, however, the gap vanishes continuously at the phase boundary, which is a characteristic features of second order phase transitions. It is worth noting that there is an order of magnitude difference in the scale of Fig. 4 and Fig. 5 as a function of $`J_\mathrm{D}`$.

Since the phase transition points were located on two special trajectories only, the phase boundary that separates the two types of Haldane-like phases is drawn by dashed line in Fig. 3 to indicate that its exact shape was not determined.

### III.3 The role of plaquette couplings

In order to consider the influence of the plaquette couplings we have chosen the special values $`KJ_\mathrm{R}=V_{\mathrm{ll}}=V_{\mathrm{dd}}`$, $`V_{\mathrm{rr}}=0`$, since the phase diagram has already been studied for this parametrization.legeza  This choice has been motivated by the possible mapping of the ladder model to the spin-1 bilinear-biquadratic chain for $`K=1`$.legeza  Here we calculate string order parameters along four lines going through the expected phase boundaries shown by solid lines in Fig. 6.

Calculations along path 1 have revealed that $`𝒪_{\mathrm{odd}}`$ is finite in the AKLT-like Haldane phase, but it vanishes in the dimerized phase. On the other hand, $`𝒪_{\mathrm{even}}`$ becomes finite in the dimerized phase, thus the two phases are well distinguished by the order parameters. Examining the trajectory (path 2) going through the phase boundary between the rung singlet like Haldane and dimerized phases, we have found that $`𝒪_{\mathrm{even}}`$ remains finite and $`𝒪_{\mathrm{odd}}`$ scales to zero everywhere. This shows that the string order parameter does not uniquely characterize the Haldane phase, it is rather an indicator of the topology of the valence bond structure, since both the rung singlet and the dimerized phases have even topological quantum numbers.

Thus when $`𝒪_{\mathrm{even}}0`$, the string order is not sufficient to classify unambiguously the ground state. Another order parameter, the so-called dimer order parameter has to be measured. Taking the difference of the local energy on neighboring bonds in the middle of the ladder this order parameter is defined as

$$D=\frac{h_{i,i+1}h_{i+1,i+2}}{\frac{1}{2}\left[h_{i,i+1}+h_{i+1,i+2}\right]}.$$
(5)
In the Haldane-type spin liquids $`D`$ vanishes, but there exists an extended region in the parameter space, where the dimer order parameter scales to finite value in the thermodynamic limit, signaling a dimerized phase with spontaneously broken translational symmetry.

Next let us consider the $`K>0`$ half-plane focusing on the $`J_\mathrm{D}=0`$ half-line and its vicinity. Calculation of the order parameters along path 3, defined by $`K=1+J_\mathrm{D}`$ with $`1<J_\mathrm{D}<1`$ revealed that $`𝒪_{\mathrm{even}}`$ is finite for $`J_\mathrm{D}<0`$, while $`𝒪_{\mathrm{odd}}`$ is finite for $`J_\mathrm{D}>0`$. At $`J_\mathrm{D}=0`$ all order parameters including $`D`$ vanish as is shown in Fig. 7. Since the string order parameters become finite for arbitrary small finite value of $`J_\mathrm{D}`$ the two types of spin liquids are separated by a critical line and the extended critical phase reported earlierlegeza  does not exist.

Finally, it can be seen on the figure that for $`J_\mathrm{D}>0.6`$ both order parameters vanish and a new phase appears. In order to analyze the question whether the new phase develops just above the $`VBS`$ point, we investigated the vicinity of the $`VBS`$ point along path 4. We found that $`𝒪_{\mathrm{odd}}`$ remains finite up to $`K1.9`$ while $`𝒪_{\mathrm{even}}`$ vanishes for the whole regime. Investigating the low lying energy spectra we found a multiply degenerate ground state. In order to decide whether this phase extends to the $`J_\mathrm{D}`$ axis we made test calculations at $`J_\mathrm{D}=0.05`$, $`K=1.5,2`$ and found a finite $`𝒪_{\mathrm{odd}}`$ but vanishing $`𝒪_{\mathrm{even}}`$ order parameters. This suggests the schematic phase boundary shown by a dashed line in the figure.

## IV Conclusions

In the present paper we have studied the relationship between various gapped Haldane and non-Haldane spin liquid phases and the two different string order parameters in the general two-leg spin ladder model, containing rung, diagonal and plaquette couplings as well. The string order parameters have been calculated numerically using the DMRG method.

It has been found that there are two kinds of Haldane-type massive spin liquids characterized by finite $`O_{\mathrm{even}}`$ or $`O_{\mathrm{odd}}`$ string order parameters, respectively. It was shown that $`𝒪_{\mathrm{odd}}`$ is finite in the AKLT-like phase, while $`𝒪_{\mathrm{even}}`$ is finite in the rung singlet like phase. Examining the transition from the AKLT-like to the rung-singlet-like phases the type of phase transition can be first or second order depending on the path in parameter space, and this shows up in the abrupt or continuous vanishing of the order parameter.

When plaquette couplings are allowed for, non-Haldane-like spin liquid phases can occur, and the string order parameters are not sufficient to characterize the phases. $`𝒪_{\mathrm{even}}`$ remains finite in the dimerized phases, too. Thus the string order parameter reflects the topology of the valence bond structure and does not characterize uniquely the Haldane-type phases. A dimer order parameter can, however, be defined, which is finite in a non-Haldane-type dimer phases but vanishes in Haldane-type spin liquids.

It was found previously that there exists a small range of the parameters where the combination of perturbing operators is irrelevant and the spectrum remains massless. Investigating the vicinity of $`J_\mathrm{D}=0,K>0`$ axis we have found that both string order parameters vanish on this critical line, but $`𝒪_{\mathrm{even}}`$ or $`𝒪_{\mathrm{odd}}`$ becomes finite immediately as we move away from this axis. Therefore, we conclude that the two kind of spin liquids are separated by a critical line and not by an extended gapless regime reported earlier. In addition, around the $`J_\mathrm{D}=1`$ line for $`K>1.8`$ we have found a new phase. This phase is characterized by multiply degenerate ground state and vanishing order parameters. The detailed analysis of this region is left for future investigation.

###### Acknowledgements.

This research was supported in part by the Hungarian Research Fund (OTKA) Grant No. T30173, F31949 and F32231. G. F. was also supported by the Bolyai Research Scholarship. Ö. L. acknowledges the hospitality of the Friedrich Alexander University of Nürnberg where part of the work was completed.


format: text

text:
# Flavor altering excitations of composite fermions

## Abstract

Past theoretical studies have considered excitations of a given flavor of composite fermions across composite-fermion quasi-Landau levels. We show that in general there exists a ladder of flavor changing excitations in which composite fermions shed none, some, or all of their vortices. The lowest energy excitations are obtained when the composite fermions do not change their flavor, whereas in the highest energy excitations they are stripped of all of their vortices, emerging as electrons in the final state. The results are relevant to the intriguing experimental discoveryPinczuk  of Hirjibehedin et al. of coexisting excitation modes of composite fermions of different flavor in the filling factor range $`1/3>\nu 1/5`$.

The quantum fluid of interacting electrons in the lowest Landau level, which manifests itself most dramatically through the fractional quantum Hall effectStormer , is characterized by the formation of a new class of fermions called composite fermionsJain  (CFs). The composite fermion is the bound state of an electron and an even number ($`2p`$) of quantized vortices, denoted by 2pCF. Composite fermions of different flavors (that is, with different numbers of attached vortices) occur in different filling factor regions, where the filling factor, $`\nu =\rho hc/eB`$ is the number of occupied Landau levels (LLs). (Here $`\rho `$ is the two-dimensional density of electrons and $`B`$ is the external magnetic field.) Specifically, 2pCFs are relevant in the range $`\frac{1}{2p1}>\nu \frac{1}{2p+1}`$.

Light scattering has proved to be a powerful tool for investigating the nature of excitations of this quantum fluid RamanPinczuk ; RamanKang ; RamanDujovne ; RamanDavies . The present work has been motivated by the recent light scattering experiment by Hirjibehedin et al.Pinczuk , which investigated the filling factor range $`1/3\nu 1/5`$, where the relevant composite fermions carry four vortices (4CFs). A remarkable aspect of the experimental results is that while new low-energy excitations appear for $`\nu <1/3`$, as expected, the old excitations from $`\nu >1/3`$ do not disappear but rather evolve continuously as the filling factor changes from $`\nu >1/3`$ to $`\nu <1/3`$; in other words, there is a coexistence of excitations of both 4CFs and 2CFs in the low filling factor region.

We begin by presenting a physical explanation. A consideration of the filling factor region $`\nu >1/3`$ (Fig. 1) provides valuable insight into this intriguing observation. Here, we have the lowest electronic LL partially occupied, and, as a result of the repulsive interaction, electrons transform into 2CFs; the lowest electronic LL thus splits into quasi-Landau levelscomment  (QLLs) of 2CFs (middle column of Fig. 1). Now, the intra-LL excitations of electrons are described as inter-QLL excitations of 2CFs. However, the system obviously still supports the inter-electronic LL excitation, namely the Kohn mode. We see in this simple case the coexistence of excitation modes at two different energy scales: the low energy $`\mathrm{\Delta }^{22}(k)`$ excitations conserve the CF flavor, whereas the high energy $`\mathrm{\Delta }^{20}(k)`$ excitations involve transformation of 2CF into 0CF (electron). In the former case, the composite fermions maintain their integrity during the transition, whereas in the latter case they leave their vortices behind in the lowest electronic Landau level, thus converting into excited electrons in the second level.

For $`\nu <1/3`$ the situation is analogous. The lowest 2CF QLL undergoes a further splitting into QLLs of 4CFs, as shown in the rightmost column of Fig. 1. Although we have only 4CFs in the ground state, there are three ladders of excitation: (i) $`\mathrm{\Delta }^{44}(k)`$ excitations, which were not present for $`\nu >1/3`$; (ii) $`\mathrm{\Delta }^{42}(k)`$ excitations involving a change of 4CFs into 2CFs, which are a continuous evolution of $`\mathrm{\Delta }^{22}(k)`$ excitations for $`\nu >1/3`$; and (iii) $`\mathrm{\Delta }^{40}(k)`$, which are similarly close cousins of $`\mathrm{\Delta }^{20}(k)`$ and $`\mathrm{\Delta }^{00}(k)`$. This is our physical explanation for the experimental result of Ref. Pinczuk,  showing a coexistence of several kinds of excitations.

To test these ideas quantitatively through a microscopic calculation, we need to construct a model that can deal with both $`\mathrm{\Delta }^{44}`$ and $`\mathrm{\Delta }^{42}`$ excitations simultaneously. The key observation that makes this possible is that the intra-QLL excitations of a given flavor of composite fermions are the inter-QLL excitations of higher order composite fermions. For example, at filling factors of the type $`\nu =n/(2n+1)`$ (with $`n`$ integer), which are described in terms of 2CFs, the intra-LL excitations of electrons (0CFs) are inter-QLL excitations of 2CFs. If we diagonalize this problem within the lowest electronic LL, we will miss the Kohn mode. In order to keep this mode, we must expand our basis space to allow at least one electron to occupy the second electronic LL; then we get not only all the $`\mathrm{\Delta }^{22}`$ excitations, but also the lowest energy $`\mathrm{\Delta }^{20}`$ excitations. This idea can be transported straightforwardly to the filling factor region $`1/3>\nu 1/5`$ by the vortex attachment trick of the CF theory, which relates $`\nu ^{}`$ to $`\nu =\nu ^{}/(2\nu ^{}+1)`$. We first construct an electron basis at $`\nu ^{}`$, wherein we keep all basis states in which either all electrons are in the lowest LL or one electron occupies the second LL. By attaching two vortices to each particle, we generate from it, following methods outlined below, a 2CF basis at $`\nu `$ shown schematically in Figs. 2a and 2b. We then diagonalize the interaction Hamiltonian in this basis. An explicit calculation for $`N=5`$ and 6 for certain filling factors in the range $`1/3>\nu 1/5`$ gives the spectrum in Fig. 3 shown by filled circles.

In order to identify the states at $`\nu `$ that correspond to the Kohn mode at $`\nu ^{}`$, we repeat the above calculation without the Kohn mode (that is, by working strictly within the lowest LL at $`\nu ^{}`$; Fig. 2a). We find that this gives a subset of the previous states, shown by empty circles in Fig. 3. (All empty circles sit on filled circles; for a given eigenstate, the energy of the filled circle is slightly lower than that of the empty circle because of the greater variational freedom in the expanded basis, but the difference is so small in most cases that the empty circles fully cover the filled ones.) The states not covered by empty circles are therefore inter-QLL excitations of 2CFs, that is, $`\mathrm{\Delta }^{42}`$ excitations.

In short, our approach is to model the $`1/3>\nu 1/5`$ region in terms of 2CFs – as opposed to 4CFs, which would be the natural choice if we were only interested in the lowest energy excitations. Then, the intra-QLL excitations of 2CFs give us the inter-QLL excitations of 4CFs, namely $`\mathrm{\Delta }^{44}`$; at the same time the inter-QLL excitations of 2CFs correspond to the $`\mathrm{\Delta }^{42}`$ excitations.

There is a remarkable separation of energy scales in Fig. 3 between the inter- and intra-QLL excitations of 2CFs, which is entirely consistent with the physics described above. A similar effect exists for the inter- and intra-LL excitations of electrons, which are governed by different energies, namely the cyclotron energy and the interaction energy. For 2CFs, however, both the inter- and intra-QLL energies are determined by the interaction energy, and the emergence of two energy scales is not obvious a priori. Fig. 3 also shows exact diagonalization results (dashes). There is an excellent agreement between the dashes and the circles, further confirming the validity of our picture. In particular, the exact spectra also show two distinct energy scales. (In fact, once we understand the physics, we can identify the $`44`$ and $`42`$ energies directly from the exact spectrum, which is what we do below.) Still higher energy states can be obtained within the framework of the CF theory by exciting more than one composite fermion to higher CF quasi-Landau levels, but our interest here is in the lowest excitation at each flavor changing ladder.

We next compare our results with experiment. In an ideal situation, inelastic light scattering only probes the excitations close to zero wave vector, but slight disorder, which breaks wave vector conservation, has been found to expose energies at critical points (extrema) in the CF exciton dispersion RamanPinczuk ; RamanKang ; RamanDujovne . In particular, the zero wave vector mode, the minimum energy neutral excitation (roton), and the excitation that contains a far separated particle hole pair of composite fermions (which gives the transport gap measured in experiments) have been observed in Raman experiments. We will concentrate on the first two, which will be labeled $`\mathrm{\Delta }(0)`$ and $`\mathrm{\Delta }(R)`$, respectively. The transport gap is roughly the same as $`\mathrm{\Delta }(R)`$ for our small systems.

The lowest energy intra- and inter-2CF QLL excitations can be identified straightforwardly from the spectra of Fig. 3. The upper panel of Fig. 4 shows the energies $`\mathrm{\Delta }^{44}(0)`$, $`\mathrm{\Delta }^{44}(R)`$, and their $`42`$ counterparts. While the $`42`$ modes evolve continuously out of the $`22`$ modes at $`\nu 1/3`$, the $`44`$ modes are new. The qualitative behavior is in agreement with experiment. There are, however, several difficulties in making a quantitative comparison with experiment.

First, the calculations are performed for a strictly two dimensional system, whereas the experimental sample is a square quantum well of width 33 nm. We incorporate finite thickness corrections into our calculation through a self consistent local density approximationOrtalano . The lower panel of Fig. 4 depicts the various energies after including modification in the effective two-dimensional interaction due to finite extent of the electron wave function in the transverse direction.

The other problems are more severe. Strictly speaking, one must calculate thermodynamic limits before comparing theoretical calculations to experiment. In the past, energy gaps have been found to vary substantially as a function of $`N`$. (For example, the transport gap at $`\nu =2/5`$ is reduced from 0.13 for $`N=4`$ to $`0.058`$ in the thermodynamic limit Manual .) However, because of the rather large dimension of the CF basis, we are not able to go to systems bigger than $`N=6`$ in a numerically stable manner. The calculation also neglects disorder, which is often a source of significant quantitative correction.

Next we compare our theoretical results for the roton energy with those in the experiment of Ref. Pinczuk,  for the filling factor range $`1/3\nu 1/5`$. For $`\mathrm{\Delta }^{44}(R)`$ the calculated values are approximately $`0.2`$ meV in this range (taking parameters of the experiment in Ref. Pinczuk, ), decreasing slightly with increasing magnetic field. The experimental values start at $`0.2`$ meV and decrease with magnetic field to $`0.11`$ meV at $`\nu =1/5`$. For $`\mathrm{\Delta }^{42}(R)`$ mode the theoretical value starts at approximately $`0.6`$ meV and increases with magnetic field to $`2.0`$ meV. The experimental value starts at $`0.45`$ meV and goes on increasing with magnetic field to $`0.7`$ meV. The theory is in qualitative agreement with experiment, both in regard to (i) the continuity of some modes across $`\nu =1/3`$ and appearance of new modes for $`\nu <1/3`$, and (ii) the observed trends as a function of the filling factor. In light of the rather small size of the numerical system, we find the level of quantitative agreement to be satisfactory.

Finally, we briefly discuss the calculational method. It involves many steps, which will only be outlined here for lack of space; more details can be found in the literatureManual ; Mandal ; PetersonJain ; MPJ . We have used the standard spherical geometryHaldane ; WuYang  which considers $`N`$ electrons on the surface of a sphere with a magnetic monopole of strength $`Q`$ placed at the center producing a radial magnetic field. The magnetic flux through the sphere is $`2Q\varphi _0`$ ($`\varphi _0=hc/e`$ is called the magnetic flux quantum); $`2Q`$ is an integer due to Dirac’s quantization condition. The single particle eigenstates are the monopole harmonicsWuYang  $`Y_{Qnm}(\mathrm{\Omega }_j)`$ where $`n=0,1,\mathrm{}`$ is the LL index, $`m=(Q+n),(Q+n)+1,\mathrm{},(Q+n)`$ labels the $`2(Q+n)+1`$ degenerate states in the $`n`$th LL, and $`\mathrm{\Omega }_j`$ represents the location on the sphere of electron $`j`$ with the usual coordinates. Since the degeneracy in successive LLs is not equal in this geometry, the ratio $`N/2Q`$ is not exactly equal to the filling factor for finite systems; the filling factor is defined as $`\nu =lim_N\mathrm{}N/2Q`$.

According to the CF theory, $`N`$ strongly interacting electrons at monopole strength $`Q`$ are transformed into $`N`$ weakly interacting composite fermions at monopole strength $`Q^{}=Qp(N1)`$. For a given $`Q`$ in the range $`1/3>\nu 1/5`$, we first consider $`Q^{}=QN+1`$, which lies in the range $`1>\nu 1/3`$. We diagonalize the electron problem at $`Q^{}`$ (allowing at most one electron in the second LL) to obtain eigenstates with definite total orbital angular momentum $`L`$. Denoting this basis by $`\{\mathrm{\Phi }_Q^{}\}`$, we construct a basis at $`Q`$ given by $`\{P_{LLL}\mathrm{\Phi }_1^2\mathrm{\Phi }_Q^{}\}`$, where $`\mathrm{\Phi }_1`$ is the wave function of one filled Landau level. $`P_{LLL}`$ denotes projection into the lowest electronic LL, which can be accomplished using the methods discussed previouslyManual . The basis at $`Q`$ obtained in this manner is not orthogonal. Following the techniques in Ref. Mandal, , we use the Gram-Schmidt procedure to obtain an orthogonal basis, and then diagonalize the interaction Hamiltonian in this basis. The various projections and Hamiltonian matrix elements require evaluation of multi-dimensional integrals, which is accomplished numerically by Monte Carlo. The uncertainty is determined by the standard deviation in the eigenvalues over different Monte Carlo runs.

A technical point in regard to our treatment of finite width should be noted. Following Ref. Xie, , we determine the interaction pseudopotentialsHaldane  in the infinite planar geometry, and feed them into exact calculations in the spherical geometry. For an infinite system, this of course gives the same result as doing the calculation entirely within the spherical geometry, but for finite systems the two methods are different. For zero thickness, the use of planar pseudopotentials, as opposed to the spherical pseudopotentials, makes only a slight deviation (within 6% for our calculations) in the excitation energies; we expect that the method works equally well for finite width systems as well.

We have only considered in this work filling factors where the ground state is incompressible. It would be of interest to investigate the filling factor dependence of the various excitations energies.

This work was supported in part by the National Science Foundation under grants no. DGE-9987589 (IGERT) and DMR-0240458. We are grateful to the High Performance Computing (HPC) Group led by V. Agarwala, J. Holmes, and J. Nucciarone, at the Penn State University ASET (Academic Services and Emerging Technologies) for assistance and computing time with the LION-XL cluster.


text:
# High Angular Resolution Measurement of Ion and Neutral Spectra as a Probe of the Magnetic Field Structure in DR21(OH)

## 1 Introduction

Measurements of magnetic fields in the environment of star forming regions are mostly carried out with observations of Zeeman splitting (Crutcher 1999), the polarimetry of dust continuum (Dotson et al. 2000; Lai et al. 2001, 2002), and the polarimetry of CO emission (Girart et al. 1999; Lai et al. 2003). The Zeeman observations only give the field strength along the line of sight, and the polarimetry of dust continuum and CO emission yield the field direction in the plane of sky. As a result, these observations do not determine the full strength and orientation of the magnetic fields in 3-D. To characterize the magnetic field environment, we need to measure directly the magnetic field strength in the plane of sky $`B_p`$. Without direct measurement of the full field strengths, the role magnetic fields play in star formation can only be assessed using indirect estimates of $`B_p`$ (the Chandrasekhar-Fermi method as used by Lai et al. 2001) and some statistical arguments (Crutcher 1999).

Houde et al. (2002) proposed an interesting new approach to probe the orientation of magnetic fields. Houde et al. (2000a,b) show that in weakly ionized plasma, even under weak fields ($`10\mu `$G), the cyclotron interaction between the ions and the field can reduce the linewidth and suppress the high velocity wings in the ion spectra when compared to coexistent neutral spectra. Their HCN and HCO+ observations with the Caltech Submillimeter Observatory (CSO) show evidence of these predicted features. Because the ion-to-neutral linewidth ratio depends largely on the inclination angle of the magnetic fields (no cyclotron interaction along the field lines), the full strength and orientation of magnetic fields can be obtained by combining the observations of Zeeman splitting, polarimetry of dust continuum and/or CO emission, and the ion-to-neutral linewidth ratio. Houde et al. (2002) have constructed the 3-D field structure of M17 using this method.

However, it is difficult to interpret the differences in the velocity dispersions of the ion/neutral spectra unambiguously as being caused by the cyclotron interaction, if the beam is large enough to include the spectral line emissions from physically and chemically different components. The intrinsic difference between the velocity structure of the ions and neutrals will also contribute to the ion-to-neutral linewidth ratio. For example, if a velocity gradient exists in a large region with HCN, but HCO+ only coexists with HCN in a small part of this region, the HCN linewidth is naturally larger than the HCO+ linewidth. Therefore, it is important to obtain high spatial resolution spectral line maps to resolve the regions with widely different physical and chemical characteristics.

In this paper, we present the comparison between H13CN and H13CO+ spectra in DR21(OH) observed with the millimeter array of the Owens Valley Radio Observatory (OVRO). DR21(OH) is an excellent object to study the magnetic field structure, because it is in the early stage of massive star formation. DR21(OH) is one of the few sources with measured average line-of-sight magnetic field strength from CN Zeeman detections (Crutcher et al. 1999) and field directions in the plane of sky mapped with BIMA using both dust and CO polarimetry (Lai et al. 2003). Deriving the full strength and 3-D structure requires knowledge of the inclination angle of the magnetic field to the line of sight, which can be obtained by measuring the ion-to-neutral linewidth ratio.

## 2 Observations and Results

The observations were made with the OVRO Millimeter Array at 3mm in the Compact, Low, and High resolution array configurations between October 2002 and May 2003. We obtained simultaneously spectral line images of HCN & HCO+ and H13CN & H13CO+; 13CS and N2H+ were observed in the opposite sidebands. The synthesized beam is 6.4′′$`\times `$6.9′′ for natural weighting and spectral resolutions are 0.434 km s-1. Among these images, the HCN and HCO+ spectra show features of optically thick emission with strong self-absorption; N2H+ has a very different spatial distribution compared to other species; H13CN, H13CO+, and 13CS have a similar spatial distribution, but 13CS has a different velocity structure from the other two. Therefore, H13CN and H13CO+ are the only neutral/ion pair that can be used to examine the cyclotron effect in our data.

Figure 1 shows the 3mm continuum and the integrated flux of H13CN and H13CO+. The two compact cores in DR21(OH), MM1 and MM2, are resolved. H13CN and H13CO+ both trace the dense cores very well. H13CO+ is slightly more extended than H13CN, which is probably due to the lower critical density of H13CO+. Therefore, H13CN and H13CO+ seem to occupy the same volume. Figure 2 shows the normalized H13CN and H13CO+ spectra overlaid on their flux ratio map. Note that H13CN has three hyperfine lines with intensity ratios of 1:5:3 for the optically thin case. The spectra show that DR21(OH) contains two dominant velocity components – one at $`5`$ km s-1 and the other at $`1`$ km s-1, which are likely to be associated with MM1 and MM2, respectively. At the positions between MM1 and MM2, the spectra are more complex and may contain an absorption feature at $``$2.5 km s-1. Overall, the spectra of H13CN and H13CO+ trace each other well, which is strong evidence for the coexistence of H13CN and H13CO+, as also suggested by their spatial intensity correlation.

A difference in the ion/neutral linewidth can arise in some cases from chemical inhomogeniety, especially, if multiple velocity features are present. In Figure 2 we show the distribution of the H13CN-to-H13CO+ flux ratio (in greyscale) which is a measure of the chemical differentiation in the core. We can safely discard the effects of chemical differentiation by excluding regions where the flux ratio changes appreciably over the small spatial scale sizes of the OVRO beam. However, chemical differentiation alone will not always cause linewidth differences. The velocity structure of the region also needs to be considered. For example, if the velocity structure is uniform across the cores such that the shape and the central velocity of the spectra stay the same, any local variation in the ion-to-neutral flux ratio will result in a weighting factor that affects only the line intensities, not the ion-to-neutral linewidth ratios for the whole core. However, in DR21(OH) at least two velocity components are detected and therefore, it is critical to isolate the regions with different dynamics when interpreting the ion-to-neutral linewidth ratio as the cyclotron effect. We see a large variation in their flux ratio (by factors of 0.5 to 4) over the entire region, but it appears to be smooth except near the southwest parts of the core. In the eastern parts of the core the variation in the flux ratio is small, and the fully resolved velocity structure is dominated by a single feature. Thus, in this region of DR21(OH), we are able to minimize any contribution to the ion/neutral linewidth differences arising from chemical inhomogeneity. In our analysis we do not include the southwest region where rapid and significant variation of ion-to-neutral flux ratio is observed. Furthermore, the ion/neutral spectra in the southwest region have low S/N.

In the eastern part of MM1 where the $``$5 km s-1 component dominates and the spectral line shapes are better defined, we can see that H13CO+ is slightly narrower than H13CN. To eliminate the complication of the H13CN hyperfine lines, we simulate H13CN spectra with all three hyperfine lines using the H13CO+ as a ‘template’. For the case of a single velocity component, if the linewidth of H13CN is smaller than half of the velocity separations between hyperfine lines, the difference spectrum between the observed H13CN and that simulated from H13CO+ is largely free from the confusion caused by the hyperfine components and thus is a useful measure of the cyclotron interaction. If indeed H13CN has larger linewidths than H13CO+, the difference spectra will have two positive peaks near the wings of the line profile as demonstrated in Figure 2d. As the line shapes of H13CN and H13CO+ are approximately Gaussians, the integration of the difference spectra over the extent of the main hyperfine component gives a measure of the difference of the dispersions in the spectra.

$$(I_nI_i)𝑑v=\sqrt{2\pi }(\sigma _n\sigma _i),$$
(1)
where $`I_n`$ and $`I_i`$ are the normalized H13CN and H13CO+ spectra and $`\sigma _n`$ and $`\sigma _i`$ are the dispersions of H13CN and H13CO+, respectively. Figure 3a shows the map of $`(\sigma _n\sigma _i)`$ which is calculated by integrating the difference spectra over the $`5`$ km s-1 component ($`v_{LSR}=9`$ to $`2.5`$ km s-1). These $`(\sigma _n\sigma _i)`$ estimated using Eq(1) are in good agreement (within 10%) with those obtained by fitting Gaussians to the observed H13CN and H13CO+ line profiles at (+10′′, +5′′). The signal-to-noise ratio (S/N) of $`(\sigma _n\sigma _i)`$ is presented in Figure 3b. We use $`(\sigma _n\sigma _i)`$ and $`\sigma _i`$ measured from H13CO+ spectra to derive $`(\sigma _i/\sigma _n)`$ in Figure 3c. It is evident from the results shown in Figure 3 that in the eastern side of MM1 the ion linewidth is smaller than the neutral linewidth. We do not present similar analysis for the $`1`$ km s-1 component, because of the low S/N of $`(\sigma _n\sigma _i)`$.

## 3 Discussion

The interferometer observations with OVRO have the advantage to resolve the velocity and chemical structures in the DR21(OH) core with high spatial resolution. In the eastern parts of DR21(OH) our observations show that H13CN and H13CO+ coexist over 15′′$`\times `$10′′ in a single velocity feature and confirm that the ion line is indeed narrower than the neutral line. We use the ion-to-neutral linewidth ratio $`(\sigma _i/\sigma _n)`$ shown in Figure 3c to infer the inclination angle of the magnetic field to the line of sight ($`\alpha `$) according to Equation (11) of Houde et al. (2002). To derive the value of ($`\alpha `$) using the ion/neutral linwidth ratio we require knowledge of the gas motions within the core. Houde et al. consider the cases with neutral flow collimation widths, $`\mathrm{\Delta }\theta `$, ($`\mathrm{\Delta }\theta =90`$ means no collimation as in a purely turbulent case). For a given value of $`(\sigma _i/\sigma _n)`$ a range of values for $`\alpha `$ are possible depending on the values for flow collimation widths (see Figure 4 in Houde et al.). For example the average ion-to-neutral linewidth ratio, $`(\sigma _i/\sigma _n)`$ of 0.82 will correspond to an inclination angle $`\alpha `$ = 76, 60and 36 for $`\mathrm{\Delta }\theta `$=20, 45and 90 respectively. The dependence of the estimates for $`\alpha `$ on $`\mathrm{\Delta }\theta `$ is less critical for $`\mathrm{\Delta }\theta `$ 90, a case when turbluent motions dominate over collimated flows. Although MM1 in DR21(OH) does power high velocity CO outflows (Lai et al. 2003), there is no evidence that the spectra of ions and neutrals observed at the rest velocity of the cores are contaminated by collimated flows. Therefore, to derive the magnetic field inclination, we have ignored collimated flows and we assume the gas is fully turbulent ($`\mathrm{\Delta }\theta =90`$). Figure 4 shows the inclination angle $`\alpha `$ ranges from 30 in the northern part of MM1 to 42 in the south, and the uncertainty in $`\alpha `$ is $``$ 4–8. The variation in $`\alpha `$ is small (12) and is comparable to the variation of the overlaid $`B_p`$ vectors derived from the dust polarization (black) and CO polarization (white). Therefore, the magnetic field directions seem to be approximately uniform across the MM1 region. The average line-of-sight magnetic field strength in MM1 $``$ $``$0.36$`\pm `$0.10 mG as measured by CN Zeeman observations (Crutcher et al. 1999). With the average $`\alpha =36`$, the field strength in the plane of sky ($`B_p`$) is 0.25$`\pm `$0.07 mG, and the full magnetic field strength is $`0.44\pm `$0.12 mG. The direction of the field is out of the plane of the sky towards the observer at an angle of 36to the line of sight, and its projection on the sky is at a position angle $`75`$.

Lai et al. (2003) estimate $`B_p`$0.9 mG for MM1 from the observed dispersion in the CO polarization position angles using the Chandrasekhar-Fermi formula. This value is a factor of 3 larger than what is derived from the cyclotron interaction data presented here. This discrepancy could be partly due to the low values for the dispersion caused by the smoothing of polarization by the finite beam size of the BIMA observations. Also, the ‘calibration factor’ of the Chandrasekhar-Fermi formula (0.5) used by Lai et al. (2003) is an average value of clouds in numerical simulations (Ostriker, Gammie, & Stone 2001), which may not be a proper value DR21(OH).

With the knowledge of the full magnetic field strength, we can evaluate the role that magnetic fields play for star formation in MM1. First, we can calculate the mass-to-magnetic flux ratio from $`M/\mathrm{\Phi }_B=1.0\times 10^{20}N(H_2)/|B|`$ in units of the critical mass-to-flux ratio (Crutcher 1999), where $`N(H_2)`$ is the column density in cm-2 and $`|B|`$ is the field strength in $`\mu `$G. With $`N(H_2)=2\times 10^{23}`$ cm-2 and B=0.44 mG, $`M/\mathrm{\Phi }_B=`$4.5. Therefore, MM1 is clearly magnetically supercritical. However, turbulence seems to provide dominant support. The H13CO+ linewidth ($``$ 2.2 km s-1) is much larger than its thermal linewidth for $`T_k`$=50 K (0.3 km s-1), indicating strong turbulence in DR21(OH). The relative importance of magnetic, thermal, and turbulent support can be evaluated by comparing their energy to the gravitational energy. We find that in MM1 the magnetic-to-gravitational energy ratio is 0.05, the thermal-to-gravitational energy ratio is 0.16, and the turbulent-to-gravitational energy ratio is 0.66. Therefore, turbulence provides the dominant support in MM1. Although our conclusion on the status of MM1 is similar to what has been suggested by Crutcher et al. (1999) with statistical argument for the full magnetic field strength, our results exclude the possiblity of the magnetic field dominance which may be true only if the field orientation in 3-D is very close to the plane of the sky.

We thank the referee for helpful suggestions. This research was conducted at the Jet Propulsion Laboratory, California Institute of Technology, with support from the National Aeronautics and Space Administration, while SPL held a National Research Council Resident Research Associateship at JPL. The OVRO Millimeter Array is supported by NSF grant AST-99-81546.


text:
# Multiepoch Sky Surveys and the Lifetime of Quasars

## 1. Introduction

The quasar lifetime ($`t_Q`$) is an important timescale for supermassive black hole growth and determination of the mechanism(s) that fuel Active Galactic Nuclei (AGN). The net time that black holes are radiating at quasar luminosities, and therefore accreting at approximately the Eddington rate, sets the importance of luminous accretion for the production of the present-day space density of dormant, supermassive black holes. If this lifetime is longer than $`t_Q>10^7`$ years, then a quasar phase would have led to the accretion of much of the mass in the largest present-day supermassive black holes. However, the low space density of quasars, even at high-redshift, then implies that a quasar phase was relatively rare and not all supermassive black holes were quasar hosts. In contrast, shorter quasar lifetimes require more of the present-day supermassive black hole population to have been quasars, but then to have accreted a smaller fraction of their total black hole mass as quasars.

Most estimates of the net quasar lifetime produce values in the range $`t_Q=10^610^8`$ years, or lower limits that are consistent with these values (for a recent review see Martini, 2003). These estimates are generally based on either integral or demographic arguments, which combine data on the present-day population of supermassive black holes and the accretion by the high-redshift quasar population (e.g. Sołtan, 1982; Haehnelt et al., 1998; Yu & Tremaine, 2002), or incorporate quasars directly into models for galaxy evolution (Kauffmann & Haehnelt, 2000). However, most demographic models only constrain the net time that a supermassive black hole is radiating above the luminosity threshold for quasars and do not address the possibility that a supermassive black hole goes through multiple, episodic quasar phases that sum to produce the net lifetime. The episodic lifetime is an important parameter for the accretion physics, as it is the timescale that an accretion disk can maintain accretion at approximately the Eddington rate. Existing constraints on the episodic lifetime suggest $`t_Q>10^410^5`$ and are primarily based on the radiative properties of quasars. These lower limits are due to the minimum lifetime required to explain the proximity effect in the Lyman$`\alpha `$ forest (Bajtlik et al., 1988; Scott et al., 2000) and the sizes of ionization-bounded narrow-line regions (Bennert et al., 2002). A recent observation of the transverse proximity effect in He II does imply a longer lifetime of $`t_Q10^7`$ years (Jakobsen et al., 2003); however there is currently only one such observation.

Measurement of both the episodic and net lifetimes of quasars are important to determine the methods that trigger quasar activity, as well as the physics of the accretion process. If net and episodic lifetimes are similar, and by implication long, quasars could be produced by the relatively rare merger of approximately equal-mass, gas-rich galaxies. This is currently the most favored mechanism invoked to explain the mass accretion rates necessary to fuel quasars, as well as their rapid decrease in space density at low redshifts (Carlberg, 1990; Barnes & Hernquist, 1992). In contrast, if the episodic lifetime of quasars is significantly less than the net lifetime, then quasar activity must be triggered by a phenomenon that is more frequent than large galaxy mergers. Some fraction of quasars could also be powered by the tidal disruption of individual stars as they pass near peribarathron (Hills, 1975; Young, Sheilds, & Wheeler, 1977), which could produce quasar luminosities for on order one year (Ulmer, 1999). The current estimates and limits on both the net and episodic lifetimes are still uncertain by several orders of magnitude due to the quality with which the current data are fit by models, and the assumptions of the models themselves. In the present work we describe a method that could improve the limit on the episodic lifetime by approximately two orders of magnitude and determine if the episodic and net lifetimes are comparable.

The most direct method to determine $`t_Q`$ is to actually measure a quasar turn on and off. This is impractical for any given quasar, given the long lifetime implied by the methods described above, but is possible for observations of an ensemble of quasars, such as are now being produced by the 2dF (Croom et al., 2002) and SDSS (York et al., 2000) surveys. The 2dF survey has produced a final catalog of 23,424 quasars, while the SDSS survey plans to complete a catalog of 100,000 quasars over the next several years. In the next section we provide a detailed calculation of how well $`t_Q`$ can be computed from these spectroscopically-identified quasar samples when compared with photometry from at least one additional epoch. We also describe how the SDSS quasar catalog could be compared to the existing Palomar Observatory Sky Survey (POSS), particularly in light of the lower-quality and comparable sensitivity of these plates. This method is illustrated through application to the SDSS Early Data Release (EDR) quasar catalog. We conclude with a discussion of the implications of these measurements or upper limits on $`t_Q`$ in the final section.

## 2. Method

We shall first examine the problem of quasar lifetimes with a simple model for the quasar luminosity: the quasar is either “on” or “off” with a characteristic on lifetime $`t_Q`$. In this scenario, the probability that any quasar will turn off over a time baseline $`\mathrm{\Delta }t`$ is $`P_{off}=\mathrm{\Delta }t(z)/t_Q`$, where $`\mathrm{\Delta }t(z)=\mathrm{\Delta }t/(1+z)`$ is the time baseline in the quasar frame from the spectroscopic identification to the second, photometric epoch. For a sample of $`N`$ quasars where only $`N_{on}`$ are observed to be on at the second epoch, the lifetime can then be directly calculated as

$$t_Q=\frac{\mathrm{\Sigma }\mathrm{\Delta }t(z_i)}{\mathrm{\Sigma }iN_{\mathrm{on}}}$$
(1)
where this assumes that a quasar will definitely be detected if it is on at the second epoch. Here we adopt the same definition of a quasar as the SDSS EDR quasar catalog: an AGN brighter than $`M_i^{}=23`$ mag (Schneider et al., 2002). For a spectroscopically-identified quasar with one additional, photometric epoch we consider a quasar as on if it is above this luminosity threshold, and off if it is fainter than this limit. The limitations of this assumption are discussed in §3.

One survey that could provide the second epoch photometric observation for the SDSS and 2dF quasar samples is the proposed Large Synoptic Survey Telescope (LSST; Tyson et al., 2002). The LSST is intended to regularly survey the entire visible sky to a depth many magnitudes fainter than the limit of the SDSS and 2dF quasar samples. This telescope will therefore easily obtain high-quality photometry of all of the quasars in the final 2dF and SDSS catalogs provided they are still on at the time of the LSST observations. We have modeled the quality with which the SDSS and 2dF samples could be used to measure the quasar lifetime by generating mock catalogs with a range of lifetimes and the same magnitude and redshift distributions as the SDSS EDR. Figure 1 shows the quasar lifetimes and uncertainty in $`t_Q`$ attainable with this second epoch LSST observation. If the LSST obtains a second baseline measurement of these 120,000 quasars ten years after their initial spectroscopic identification, then this will measure $`t_Q`$ to within a factor of two if $`t_Q<10^5`$ years, or otherwise set a lower limit to the quasar lifetime if no quasars have been observed to turn off. The upper and lower limits to the quasar lifetime are set by the length of the time baseline and the probability that no quasars will actually be observed to turn off over this baseline, respectively. For long values of the lifetime only a small number of quasars will be observed to have turned off; we have calculated the one-sided, $`1\sigma `$ confidence limits on the number of quasars that are off using Poisson statistics as described in Gehrels (1986).

For second-epoch survey material that does not extend many magnitudes below the limits of the spectroscopic sample, the sensitivity limit, quasar variability, and the quality of the photometric calibration could all produce a false signal of off quasars. These additional effects can be accounted for by assigning a weight $`w_i`$ to each quasar that quantifies the probability that it could be detected at the second epoch, given an assumed magnitude distribution based on the spectroscopic epoch, a parameterization of quasar variability, and the sensitivity of the second epoch photometric survey:

$$w_i=𝑑m_{ij}P_m(m_{ij}|\mu _i)P_d(m_{ij})$$
(2)
The expected magnitude distribution for the quasar $`P_m(m_{ij}|\mu _i)`$ can be calculated from a Gaussian distribution centered on its assumed mean magnitude $`\mu _i`$ at the time of identification, as well as the photometric uncertainty and quasar variability. Due to a variant of the Eddington bias, this mean magnitude will on average be fainter than the magnitude $`m_i`$ at the time of spectroscopic selection because of quasar variability $`\sigma _V`$, the photometric uncertainty $`\sigma _P`$, and the slope of the quasar number-magnitude relation. The variability of quasars is commonly parametrized with a structure function (Simonetti et al., 1985; Hughes et al., 1992) and long term monitoring of quasars shows that the amplitude of the variability increases with time (Hawkins, 2002). The width of the distribution of likely quasar magnitudes at the second epoch is thus the quadrature sum of the photometric uncertainty and the time-dependent variability: $`\sigma (\mathrm{\Delta }t)=\sqrt{\sigma _P^2+\sigma _V^2(\mathrm{\Delta }t)}`$. This magnitude distribution is then convolved with the detection probability $`P_d(m_{ij})`$ for a quasar with magnitude $`m_{ij}`$ at epoch $`j`$. With this weighting function the total number of expected quasars in the second epoch becomes the sum of these weights $`\mathrm{\Sigma }w_i`$ and the total time baseline available becomes $`\mathrm{\Sigma }w_i\mathrm{\Delta }t(z_i)`$. The equation for the lifetime is rewritten:

$$t_Q=\frac{\mathrm{\Sigma }w_i\mathrm{\Delta }t(z_i)}{\mathrm{\Sigma }w_iN_{\mathrm{on}}}.$$
(3)
A particularly useful set of archival observations are the POSS-I plates obtained in the 1950s. All of these observations have been digitized, are publicly available on CDROM, and provide the longest time baseline for the SDSS quasar catalog. As above for the LSST calculation, we have used the redshift and magnitude distribution of the SDSS EDR quasars to model the sensitivity of the SDSS $`+`$ POSS-I to the quasar lifetime. For regions which were not covered by the POSS-I survey, we have used the more recent POSS-II data from the 1980s. While the 23,000 2dF quasars will also be valuable when combined with future LSST observations, they are not suitable for comparison with the archival plate material as the spectroscopic quasar candidates were identified from these plates.

The first step in this analysis was the photometric calibration of the POSS plates, which was carried out with USNO stars selected from the same plates as the EDR quasars, to have similar colors (Richards et al., 2002), and small proper motions. This produced a catalog of $`9800`$ stars that were used to calibrate the POSS plates directly onto the SDSS magnitude system. The rms uncertainty in the photometric calibration is $`0.4`$ mag. These same stars were then used to determine the detection efficiency as a function of magnitude by measuring the number of stars identified with SExtractor (Bertin & Arnouts, 1996) in $`0.5`$ mag bins. The USNO catalog is complete to $`V=21`$, which is fainter than the POSS-I limit and the magnitude for SDSS quasar selection. The EDR quasars were detected using the same SExtractor parameters and requiring that the source on the POSS-I plate be within $`4^{\prime \prime }`$ of the SDSS coordinates. Finally, quasar variability was parameterized with the structure function calculated by Hawkins (2002). This structure function, which is similar to that recently derived from the SDSS EDR quasars by de Vries, Becker, & White (2003), was used to set the width of the expected quasar variability distribution as a function of time baseline in the quasar frame. Figure 2 shows the quasar lifetimes and uncertainty in $`t_Q`$ attainable with a combination of the 100,000 SDSS quasars and the POSS epoch. These data will thus in principle be able to measure the quasar lifetime if $`t_Q<10^{5.5}`$, although in practice the sensitivity will be somewhat less due to uncertainties in quasar variability and the photometric calibration of the plates. Application of this method to the SDSS EDR sample sets a lower limit of $`t_Q>20,000`$ years on the episodic quasar lifetime.

## 3. Discussion

The present lower limit of $`t_Q>20,000`$ years from the SDSS EDR quasar sample is in agreement with the lower limits to $`t_Q`$ from the radiative arguments outlined in the Introduction. This limit can be improved by approximately the factor of 30 increase in sample size expected from the final SDSS quasar sample. The potential measurement or limit of $`t_Q>10^{5.5}`$ years from the full sample will either measure $`t_Q`$ or produce an improved lower limit to the episodic lifetime of quasars. Future observations of the full quasar sample with the LSST and a ten-year baseline could push the lower limit to $`t_Q10^6`$ years through use of the combined baseline provided by POSS and LSST data. These observations will also be able to set interesting limits on the rate at which tidal disruptions of stars produce quasar luminosities.

The measurement of $`t_Q`$ via LSST observations of the complete SDSS and 2dF quasar samples will provide a comparable constraint on $`t_Q`$ to many of the existing, indirect estimates from demographic arguments, which typically find values of $`t_Q10^7`$ years (Haehnelt et al., 1998; Kauffmann & Haehnelt, 2000; Yu & Tremaine, 2002). If the sky-survey measurement or lower limit and these other estimates of $`t_Q`$ remain comparable in the face of improved data, this will be a strong validation of the assumptions of radiative efficiency, black hole mass function, and Eddington fraction that are key components of the demographic models.

Observations with the LSST will also provide superb measurements of quasar variability, which motivates a more careful consideration of our definition of the episodic quasar lifetime. A particularly important question is what is meant by a quasar that is off, relative to a quasar that is instead experiencing extreme variability. Here we have defined a quasar as off if it becomes fainter than the absolute magnitude $`M=23`$ used to define the quasar sample. This definition is motivated by the small number of available epochs for the quasars in the POSS data and the poor photometric precision of the plates, which make it difficult to characterize the variability of individual quasars. Our definition, however, allows for the spurious identification of off quasars due to misclassified, extremely variable objects, or simply variable quasars near the luminosity limit that defines the quasar class. In principle, a quasar should not be classified as off if it has temporarily faded to become a lower-luminosity AGN and then subsequently increases in brightness to become a quasar again. For quasars with only a small number of epochs, particularly of poor photometric precision, this point will be hard to circumvent.

The multiple-epoch and greater precision observations with the LSST, however will allow this point to be addressed in two ways. First, extremely variable quasars can be identified and not included in the sample through a variability cut for the lifetime estimate. A second approach is to include even these extremely variable quasars, yet require that they fade several magnitudes below the defining luminosity of a quasar for some number of years before they are formally classified as off. This later approach is more appealing as eliminating the most variable segment of the quasar population may eliminate quasars that are near the end of their lives. If the quasar has turned off, whether completely or is simply temporarily quiescent, then this will offer an unprecedented opportunity to study the quasar host galaxy.

The luminosity and variability evolution of quasars can be considered in greater detail with the multiple-epoch LSST observations. We have adopted the simple assumption that quasars undergo no luminosity evolution while they are on, other than their known variability about some constant mean. Many models for quasars, however, adopt an exponentially decaying luminosity (e.g. Haehnelt et al., 1998). If the mean luminosity of a quasar decays exponentially over an $`e`$folding timescale $`t_{Q}^{}{}_{}{}^{}`$, rather than experiencing a simple shut off after $`t_Q`$, then the SDSS and LSST photometry could be used to search for a mean decrease in the luminosity of the quasar population, provided it is relatively short. For an assumed photometric uncertainty of $`0.02`$ mag and a baseline of ten years, a dataset of 100,000 quasars could in principle measure a change in mean luminosity for a characteristic decay time of close to $`t_{Q}^{}{}_{}{}^{}10^5`$ years, although in practice quasar variability will limit the sensitivity of such a measurement to much shorter characteristic decay timescales.

Longer lifetimes, whether in the simple on–off model or with more complex luminosity evolution, could naturally be even better investigated with yet larger samples. While a significantly larger spectroscopically-identified sample is unlikely in the near future, a larger color-selected sample is possible. Based on studies with the color-selection algorithms used for the SDSS quasar sample, there are regions of color space that are almost exclusively occupied by quasars (Richards et al., 2002). Thus large numbers of quasars could be identified through color-selection alone, either in the SDSS or with LSST, and then monitored over a number of years. While the remaining contaminating population must be considered, these are almost exclusively non-variable objects such as white dwarfs, and many of these contaminants could be identified via their proper motions.

Completion of the SDSS should also produce improved estimates of the net lifetime of quasars from other techniques such as quasar clustering (Martini & Weinberg, 2001; Haiman & Hui, 2001). A preliminary comparison of the 2dF clustering measurements to the models in Martini & Weinberg (2001) suggest the quasar lifetime may be as short as $`t_Q10^6`$ years. If refined estimates based on SDSS measurements continue to point toward a shorter lifetime, then the net lifetimes determined from clustering will be comparable to the range of lifetimes probed by the method outlined here and these two methods could determine if the net and episodic lifetimes are comparable. If they are, then supermassive black holes only experience one quasar phase and the mechanism that ignites and fuels quasar activity is a relatively rare event. In contrast, if the episodic lifetime is much shorter than the net lifetime, then quasar activity must be triggered by additional mechanisms than just the relatively rare mergers of gas-rich galaxies.

We acknowledge helpful comments from Pat Hall, Gordon Richards, Martin Gaskell, and the referee, Elaine Sadler. PM was supported by a Carnegie Starr Fellowship and DPS was partially supported by National Science Foundation Grant NSF-03007582. Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The National Geographic Society - Palomar Observatory Sky Atlas (POSS-I) was made by the California Institute of Technology with grants from the National Geographic Society. The Second Palomar Observatory Sky Survey (POSS-II) was made by the California Institute of Technology with funds from the National Science Foundation, the National Geographic Society, the Sloan Foundation, the Samuel Oschin Foundation, and the Eastman Kodak Corporation.


text:
# Incident-Energy Dependence of the Effective Temperature in Heavy-Ion Collisions

## I Introduction

Experimental data on transversal momentum distributions of kaons produced in central Pb+Pb or Au+Au collisions show an anomalous dependence on the collision energy ags ; sps ; rhic . The effective temperature $`T^{}`$, or slope parameter of the transversal-momentum spectra note1 , increases with energy in AGS and RHIC energy domains. However, in the SPS energies the effective temperature keeps approximately constant. Recently, it was argued MG2a  that this behaviour might be caused by a modification of the equation of state in the transition region between confined and deconfined matter, as suggested by Van Hove Hov  a long time ago, and could be considered as new signal of deconfinement in the SPS energy domain.

Our main object in the present work is to study the behaviour of the effective temperature for K+ in several energy domains. For this purpose, we apply the recently developed SPheRIO Aguiar2001 ; Aguiar2002  code for hydrodynamics in 3+1 dimensions, using both Landau-type compact initial conditions and spatially more spread ones. For the latter we used the average over many events generated by NeXus nexus1 ; nexus2 . We show that initial conditions given in small volume, like Landau-type ones, are unable to reproduce the effective temperature together with other data (multiplicities and rapidity distributions). These quantities can be reproduced altogether only when using a large initial volume with an appropriate velocity distribution.

Besides, it seems that the increase of $`T^{}`$ in the RHIC energy domains is caused mainly by the larger expansion time in the hadronic phase for higher incidente energy, which implies lower freeze-out temperature as $`\sqrt{s}`$ increases HN ; NNOP . It seems that, within our analyses with NeXus initial conditions, the RHIC incident energies are not enough to produce noticeable increase in the transverse acceleration during the quark-gluon plasma (QGP) phase.

## II Hydrodynamics in 3+1 dimensions

The hydrodynamical code in 3+1 dimensions used here is based on the technique called *Smoothed Particle Hydrodynamics* (SPH) Lucy1977 ; Gingold1977 . This is a method which uses Lagrangian coordinates. The fluid is represented by small volumes called SPH-particles and the equations for the fluid evolution become a system of ordinary differencial equations for the SPH-particles, in this representation.

The numerical code SPheRIO note2  is a suitable implementation of this method for the relativistic nuclear collisions. An advantage of this technique is the convenience in the study of problems where the geometry is highly irregular, as is the case of non-central relativistic collisions. A SPH-particle has attached to it conserved quantities; in the present version of SPheRIO, the entropy and the baryonic number are the quantities which are kept constant during the fluid evolution. In the SPH representation, the entropy density and the baryonic density are parametrized as:

$`s(\tau ,𝒙)`$ $`=`$ $`{\displaystyle \frac{1}{\tau \gamma }}{\displaystyle \underset{i}{}}\nu _i^sW(𝒙𝒙_i;h),`$ (1)
$`n_b(\tau ,𝒙)`$ $`=`$ $`{\displaystyle \frac{1}{\tau \gamma }}{\displaystyle \underset{i}{}}\nu _i^bW(𝒙𝒙_i;h),`$ (2)
where $`\nu _i^s`$ ($`\nu _i^b`$) is the entropy (baryonic number) of the $`i`$-th SPH-particle, $`\gamma =1/\sqrt{(1v_x^2v_y^2\tau ^2v_\eta ^2)}`$, and $`W`$ is the interpolating kernel with width $`h`$. Here we use the hyperbolic coordinates $`\tau \sqrt{(t^2z^2)}`$, $`x`$, $`y`$ and $`\eta 1/2\mathrm{ln}[(t+z)/(tz)]`$ which are convenient for a system in rapid longitudinal expansion. The equations for the SPH-particles are given by:

$`{\displaystyle \frac{d𝒙_i}{d\tau }}`$ $`=`$ $`𝒗_i,`$ (3)
$`{\displaystyle \frac{d}{d\tau }}\left(\gamma \nu ^s{\displaystyle \frac{P+\epsilon }{s}}𝝎\right)_i`$ $`=`$ $`{\displaystyle \underset{j}{}}\nu _i^s\nu _j^s\left[{\displaystyle \frac{P_i}{\tau \gamma _i^2s_i^2}}+{\displaystyle \frac{P_j}{\tau \gamma _j^2s_j^2}}\right]\mathbf{}_iW(𝒙_i𝒙_j;h),`$ (4)
where $`𝝎=(v_x,v_y,\tau ^2v_\eta )`$. The pressure, $`P`$, and the energy density, $`\epsilon `$, are related to the entropy density $`s`$ and the baryonic density $`n_b`$ via equation of state.

### II.1 Decoupling criterion

The distribution of final particles is obtained by using Cooper-Frye’s prescription CF :

$$E\frac{d^3N}{dp^3}=_{\mathrm{\Sigma }_{fo}}𝑑\sigma pf(pu),$$
(5)
where $`f(pu)`$ is the distribution function (Bose-Einstein or Fermi-Dirac) and the integral is evaluated on the freeze-out surface $`\mathrm{\Sigma }_{fo}`$. In the SPH representation, Eq. (5) is written as

$$E\frac{d^3N}{dp^3}=\underset{i}{}\frac{\nu _i^s\widehat{n}_ip}{s_i|\widehat{n}_iu_i|}f(pu_i),$$
(6)
where the sum is over all SPH-particles on the freeze-out surface and $`\widehat{n}`$ is the normal four-vector to this surface, given by $`\widehat{n}_\mu (T/\tau ,T/x,T/y,T/\eta )`$.

## III Equations of state

Let us consider here two sets of equations of state. The first one type-I (EOS-I) has a first-order phase transition between an ideal gas of massless quarks (u, d and s) and gluons and an ideal hadron gas (baryon number is assumed to be zero), where the pressure, the energy density and the entropy density are given by:

$`P_q`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{90}}g_qT^4B,`$ (7)
$`\epsilon _q`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{30}}g_qT^4+B,`$ (8)
$`s_q`$ $`=`$ $`{\displaystyle \frac{2\pi ^2}{45}}g_qT^3,`$ (9)
$`P_h`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{90}}g_hT^4,`$ (10)
$`\epsilon _h`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{30}}g_hT^4,`$ (11)
$`s_h`$ $`=`$ $`{\displaystyle \frac{2\pi ^2}{45}}g_hT^3.`$ (12)
Here $`g_h=16`$ is an effective parameter (see MG2b  for details), $`g_q=95/2`$, and $`B`$ is the bag model parameter.

The second one, type-II (EOS-II), somewhat more realistic than the previous one, considers a first-order phase transition between a QGP and a hadronic resonance gas (baryon number is taking into account). In the QGP we consider an ideal gas of massless quarks (u, d, s) and gluons. The thermodynamical quantities are given by

$`P_q`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{90}}g_qT^4+{\displaystyle \frac{3}{2}}\left(T^2\mu _q^2+{\displaystyle \frac{\mu _q^4}{2\pi ^2}}\right)B,`$ (13)
$`\epsilon _q`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{30}}g_qT^4+{\displaystyle \frac{9}{2}}\left(T^2\mu _q^2+{\displaystyle \frac{\mu _q^4}{2\pi ^2}}\right)+B,`$ (14)
$`s_q`$ $`=`$ $`{\displaystyle \frac{2\pi ^2}{45}}g_qT^3+3T\mu _q^2,`$ (15)
$`n_q`$ $`=`$ $`3\left(T^2\mu _q^2+{\displaystyle \frac{\mu _q^3}{\pi ^2}}\right).`$ (16)
The hadronic phase is composed of resonances with mass below 2.5 GeV/c2, where volume correction is taken into account. The corrected-volume pressure, $`P^{ex}`$, is written in function of the pressure for an ideal gas $`P^{id}`$:

$$\underset{i}{}P_i^{ex}(T,\mu _i)=\underset{i}{}P_i^{id}(T,\stackrel{~}{\mu }_i),$$
(17)
where $`\stackrel{~}{\mu }_i=\mu _iV_i_jP_j^{id}(T,\stackrel{~}{\mu }_j)`$, ($`V_i`$: volume of the hadron $`i`$). Others thermodynamical quantities are given following the relation:

$$\chi _i^{ex}(T,\mu _i)=\frac{\chi _i^{id}(T,\stackrel{~}{\mu }_i)}{1+_jV_jn_j^{id}(T,\stackrel{~}{\mu }_j)},$$
(18)
where $`\chi `$ represents $`n`$, $`\epsilon `$ or $`s`$.

In both equations of state, the transition temperature is assumed to be 160 MeV.

## IV Initial conditions and results

In order to study the behaviour of the effective temperature, we began with Landau-type initial conditions, where the fluid is at rest at $`t=0`$ fm/c, and the matter is localized in a Lorentz contracted sphere. The initial energy density is assumed to be constant and given by: $`\epsilon _0=E/V`$, where $`E=\zeta (\sqrt{s}2m)A`$ is the available kinetical energy in the collision, and $`V=\gamma _0^1V_0=\gamma _0^1(4/3)\pi R^3`$ is the volume of the contracted incident nuclei when superposed. Here $`\zeta `$ is the inelasticity parameter and $`\gamma _0=\sqrt{s}/2m`$. These initial conditions (together with EOS-I) can reproduce both the particle multiplicities and the rapidity distributions fairly well (Ref. MG2b ). However it has an inconvenience that the initial energy density is too high. As a consequence, the time of expansion becomes very large and transverse expansion is considerable. The result is a very large effective temperature and a non-exponential shape of the spectra in transverse momentum.

A solution we found to this problem was to give the fluid an initial longitudinal-velocity distribution. We did not change the shape of the region where the fluid is formed, which remained a Lorentz contracted sphere. The initial longitudinal velocity is not constant, but proportional to $`z`$:

$`u^\mu (z)`$ $`=`$ $`(\gamma (z),0,0,\gamma (z)v),`$ (19)
$`u^z(z)`$ $`=`$ $`{\displaystyle \frac{z}{z_m}}\zeta \sqrt{\gamma _0^21},`$ (20)
$`\gamma (z)`$ $`=`$ $`\sqrt{\left({\displaystyle \frac{z}{z_m}}\right)^2\zeta ^2(\gamma _0^21)+1},`$ (21)
where $`z_m=R/\gamma _0`$. As a consequence of this change, the initial energy density $`\epsilon _0`$, which we took constant and obtained solving the equation

$$E=\underset{V:r^2+\gamma _0^2z^2=R^2}{}\left[(\epsilon _0+p)\gamma ^2(z)p\right]𝑑V,$$
(22)
became much smaller and then the transverse expansion compatible with data, as shown in Table 1 and Fig. 1. In this simulation we consider freeze-out temperatures smoothly increasing with energy, until SPS domains. For RHIC energies, we show two cases: a decreasing freeze-out temperature and a increasing one (bold-faced types in Table 1 and triangle in Fig. 1).

However, since $`\epsilon _0`$ became much smaller than the previous case, so did the entropy density $`s_0`$ and the total multiplicity, although the rapidity distributions remained more or less the same in the shape.

To reproduce all quantities (multiplicities, rapidity distributions and transverse-momentum distributions) we have to take a large initial volume, with some appropriate longitudinal-velocity distribution. To do this, we use the NeXus event generator, which produces the event-by-event initial conditions at time $`\tau =1`$ fm/c. Instead of using these fluctuating initial conditions, here we smooth them out by averaging 30 random events for each collision energy. For this case we use the EOS-II. The results are shown in Table 2 and Fig. 2. Besides reproducing the effective temperature data quite well, it was shown that this initial conditions give good results also of multiplicity and rapidity distributions HGSAKPTO .

## V Conclusions

We conclude that the initial conditions given in small volume like Landau-type discussed above are not appropriated to describe the multiplicity data, rapidity distributions and the transverse-momentum distributions at the same time. If we try to obtain the correct multiplicities (together with the rapidity distributions), the effective temperature becomes too large (the first case discussed in Sec. IV, where the fluid is at rest) . On the other side, if we try to get the correct effective temperature and the shape of rapidity distributions as we did, the multiplicities become small, unless an additional entropy is generated during the expansion (the second case discussed in Sec. IV). It is clear that these conclusions about the initial conditions do not change if we replace EOS-I by EOS-II.

Nexus type initial conditions, with spatially spread energy and velocity distributions, reproduce well all the main characteristics of data, namely particle multiplicities, rapidity distributions and $`p_T`$ spectra. In this case, the increase of $`T^{}`$ in RHIC energy domain seems to be due mostly to the larger expansion in the hadronic phase than to that in the QGP phase.

## Acknowledgments

This work was partially supported by FAPESP (Contract Nos. 2001/09861-1 and 2000/04422-7).