text:
# Ehrhart polynomials, simplicial polytopes, magic squares and a conjecture of Stanley

## 1. Introduction

A *magic square* is a square matrix with nonnegative integer entries having all line sums equal to each other, where a line is a row or a column. Let $`H_n(r)`$ be the number of $`n\times n`$ magic squares with line sums equal to $`r`$. The problem to determine $`H_n(r)`$ appeared early in the twentieth century . Since then it has attracted considerable attention within areas such as combinatorics, combinatorial and computational commutative algebra, discrete and computational geometry, probability and statistics . It was conjectured by Anand, Dumir and Gupta  and proved by Ehrhart  and Stanley  (see also \[23, Section I.5\] and \[24, Section 4.6\]) that for any fixed positive integer $`n`$, the quantity $`H_n(r)`$ is a polynomial in $`r`$ of degree $`(n1)^2`$. More precisely, the following theorem holds.

###### Theorem 1.1.

(Stanley ) For any positive integer $`n`$ we have

(1)
$$\underset{r0}{}H_n(r)t^r=\frac{h_0+h_1t+\mathrm{}+h_dt^d}{(1t)^{(n1)^2+1}},$$

where $`d=n^23n+2`$ and the $`h_i`$ are nonnegative integers satisfying $`h_0=1`$ and $`h_i=h_{di}`$ for all $`i`$.

It is the first conjecture stated in  (see Section I.1 there) that the integers $`h_i`$ appearing in (1) satisfy further the inequalities

(2)
$$h_0h_1\mathrm{}h_{d/2}.$$

In this paper we prove this conjecture by showing that $`(h_0,h_1,\mathrm{},h_d)`$ is equal to the $`h`$-vector of a $`d`$-dimensional simplicial polytope. Such vectors are known to be symmetric and unimodal and are characterized by McMullen’s $`g`$-theorem ; see  and Section 2.

A few comments on Stanley’s conjecture and the proof given in this paper are in order. It is known that the sequence $`(h_0,h_1,\mathrm{},h_d)`$ of Theorem 1.1 is a *Gorenstein sequence*, meaning it is the $`h`$-vector of a standard, graded, Gorenstein commutative ring; see for instance \[23, Section I.13\]. It is an important open problem to characterize Gorenstein sequences; see \[23, Section II.6\]. In this direction it was originally conjectured by Stanley  that a sequence $`(h_0,h_1,\mathrm{},h_d)`$ is Gorenstein if and only if it satisfies the conditions of the $`g`$-theorem but a counterexample was later given in . Our result (Corollary 3.6) and its generalization to the enumeration of magic labelings of regular bipartite graphs (Corollary 4.5) give an instance in which Stanley’s original conjecture turns out to be true. Another such instance, in which the entries of $`(h_0,h_1,\mathrm{},h_d)`$ count linear extensions of a naturally labeled poset by the number of descents, was given recently by Reiner and Welker . The polynomial $`H_n(r)`$ is the Ehrhart polynomial (see Section 2) of the Birkhoff polytope of doubly stochastic $`n\times n`$ matrices. We will state our main result in the context of Ehrhart polynomials of integer polytopes (Theorem 3.5) as well as that of enumerating solutions to systems of linear homogeneous Diophantine equations (Corollary 4.1) and will show that both situations of Theorem 1.1 and  appear as special cases.

This paper was largely motivated by the work  of Reiner and Welker. I am grateful to Volkmar Welker for encouraging discussions and to Jesús DeLoera, Victor Reiner, Francisco Santos and Richard Stanley for helpful suggestions.

## 2. Background

In this section we review some basic definitions and background on convex polytopes and their face numbers, triangulations and Ehrhart polynomials. We refer the reader to the texts by Stanley , Sturmfels  and Ziegler  for more information on these topics. We denote by $``$ the set of nonnegative integers.

Face enumeration. Given a finite (abstract or geometric) simplicial complex $`\mathrm{\Delta }`$ of dimension $`d1`$, let $`f_i`$ denote the number of $`i`$-dimensional faces of $`\mathrm{\Delta }`$, so that $`(f_0,f_1,\mathrm{},f_{d1})`$ is the $`f`$-*vector* of $`\mathrm{\Delta }`$. The polynomial

(3)
$$\underset{i=0}{\overset{d}{}}f_{i1}(x1)^{di}=\underset{i=0}{\overset{d}{}}h_ix^{di},$$

where $`f_1=1`$ unless $`\mathrm{\Delta }`$ is empty, is the $`h`$-*polynomial* of $`\mathrm{\Delta }`$, denoted $`h(\mathrm{\Delta },t)`$. The $`h`$-*vector* of $`\mathrm{\Delta }`$ is the sequence $`(h_0,h_1,\mathrm{},h_d)`$ defined by (3).

A *polytopal complex* $``$ \[28, Section 8.1\] is a finite, nonempty collection of convex polytopes such that (i) any face of a polytope in $``$ is also in $``$ and (ii) the intersection of any two polytopes in $``$ is either empty or a face of both. The elements of $``$ are its *faces* and those of dimension $`0`$ are its *vertices*. The *dimension* of $``$ is the maximum dimension of a face. The complex $``$ is *pure* if all maximal faces of $``$ have the same dimension. The collection $`(P)`$ of all faces of a polytope $`P`$ and the collection $`(P)`$ of its proper faces are pure polytopal complexes called the *face complex* and *boundary complex* of $`P`$, respectively. Thus $`P`$ is *simplicial* if $`(P)`$ is a simplicial complex. The $`h`$-vectors of boundary complexes of simplicial polytopes are characterized by McMullen’s $`g`$-*theorem*  \[23, Section III.1\] \[28, Section 8.6\] as follows. A sequence $`(g_0,g_1,\mathrm{},g_{\mathrm{}})`$ of nonnegative integers is said to be an $`M`$-*vector* if

1. $`g_0=1`$ and
2. $`0g_{i+1}g_i^{(i)}`$ for $`1i\mathrm{}1`$,

where $`0^{(i)}=0`$ and

$$n^{(i)}=\left(\genfrac{}{}{0pt}{}{k_i+1}{i+1}\right)+\left(\genfrac{}{}{0pt}{}{k_{i1}+1}{i}\right)+\mathrm{}+\left(\genfrac{}{}{0pt}{}{k_j+1}{j+1}\right)$$

for the unique representation

$$n=\left(\genfrac{}{}{0pt}{}{k_i}{i}\right)+\left(\genfrac{}{}{0pt}{}{k_{i1}}{i1}\right)+\mathrm{}+\left(\genfrac{}{}{0pt}{}{k_j}{j}\right)$$

with $`k_i>k_{i1}>\mathrm{}>k_jj1`$, if $`n1`$. A sequence $`(h_0,h_1,\mathrm{},h_d)`$ of nonnegative integers is the $`h`$-vector of the boundary complex of a $`d`$-dimensional simplicial polytope if and only if

1. $`h_i=h_{di}`$ for all $`i`$ and
2. $`(h_0,h_1h_0,\mathrm{},h_{d/2}h_{d/21})`$ is an $`M`$-vector.

In particular $`(h_0,h_1,\mathrm{},h_d)`$ is symmetric and unimodal and hence satisfies the inequalities (2), known as the *Generalized Lower Bound Theorem* for simplicial polytopes.

Triangulations and Ehrhart polynomials. A *triangulation* of a polytopal complex $``$ is a geometric simplicial complex $`\mathrm{\Delta }`$ with vertices those of $``$ and underlying space equal to the union of the faces of $``$, such that every maximal face of $`\mathrm{\Delta }`$ is contained in a face of $``$. A triangulation of the face complex $`(P)`$ of a polytope $`P`$ is simply called a triangulation of $`P`$.

For any set $`\sigma `$ consisting of vertices of the polytopal complex $``$ we denote by $`\sigma `$ the subcomplex of faces of $``$ which do not contain any of the vertices in $`\sigma `$ and write $`v`$ for $`\sigma `$ if $`\sigma `$ consists of a single vertex $`v`$. Given a linear ordering $`\tau =(v_1,v_2,\mathrm{},v_p)`$ of the set of vertices of $``$ we define the *reverse lexicographic triangulation* or *pulling triangulation* $`\mathrm{\Delta }()=\mathrm{\Delta }_\tau ()`$ with respect to $`\tau `$   \[27, p. 67\] as $`\mathrm{\Delta }()=\{v\}`$ if $``$ consists of a single vertex $`v`$ and

$$\mathrm{\Delta }()=\mathrm{\Delta }(v_p)\underset{F}{}\{\mathrm{conv}(\{v_p\}G):G\mathrm{\Delta }((F))\{\mathrm{}\}\}$$

otherwise, where the union runs through the facets $`F`$ not containing $`v_p`$ of the maximal faces of $``$ which contain $`v_p`$ and $`\mathrm{\Delta }(v_p)`$ and $`\mathrm{\Delta }((F))`$ are defined with respect to the linear orderings of the vertices of $`v_p`$ and $`F`$, respectively, induced by $`\tau `$. Equivalently, for $`i_0<i_1<\mathrm{}<i_t`$ the set $`\{v_{i_0},v_{i_1},\mathrm{},v_{i_t}\}`$ is the vertex set of a maximal simplex of $`\mathrm{\Delta }_\tau ()`$ if there exists a maximal flag $`F_0F_1\mathrm{}F_t`$ of faces of $``$ such that $`v_{i_j}`$ is the last vertex of $`F_j`$ with respect to $`\tau `$ for all $`j`$ and $`v_{i_j}`$ is not a vertex of $`F_{j1}`$ for $`j1`$. A different way to define $`\mathrm{\Delta }_\tau ()`$ is the following. For any vertex $`v`$ of $``$ let

$$\mathrm{pull}_v()=(v)\underset{F}{}\{\mathrm{conv}(\{v\}G):G(F)\{\mathrm{}\}\},$$

where the union runs through the facets $`F`$ not containing $`v`$ of the maximal faces of $``$ which contain $`v`$. If $`_0=`$ and $`_i=\mathrm{pull}_{v_{pi+1}}(_{i1})`$ for $`1ip`$ then $`_p`$ is a triangulation of $``$ which coincides with $`\mathrm{\Delta }_\tau ()`$. It follows from \[12, Theorem 2.5.23\] (see also \[6, p. 80\]) that if $``$ is the boundary complex of a polytope $`P`$ then $`\mathrm{pull}_v()`$ is the boundary complex of another polytope, obtained from $`P`$ by moving its vertex $`v`$ beyond the hyperplanes supporting exactly those facets of $`P`$ which contain $`v`$. This observation implies the following lemma.

###### Lemma 2.1.

The reverse lexicographic triangulation of the boundary complex of a polytope with respect to any ordering of its vertices is abstractly isomorphic to the boundary complex of a simplicial polytope of the same dimension.

A convex polytope $`P^q`$ is said to be a *rational* or an *integer* polytope if all its vertices have rational or integer coordinates, respectively. It is called a *0-1 polytope* if all its vertices are 0-1 vectors in $`^q`$. If $`P`$ is rational then the function defined for nonnegative integers $`r`$ by the formula

$$\mathrm{Ehrhart}(P,r)=\mathrm{\#}(rP^q)$$

is a quasi-polynomial in $`r`$, called the *Ehrhart quasi-polynomial* of $`P`$ \[24, Section 4.6\]. If $`P`$ is an integer polytope then this quasi-polynomial is actually a polynomial in $`r`$. Let $`A^q`$ be the affine span of the integer polytope $`P`$. A triangulation $`\mathrm{\Delta }`$ of $`P`$ is called *unimodular* if the vertex set of any maximal simplex of $`\mathrm{\Delta }`$ is a basis of the affine integer lattice $`A^q`$. We denote by $`\mathrm{\Delta }_\tau `$ the reverse lexicographic triangulation of an arbitrary polytope $`P`$ with respect to the ordering $`\tau `$ of its vertices. Following  we call such an ordering of the vertices of an integer polytope $`P`$ *compressed* if $`\mathrm{\Delta }_\tau `$ is unimodular and call $`P`$ itself *compressed* if so is any linear ordering of its vertices. The following lemma holds for any unimodular triangulation of $`P`$, although we will not need this fact here.

###### Lemma 2.2.

(\[22, Corollary 2.5\]) If $`P`$ is an $`m`$-dimensional integer polytope in $`^q`$ and $`\tau `$ is a compressed ordering of its vertices then

$$\underset{r0}{}\mathrm{Ehrhart}(P,r)t^r=\frac{h(\mathrm{\Delta }_\tau ,t)}{(1t)^{m+1}}.$$

If $`P^m`$ is an $`m`$-dimensional polytope and $`V`$ is any linear subspace of $`^m`$ then the *quotient polytope* $`P/V^m/V`$ is the image of $`P`$ under the canonical surjection $`^m^m/V`$. This is a convex polytope in $`^m/V`$ linearly isomorphic to the image $`\pi (P)`$ of $`P`$ under any linear surjection $`\pi :^m^{mdimV}`$ with kernel $`V`$. Recall that the simplicial join $`\mathrm{\Delta }_1\mathrm{\Delta }_2`$ of two abstract simplicial complexes $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$ on disjoint vertex sets has faces the sets of the form $`\sigma _1\sigma _2`$, where $`\sigma _1\mathrm{\Delta }_1`$ and $`\sigma _2\mathrm{\Delta }_2`$ and that $`h(\mathrm{\Delta }_1\mathrm{\Delta }_2,t)=h(\mathrm{\Delta }_1,t)h(\mathrm{\Delta }_2,t)`$. The following proposition is essentially Proposition 3.12 in .

###### Proposition 2.3.

Let $`P`$ be an $`m`$-dimensional polytope in $`^m`$ having a triangulation abstractly isomorphic to $`\sigma \mathrm{\Delta }`$, where $`\sigma `$ is the vertex set of a simplex not contained in the boundary of $`P`$. Let $`V`$ be the linear subspace of $`^m`$ parallel to the affine span of $`\sigma `$.

The boundary complex of the quotient polytope $`P/V^m/V`$ is abstractly isomorphic to $`(P)\sigma `$ and inherits a triangulation abstractly isomorphic to $`\mathrm{\Delta }`$.

Two compressed polytopes. (a) A real $`n\times n`$ matrix is said to be *doubly stochastic* if all its entries are nonnegative and all its rows and columns sum to 1. The set $`P`$ of all real doubly stochastic $`n\times n`$ matrices is a convex polytope in $`^{n\times n}`$ of dimension $`(n1)^2`$, called the *Birkhoff polytope* \[28, Example 0.12\]. It follows from the classical Birkhoff-von Neumann theorem that the vertices of $`P`$ are the $`n\times n`$ permutation matrices, so that $`P`$ is a 0-1 polytope. The Birkhoff polytope was shown to be compressed by Stanley \[22, Example 2.4 (b)\] (see also \[27, Corollary 14.9\]).

(b) Let $`\mathrm{\Omega }`$ be a poset (short for partially ordered set) on the ground set $`[m]:=\{1,2,\mathrm{},m\}`$. Recall that an (order) ideal of $`\mathrm{\Omega }`$ is a subset $`I\mathrm{\Omega }`$ for which $`i<_\mathrm{\Omega }j`$ and $`jI`$ imply that $`iI`$. Let $`\mathrm{\Omega }^0`$ be the poset obtained from $`\mathrm{\Omega }`$ by adjoining a minimum element $`0`$. The *order polytope*  of $`\mathrm{\Omega }`$, denoted $`O(\mathrm{\Omega })`$, is the intersection of the hyperplane $`x_0=1`$ in $`^{m+1}`$ with the cone defined by the inequalities $`x_ix_j`$ for $`i<j`$ in $`\mathrm{\Omega }^0`$ and $`x_i0`$ for all $`i`$. Thus $`O(\mathrm{\Omega })`$ is an $`m`$-dimensional convex polytope. The vertices of $`O(\mathrm{\Omega })`$ are the characteristic vectors of the nonempty ideals of $`\mathrm{\Omega }^0`$ \[25, Corollary 1.3\] so, in particular, $`O(\mathrm{\Omega })`$ is a 0-1 polytope (see \[25, Theorem 1.2\] for a complete description of the facial structure of $`O(\mathrm{\Omega })`$). Order polytopes were shown to be compresed by Ohsugi and Hibi \[13, Example 1.3 (b)\].

## 3. Special simplices

Throughout this section $`P`$ denotes an $`m`$-dimensional convex polytope in $`^q`$ with face complex $`(P)`$. Let $`\mathrm{\Sigma }`$ be a simplex spanned by $`n`$ vertices of $`P`$. We call $`\mathrm{\Sigma }`$ a *special simplex* in $`P`$ if each facet of $`P`$ contains exactly $`n1`$ of the vertices of $`\mathrm{\Sigma }`$. Note that, in particular, $`\mathrm{\Sigma }`$ is not contained in the boundary of $`P`$.

###### Example 3.1.

Let $`P`$ be the polytope of real doubly stochastic $`n\times n`$ matrices. If $`v_1,v_2,\mathrm{},v_n`$ are the $`n\times n`$ permutation matrices corresponding to the elements of the cyclic subgroup of the symmetric group generated by the cycle $`(\mathrm{1\; 2}\mathrm{}n)`$ (or any $`n`$ permutation matrices with pairwise disjoint supports) then $`v_1,v_2,\mathrm{},v_n`$ are the vertices of a special simplex in $`P`$. Indeed, each facet of $`P`$ is defined by an equation of the form $`x_{ij}=0`$ in $`^{n\times n}`$ and misses exactly one of $`v_1,v_2,\mathrm{},v_n`$.

###### Example 3.2.

Let $`\mathrm{\Omega }`$ be a poset on the ground set $`[m]:=\{1,2,\mathrm{},m\}`$ which is graded of rank $`n2`$ (we refer to \[24, Chapter 3\] for basic background and terminology on partially ordered sets) and $`P=O(\mathrm{\Omega })`$ be the order polytope of $`\mathrm{\Omega }`$ in $`^{m+1}`$. Let $`\mathrm{\Omega }^0`$ be the poset obtained from $`\mathrm{\Omega }`$ by adjoining a minimum element $`0`$ and for $`1in`$ let $`v_i`$ be the characteristic vector of the ideal of elements of $`\mathrm{\Omega }^0`$ of rank at most $`i1`$, so that $`v_i`$ is a vertex of $`P`$. Since a facet of $`P`$ is defined either by an equation of the form $`x_i=x_j`$ with $`i<j`$ in $`\mathrm{\Omega }^0`$ and $`i,j`$ in successive ranks or by one of the form $`x_i=0`$ for $`i\mathrm{\Omega }^0`$ of rank $`n1`$, it follows that $`v_1,v_2,\mathrm{},v_n`$ are the vertices of a special simplex in $`P`$.

###### Lemma 3.3.

Suppose that $`v_1,v_2,\mathrm{},v_n`$ are the vertices of a special simplex in $`P`$. If $`F`$ is a face of $`P`$ of codimension $`k`$ for some $`1kn1`$ and $`F`$ does not contain any of $`v_1,v_2,\mathrm{},v_k`$ then $`F`$ must contain $`v_i`$ for all $`k+1in`$.

###### Proof.

Let $`\mathrm{\Sigma }`$ be the special simplex with vertices $`v_1,v_2,\mathrm{},v_n`$. Any codimension $`k`$ face of a polytope can be written as the intersection of $`k`$ facets, so we can write $`F=F_1F_2\mathrm{}F_k`$ where the $`F_j`$ are facets of $`P`$. For each $`1ik`$ we have $`v_iF`$ and hence $`v_iF_j`$ for some $`j=j_i`$. Since $`\mathrm{\Sigma }`$ is special the integers $`j_1,j_2,\mathrm{},j_k`$ are all distinct and hence for each $`1jk`$ we have $`v_iF_j`$ for some $`1ik`$, which in turn implies that $`v_iF_j`$ for all $`k+1in`$. It follows that $`v_iF_1F_2\mathrm{}F_k=F`$ for all $`k+1in`$. ∎

###### Lemma 3.4.

Suppose that $`\tau =(v_p,v_{p1},\mathrm{},v_1)`$ is an ordering of the vertices of $`P`$ such that $`\sigma =\{v_1,v_2,\mathrm{},v_n\}`$ is the vertex set of a special simplex in $`P`$. Let $`\mathrm{\Delta }`$ be the abstract simplicial complex on $`\{v_{n+1},\mathrm{},v_p\}`$ defined by the reverse lexicographic triangulation of $`(P)\sigma `$ with respect to $`(v_p,v_{p1},\mathrm{},v_{n+1})`$.

1. The reverse lexicographic triangulation $`\mathrm{\Delta }_\tau `$ of $`P`$ is abstractly isomorphic to the simplicial join $`\sigma \mathrm{\Delta }`$.
2. $`\mathrm{\Delta }`$ is abstractly isomorphic to the boundary complex of a simplicial polytope of dimension $`mn+1`$.

###### Proof.

(i) Let $`\sigma _i=\{v_1,\mathrm{},v_i\}`$ for $`0in`$, so that $`\sigma _0=\mathrm{}`$ and $`\sigma _n=\sigma `$, and let $`\mathrm{\Delta }_i`$ denote the abstract simplicial complex on the set $`\{v_{i+1},\mathrm{},v_p\}`$ defined by the reverse lexicographic triangulation of $`(P)\sigma _i`$ with respect to the ordering $`(v_p,v_{p1},\mathrm{},v_{i+1})`$. To prove that $`\mathrm{\Delta }_0=\sigma _n\mathrm{\Delta }_n`$, which is the assertion in the lemma, we will prove that $`(P)\sigma _i`$ is pure $`(mi)`$-dimensional and that $`\mathrm{\Delta }_0=\sigma _i\mathrm{\Delta }_i`$ for all $`0in`$ by induction on $`i`$. This is obvious for $`i=0`$ so let $`1in`$. By induction, any maximal face $`F`$ of $`(P)\sigma _{i1}`$ is a codimension $`i1`$ face of $`P`$. Since $`F`$ does not contain any of the vertices $`v_1,\mathrm{},v_{i1}`$, by Lemma 3.3 we have $`v_iF`$. This implies that $`(P)\sigma _i`$ is pure $`(mi)`$-dimensional and that $`\mathrm{\Delta }_{i1}=v_i\mathrm{\Delta }_i`$. The last equality and the induction hypothesis $`\mathrm{\Delta }_0=\sigma _{i1}\mathrm{\Delta }_{i1}`$ imply that $`\mathrm{\Delta }_0=\sigma _i\mathrm{\Delta }_i`$, which completes the induction.

(ii) Let $`V`$ be the linear subspace of $`^q`$ parallel to the affine span of the vertices in $`\sigma `$ and $`P/V`$ be the corresponding quotient polytope of $`P`$, so that $`P/V`$ has dimension $`mn+1`$. Part (i) and Proposition 2.3 imply that $`\mathrm{\Delta }`$ is abstractly isomorphic to a reverse lexicographic triangulation of the boundary complex of $`P/V`$. This is in turn isomorphic to the boundary complex of a simplicial polytope of dimension $`mn+1`$ by Lemma 2.1. ∎

The following theorem is the key to the results in this paper.

###### Theorem 3.5.

Suppose that $`P`$ is an integer polytope and $`\tau =(v_p,v_{p1},\mathrm{},v_1)`$ is an ordering of its vertices such that:

1. $`\tau `$ is compressed and
2. $`\{v_1,v_2,\mathrm{},v_n\}`$ is the vertex set of a special simplex in $`P`$.

Then

$$\underset{r0}{}\mathrm{Ehrhart}(P,r)t^r=\frac{h(t)}{(1t)^{m+1}}$$

where $`h(t)=h_0+h_1t+\mathrm{}+h_dt^d`$ is the $`h`$-polynomial of the boundary complex of a simplicial polytope $`Q`$ of dimension $`d=mn+1`$, so that $`h(t)`$ satisfies the conditions in the $`g`$-theorem.

In particular $`h_i=h_{di}`$ for all $`i`$ and $`1=h_0h_1\mathrm{}h_{d/2}`$.

Moreover, $`Q`$ can be chosen so that its boundary complex is abstractly isomorphic to the reverse lexicographic triangulation of $`(P)\{v_1,\mathrm{},v_n\}`$ with respect to the ordering $`(v_p,v_{p1},\mathrm{},v_{n+1})`$.

###### Proof.

Let $`\sigma =\{v_1,v_2,\mathrm{},v_n\}`$ and let $`\mathrm{\Delta }`$ denote the reverse lexicographic triangulation of $`(P)\sigma `$ with respect to the ordering $`(v_p,v_{p1},\mathrm{},v_{n+1})`$. Lemma 2.2 guarantees that the proposed equation holds with $`h(t)=h(\mathrm{\Delta }_\tau ,t)`$. Part (i) of Lemma 3.4 implies that

$$h(\mathrm{\Delta }_\tau ,t)=h(\sigma \mathrm{\Delta },t)=h(\sigma ,t)h(\mathrm{\Delta },t)=h(\mathrm{\Delta },t),$$

since face complexes of simplices have $`h`$-polynomial equal to $`1`$, and the result follows from part (ii) of the same lemma. ∎

We now apply Theorem 3.5 to the Birkhoff polytope and to order polytopes of graded posets. Observe that our theorem does not apply to all integer polytopes since 0-1 polytopes with no regular unimodular triangulations are known to exist .

Magic squares and the Birkhoff polytope. Let $`P`$ be the polytope of real doubly stochastic $`n\times n`$ matrices. Observe that the polynomial $`\mathrm{Ehrhart}(P,r)`$ coincides with the function $`H_n(r)`$ of Theorem 1.1. Since $`P`$ is a compressed integer polytope of dimension $`(n1)^2`$, Theorem 3.5 and Example 3.1 imply immediately the following corollary.

###### Corollary 3.6.

For any positive integer $`n`$ we have

$$\underset{r0}{}H_n(r)t^r=\frac{h(t)}{(1t)^{(n1)^2+1}}$$

where $`h(t)=h_0+h_1t+\mathrm{}+h_dt^d`$ is the $`h`$-polynomial of the boundary complex of a simplicial polytope of dimension $`d=n^23n+2`$, so that $`h(t)`$ satisfies the conditions in the $`g`$-theorem.

In particular $`h_i=h_{di}`$ for all $`i`$ and $`1=h_0h_1\mathrm{}h_{d/2}`$.

In view of the last statement in Corollary 4.1, the polytope in the previous corollary can be constructed by pulling in an arbitrary order the vertices of the quotient of $`P`$ with respect to the affine span of the vertices $`v_1,v_2,\mathrm{},v_n`$, chosen explicitly as in Example 3.1.

Eulerian polynomials and equatorial spheres. Let $`\mathrm{\Omega }`$ be a graded poset on the ground set $`[m]:=\{1,2,\mathrm{},m\}`$ of rank $`n2`$. Let $`\mathrm{\Omega }_i`$ be the set of elements of $`\mathrm{\Omega }`$ of rank $`i1`$ for $`1in1`$ and $`(\mathrm{\Omega })`$ be the set of *linear extensions* of $`\mathrm{\Omega }`$, meaning the set of permutations $`w=(w_1,w_2,\mathrm{},w_m)`$ of $`[m]`$ for which $`w_i<_\mathrm{\Omega }w_j`$ implies $`i<j`$. We assume that $`\mathrm{\Omega }`$ is *naturally labeled*, meaning that the identity permutation $`(1,2,\mathrm{},m)`$ is a linear extension. The $`\mathrm{\Omega }`$-*Eulerian polynomial* is defined as

$$W(\mathrm{\Omega },t)=\underset{w(\mathrm{\Omega })}{}t^{\mathrm{des}(w)}$$

where

$$\mathrm{des}(w)=\mathrm{\#}\{i[m1]:w_i>w_{i+1}\}$$

is the number of descents of $`w`$. Following  we call a function $`g:\mathrm{\Omega }`$ *equatorial* if $`\mathrm{min}_{a\mathrm{\Omega }}g(a)=0`$ and for each $`2in1`$ there exist $`a_{i1}\mathrm{\Omega }_{i1}`$ and $`a_i\mathrm{\Omega }_i`$ such that $`a_{i1}<_\mathrm{\Omega }a_i`$ and $`g(a_{i1})=g(a_i)`$. An ideal $`I`$ or, more generally, a strictly increasing chain of ideals $`I_1I_2\mathrm{}I_k`$ in $`\mathrm{\Omega }`$ is *equatorial* if the characteristic function $`\chi _I`$ of $`I`$ or the sum $`\chi _{I_1}+\chi _{I_2}+\mathrm{}+\chi _{I_k}`$, respectively, is equatorial. The *equatorial complex* $`\mathrm{\Delta }_{eq}(\mathrm{\Omega })`$, introduced in , is the abstract simplicial complex on the vertex set of equatorial ideals of $`\mathrm{\Omega }`$ whose simplices are the equatorial chains of ideals in $`\mathrm{\Omega }`$.

The following theorem is proved in Corollary 3.8 and Theorem 3.14 of .

###### Theorem 3.7.

(Reiner–Welker ) Let $`\mathrm{\Omega }`$ be a naturally labeled, graded poset on $`[m]`$ having $`n1`$ ranks. The equatorial complex $`\mathrm{\Delta }_{eq}(\mathrm{\Omega })`$ is abstractly isomorphic to the boundary complex of a simplicial polytope of dimension $`d=mn+1`$ which has $`h`$-polynomial equal to the $`\mathrm{\Omega }`$-Eulerian polynomial $`W(\mathrm{\Omega },t)`$.

Hence $`W(\mathrm{\Omega },t)`$ satisfies the conditions in the $`g`$-theorem and, in particular, it has symmetric and unimodal coefficients.

Let $`P`$ be the order polytope of $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }^0`$ be the poset obtained from $`\mathrm{\Omega }`$ by adjoining a minimum element $`0`$. Recall that the vertices of $`P`$ are the characteristic vectors of the nonempty ideals of $`\mathrm{\Omega }^0`$. The order polytope comes with its canonical triangulation  \[15, Proposition 2.1\], which is a unimodular triangulation with maximal simplices bijecting to the linear extensions of $`\mathrm{\Omega }`$. This canonical triangulation is in fact the reverse lexicographic triangulation of $`O(\mathrm{\Omega })`$ with respect to any ordering $`(u_p,u_{p1},\mathrm{},u_1)`$ of its vertices such that $`i<j`$ whenever the ideal of $`\mathrm{\Omega }^0`$ defined by $`u_i`$ is strictly contained in that defined by $`u_j`$. We will use the following lemma.

###### Lemma 3.8.

Let $`v_i`$ be the characteristic vector of the ideal of elements of $`\mathrm{\Omega }^0`$ of rank at most $`i1`$ for $`1in`$. Let $`\sigma =\{v_1,\mathrm{},v_n\}`$ and $`\tau =(v_p,\mathrm{},v_{n+1})`$ be an ordering of the remaining vertices of $`P`$ such that $`i<j`$ whenever $`i,jn+1`$ and the ideal defined by $`v_i`$ is strictly contained in that defined by $`v_j`$.

The equatorial complex $`\mathrm{\Delta }_{eq}(\mathrm{\Omega })`$ is the abstract simplicial complex defined by the reverse lexicographic triangulation of $`(P)\sigma `$ with respect to $`\tau `$.

###### Proof.

Let $``$ denote the face complex of $`P`$ and let $`x_{\widehat{1}}=0`$ by convention. The maximal faces of $`\sigma `$ are the faces of $`P`$ defined by systems of equations of the form $`x_{i_s}=x_{j_s}`$ for $`0sn1`$ where (i) $`i_0=0`$ and $`j_{n1}=\widehat{1}`$, (ii) $`i_s\mathrm{\Omega }_s`$ for $`1sn1`$, $`j_s\mathrm{\Omega }_{s+1}`$ for $`0sn2`$ and $`i_s<_\mathrm{\Omega }j_s`$ for $`1sn2`$ and (iii) if $`j_s=i_{s+1}`$ for consecutive values $`s=a,a+1,\mathrm{},b1`$ of $`s`$ then the interval $`[i_a,j_b]`$ in $`\widehat{\mathrm{\Omega }}`$ consists only of the elements of the chain $`i_a<i_{a+1}<\mathrm{}<i_b<j_b`$. The statement of the Lemma follows from the description of the maximal faces of a reverse lexicographic triangulation $`\mathrm{\Delta }()`$ (see Section 2) and that of the maximal faces of $`\mathrm{\Delta }_{eq}(\mathrm{\Omega })`$ (see \[15, Proposition 3.5\]). We omit the details. ∎

*Proof of Theorem 3.7.* Let $`P`$ be the order polytope of $`\mathrm{\Omega }`$, as before. Observe that $`\mathrm{Ehrhart}(P,r)`$ is equal to the number of order reversing maps $`\rho :\mathrm{\Omega }\{0,1,\mathrm{},r\}`$. It follows from \[24, Theorem 4.5.14\] that

(4)
$$\underset{r0}{}\mathrm{Ehrhart}(P,r)t^r=\frac{W(\mathrm{\Omega },t)}{(1t)^{m+1}}.$$

Let the vertices $`v_1,v_2,\mathrm{},v_p`$ of $`P`$ and $`\tau =(v_p,\mathrm{},v_{n+1})`$ be as in Lemma 3.8. We checked in Example 3.2 that $`v_1,v_2,\mathrm{},v_n`$ are the vertices of a special simplex in $`P`$. Since $`P`$ is a compressed integer polytope (see Section 2) Theorem 3.5 applies and we have

$$\underset{r0}{}\mathrm{Ehrhart}(P,r)t^r=\frac{h(t)}{(1t)^{m+1}},$$

where $`h(t)=h_0+h_1t+\mathrm{}+h_dt^d`$ is the $`h`$-polynomial of a simplicial polytope of dimension $`d=mn+1`$ having, in view of Lemma 3.8, boundary complex abstractly isomorphic to $`\mathrm{\Delta }_{eq}(\mathrm{\Omega })`$. Comparison with (4) yields $`h(t)=W(\mathrm{\Omega },t)`$ and completes the proof. ∎

## 4. Rational polyhedral cones

In this section we state several corollaries of Theorem 3.5, including a generalization of Corollary 3.6 to magic labelings of bipartite graphs. Let $`\mathrm{\Phi }`$ be a $`p\times q`$ integer matrix with rank $`p`$ and $`E_\mathrm{\Phi }`$ be the monoid of vectors $`x^q`$ satisfying the homogeneous system of linear equations $`\mathrm{\Phi }x=0`$. We assume that $`E_\mathrm{\Phi }`$ is nonzero, so that $`qp+1`$, and set $`m=qp1`$. We denote by $`𝒞_\mathrm{\Phi }`$ the cone of vectors $`x=(x_1,x_2,\mathrm{},x_q)^q`$ satisfying $`\mathrm{\Phi }x=0`$ and $`x_i0`$ for all $`i`$. Let $`L(x)=a_1x_1+a_2x_2+\mathrm{}+a_qx_q`$ be a linear functional on $`^q`$ with $`a_i`$ for all $`i`$ such that the set

(5)
$$P=\{x𝒞_\mathrm{\Phi }:L(x)=1\}$$

is nonempty and bounded. Thus $`P`$ is an $`m`$-dimensional convex polytope. Observe that

(6)
$$\mathrm{Ehrhart}(P,r)=\mathrm{\#}\{xE_\mathrm{\Phi }:L(x)=r\}$$

for any nonnegative integer $`r`$. Let $`\overline{E}_\mathrm{\Phi }`$ be the submonoid of $`E_\mathrm{\Phi }`$ consisting of those elements having positive coordinates. We say that $`\beta \overline{E}_\mathrm{\Phi }`$ is the *unique minimal element* of $`\overline{E}_\mathrm{\Phi }`$ if we have $`\beta \gamma `$ coordinatewise for all $`\gamma \overline{E}_\mathrm{\Phi }`$.

###### Corollary 4.1.

Let $`P`$ be an $`m`$-dimensional polytope, as in (5). If

1. $`P`$ is an integer polytope,
2. $`\overline{E}_\mathrm{\Phi }`$ has a unique minimal element $`\beta `$ and
3. there exists a compressed ordering $`\tau =(v_p,v_{p1},\mathrm{},v_1)`$ of the vertices of $`P`$ such that $`v_1+v_2+\mathrm{}+v_n=\beta `$ for some $`n`$

then the conclusion of Theorem 3.5 holds.

###### Proof.

In view of Theorem 3.5 it suffices to show that $`\{v_1,v_2,\mathrm{},v_n\}`$ is the vertex set of a special simplex in $`P`$. Let $`\beta =(\beta _1,\beta _2,\mathrm{},\beta _q)`$ and let $`F`$ be a facet of $`P`$, so that $`F`$ is defined by an equation of the form $`x_k=0`$ for some $`1kq`$. We need to show that exactly one of $`v_1,v_2,\mathrm{},v_n`$ has positive $`k`$th coordinate. Clearly at least one of $`v_1,v_2,\mathrm{},v_n`$ has this property, since $`\beta _i>0`$ for all $`i`$. Assume on the contrary that at least two $`v_j`$ have positive $`k`$th coordinate, say $`v_1`$ and $`v_2`$, so that $`1\gamma _k<\beta _k`$ if $`v_1=(\gamma _1,\gamma _2,\mathrm{},\gamma _q)`$. Since $`F`$ is a facet of $`P`$ there exists a point $`x=(x_1,x_2,\mathrm{},x_q)`$ in the affine span of $`P`$, which we may assume to be rational, satisfying $`x_k<0`$ and $`x_i>0`$ for all $`ik`$. By replacing $`x`$ with a suitable integer multiple we find a point $`\alpha =(\alpha _1,\alpha _2,\mathrm{},\alpha _q)^q`$ satisfying $`\mathrm{\Phi }\alpha =0`$, $`\alpha _k<0`$ and $`\alpha _i>0`$ for all $`ik`$. We may choose a nonnegative integer $`t`$ so that $`0<\alpha _k+\beta _k+t\gamma _k<\beta _k`$ (with $`t=0`$ if $`\alpha _k+\beta _k>0`$). Then $`\alpha +\beta +tv_1`$ is in $`\overline{E}_\mathrm{\Phi }`$ and has $`k`$th coordinate strictly less than $`\beta _k`$, which contradicts the minimality of $`\beta `$. ∎

The condition that $`\overline{E}_\mathrm{\Phi }`$ has a unique minimal element is satisfied if $`(1,1,\mathrm{},1)^q`$ is in $`E_\mathrm{\Phi }`$ and is known to hold if and only if the semigroup ring $`R_\mathrm{\Phi }=𝐤[E_\mathrm{\Phi }]`$, generated over a field $`𝐤`$ by the monomials corresponding to elements of $`E_\mathrm{\Phi }`$, is Gorenstein; see for instance \[23, Section I.13\].

###### Corollary 4.2.

Let $`P`$ be an $`m`$-dimensional polytope, as in (5). If

1. $`P`$ is a compressed integer polytope,
2. $`R_\mathrm{\Phi }`$ is Gorenstein and
3. $`E_\mathrm{\Phi }`$ is generated as a monoid by the vertices of $`P`$

then the conclusion of Theorem 3.5 holds where, in the statement of the theorem, $`n`$ has the value $`L(\beta )`$ for the unique minimal element $`\beta `$ of $`\overline{E}_\mathrm{\Phi }`$.

###### Proof.

Let $`\beta `$ be the unique minimal element of $`\overline{E}_\mathrm{\Phi }`$, whose existence is guaranteed by (ii). Then (iii) implies that $`\beta =v_1+v_2+\mathrm{}+v_n`$ for some vertices $`v_1`$, $`v_2,\mathrm{},v_n`$ of $`P`$, which must be pairwise distinct. Because of (i) any ordering $`(v_p,v_{p1},\mathrm{},v_1)`$ of the vertices of $`P`$ satisfies the assumptions of Corollary 4.1. The result follows from this corollary observing that $`L(\beta )=n`$. ∎

General conditions on $`\mathrm{\Phi }`$ and $`P`$ which guarantee assumptions (i) and (iii) of Corollary 4.2 were given by Ohsugi and Hibi .

###### Corollary 4.3.

Let $`P`$ be an $`m`$-dimensional polytope, as in (5). If $`P`$ is a 0-1 polytope and $`R_\mathrm{\Phi }`$ is Gorenstein then the conclusion of Corollary 4.2 holds.

###### Proof.

Since $`P`$ is a 0-1 polytope, Theorem 1.1 and Lemma 2.1 in  imply, respectively, that $`P`$ is compressed and that $`E_\mathrm{\Phi }`$ is generated as a monoid by the integer points of $`P`$, which are exactly the vertices of $`P`$. Thus the assumptions of Corollary 4.2 hold. ∎

###### Remark 4.4.

Recall that a matrix is called *totally unimodular* if all its subdeterminants are equal to $`0,1`$ or $`1`$. Suppose that the linear functional $`L(x)=a_1x_1+a_2x_2+\mathrm{}+a_qx_q`$ satisfies $`a_i`$ for all $`i`$. It follows from the results of  that $`P`$ is an integer polytope if the matrix obtained from $`\mathrm{\Phi }`$ by adding the row $`(a_1,a_2,\mathrm{},a_q)`$ is totally unimodular. Hence, in view of Corollary 4.3, this statement and the assumptions that (i) $`P[0,1]^q`$ and (ii) $`R_\mathrm{\Phi }`$ is Gorenstein imply the conclusion of Corollary 4.1.

Magic labelings of graphs. Let $`G`$ be a graph (multiple edges and loops allowed) with $`p`$ vertices and $`q`$ edges and edge set $``$. A *magic labeling*  of $`G`$ of *index* $`r`$ is an assignment $`\mathrm{}:`$ of nonnegative integers to the edges of $`G`$ such that for each vertex $`v`$ of $`G`$ the sum of the labels of all edges incident to $`v`$ is equal to $`r`$, in other words,

$$\underset{e:ve}{}\mathrm{}(e)=r.$$

Let $`^{}`$ denote the real vector space with basis $``$ and let $`x_e`$ be the linear functional on $`^{}`$ dual to the basis element $`e`$. Let $`\mathrm{\Phi }`$ be a full rank, $`p^{}\times q`$ integer matrix with kernel the set of $`x^{}`$ satisfying the linear system of equations of the form

(7)
$$\underset{e:ve}{}x_e=\underset{e^{}:v^{}e^{}}{}x_e^{}$$

in $`^{}`$, where $`v,v^{}`$ are vertices of $`G`$. If $`L(x)`$ is any of the functionals on $`^{}`$ in (7) and $`P`$ is as in (5) then $`\mathrm{Ehrhart}(P,r)`$ counts the number $`H_G(r)`$ of magic labelings of $`G`$ of index $`r`$. It follows from \[17, Proposition 2.9\] and either \[22, Theorem 2.3\] (applied as in \[22, Example 2.4 (b)\] in the case of the Birkhoff polytope) or \[13, Theorem 1.1\] that conditions (i) and (iii) of Corollary 4.2 are both satisfied if the graph $`G`$ is bipartite (or, more generally, if it satisfies condition (iii) of \[17, Proposition 2.9\]). Clearly condition (ii) is satisfied if $`G`$ is regular, since then $`(1,1,\mathrm{},1)E_\mathrm{\Phi }`$. Assuming further that $`G`$ is connected we have $`p^{}=p2`$. The following corollary specializes to Corollary 3.6 when $`G`$ is the complete bipartite graph on two sets of vertices, each of size $`n`$.

###### Corollary 4.5.

For $`n1`$ and for any connected regular bipartite graph $`G`$ with $`p`$ vertices and $`q=np/2`$ edges we have

$$\underset{r0}{}H_G(r)t^r=\frac{h(t)}{(1t)^{m+1}}$$

where $`m=qp+1`$ and $`h(t)=h_0+h_1t+\mathrm{}+h_dt^d`$ is the $`h`$-polynomial of the boundary complex of a simplicial polytope of dimension $`d=mn+1`$, so that $`h(t)`$ satisfies the conditions in the $`g`$-theorem.

In particular $`h_i=h_{di}`$ for all $`i`$ and $`1=h_0h_1\mathrm{}h_{d/2}`$.


format: text

text:
# Application of the Variability → Luminosity Indicator to X-Ray Flashes

## 1. Introduction

X-ray flashes (XRFs) are a new type of fast transient source observed with the BeppoSAX Wide Field Cameras (WFCs) at a rate of about four per year (Heise et al. 2001). The rate for XRFs is thus $``$ 40% of the rate for gamma-ray bursts (GRBs).

XRFs are distinguished from Galactic transient sources by their isotropic spatial distribution and short ($``$ 10-100 sec) durations. Furthermore, they are distinguished from GRBs based on the fact that the are not detected above 40 keV by the BeppoSAX GRB Monitor – implying much larger $`L_x/L_\gamma `$ ratios than GRBs. Actually, the distribution of $`L_x/L_\gamma `$ ratios overlaps considerably with that of GRBs (Heise et al. 2001). In other respects, such as duration, temporal structure, spectrum and spectral evolution, XRFs exhibit properties that are qualitatively similar to the X-ray properties of GRBs. These similarities have led to the suggestion that the XRFs are in fact very “X-ray rich” GRBs (Heise et al. 2001).

However, no host galaxy nor redshift has yet been reported for an XRF. Therefore, it is not known whether XRFs are very high redshift GRBs, normal GRBs viewed off the axis of the jet, or GRBs that have unusually small energies.

## 2. Results

A possible Cepheid-like luminosity estimator for GRBs was suggested by Ramirez-Ruiz and Fenimore (2000) and has been developed further by Reichart et al. (2000; see also Reichart & Lamb 2001a,b). These authors have shown that there exists a correlation between a measure $`V`$ of the variability of the burst time history and the intrinsic isotropic-equivalent peak photon energy luminosity $`L`$ of the burst for the 19 GRBs for which some redshift information (either a redshift measurement or a redshift limit) exists. We show this correlation in Figure 1.

In this paper, we use the possible Cepheid-like luminosity estimator for the long-duration GRBs developed by Reichart et al. (2001) to estimate the intrinsic isotropic-equivalent peak photon energy luminosity $`L`$, and thus the redshift $`z`$, of 10 XRFs detected by the BeppoSAX WFCs (Heise et al. 2001) and the Burst and Transient Source Experiment (BATSE) on the Compton Gamma-Ray Observatory (Kippen et al. 2002). There is a significant limitation to what we can say, since these XRFs did not trigger BATSE and therefore we must use BATSE survey data for these ten events. The BATSE survey data is 1024 msec time resolution data, whereas all previous applications of our proposed variability measure have used 64 msec time resolution data. This means that, technically, we can only establish lower limits on the luminosities, and therefore the redshifts, of the ten XRFs. However, our experience applying our variability measure to $`900`$ GRBs in the BATSE 4B catalog (Paciesas et al. 1999) suggests that this may not be a severe limitation.

Table 1 lists our results. Figure 2 compares the distribution in the ($`z,L`$)-plane of the lower limits in luminosity $`L`$ and redshift $`z`$ for the ten XRFs with the distribution in the ($`z,L`$)-plane of the luminosities and redshifts of $`900`$ GRBs from the BATSE 4B catalog (Paciesas et al. 1999) estimated (Reichart & Lamb 2001a,b) using the variability measure $`V`$ of Reichart et al. (2001). The diagonal solid line shows the approximate BATSE detection threshold. Figure 2 shows that the distribution of the lower limits for the XRFs is indistinguishable from the distribution of the estimates for the GRBs.

The fact that all ten XRFs lie below but near the BATSE detection threshold is expected because none of the ten events triggered the BATSE instrument. This fact also indicates that the variability measures $`V`$ (and therefore the estimates of the luminosities $`L`$ and redshifts $`z`$) which we would obtain for the ten XRFs, using 64 msec time resolution data, would not differ much from the lower bounds which we have obtained using BATSE survey data, which has 1024 msec time resolution. Otherwise, the ten XRFs would lie above the BATSE detection threshold, contradicting the fact that they did not trigger the BATSE instrument.

These results suggest that the variability $``$ luminosity indicator probably works for XRFs, and that the lower bounds on the luminosities $`L`$ and redshifts $`z`$ that we have obtained for the ten XRFs do not differ much from the estimates of these quantities that we would find using higher time resolution data.

Figure 3 compares the cumulative distribution of the lower limits on $`z`$ for the XRFs with the cumulative distribution of the estimates of $`z`$ for the GRBs. The two distributions appear similar, and a KS test confirms that there is no significant difference between them.

## 3. Conclusions

Our results suggest that the variability $``$ luminosity indicator probably works for XRFs. Assuming this to be so, we find that the luminosity and redshift distributions of XRFs are consistent with those of long-duration gamma-ray bursts (GRBs), and therefore most XRFs are probably not very high redshift GRBs. The fact that XRFs and GRBs have similar luminosity and redshift distributions suggests that XRFs and long-duration GRBs are produced by a similar mechanism.

## References

Heise, J., in’t Zand, J. J. M., Kippen, R. M., and Woods, P. M. 2001, in Proceedings of the Rome Workshop on Gamma-Ray Bursts in the Afterglow Era, ed. E. Costa, F. Frontera & J. Hjorth, 16

Kippen, R. M., Woods, P. M., Heise, J., in’t Zand, J. J. M., Briggs, M.S. and Preece, R. D. 2002, in Proceedings of the Woods Hole Workshop on Gamma-Ray Bursts, ed. R. Vanderspek & G. R. Ricker (New York: AIP), in press

Paciesas, W. S., et al. 1999, ApJS, 122, 465

Ramirez-Ruiz, E. & Fenimore, E. E. 2000, ApJ, submitted (astro-ph/0004176)

Reichart, D., Lamb, D., Fenimore, E., Ramirez-Ruiz, E., Cline, T. & Hurley, K., et al. 2001, ApJ, 552, 57

Reichart, D. E. and Lamb, D. Q. 2001a, in Gamma-Ray Bursts in the Afterglow Era, ed. E. Costa, F. Frontera, and J. Hjorth (Springer-Verlag; Berlin), 233

Reichart, D. E. and Lamb, D. Q. 2001b, in Relativistic Astrophysics, AIP Conference Proceedings No. 586, ed. J. C. Wheeler and H. Martel (AIP: New York), p. 599


text:
# The X-ray Variability of NGC 4945: Characterizing the Power Spectrum through Light Curve Simulations

## 1 Introduction and Observations

X-ray variability is a prevailing feature of active galactic nuclei (AGN). With long data sets from monitoring campaigns by X-ray missions such as *Exosat*, *Ginga* and especially *RXTE*, it became possible to study the shape of the power density spectrum (PDS) over many decades of temporal frequency. Locally, the PDS is well-described by a power law, but globally, a suppression of power on long time scales has been found in a number of AGN (Papadakis & McHardy, 1995; Edelson & Nandra, 1999; Chiang et al., 2000; Nowak & Chiang, 2000; Pounds et al., 2001). Because of the direct proportionality between the mass and the Schwarzschild radius of a black hole, the timescale at which this suppression becomes important is expected to obey a linear relationship with the mass. Evidence for such a correlation has indeed started to emerge (Markowitz et al., 2003; McHardy et al., 2003) and can be extended over several orders of magnitude to galactic black hole candidates in their low/hard states, suggesting that the X-ray emission mechanism is similar for both classes of objects.

Our method for analyzing unevenly sampled light curves is based on simulation techniques used by a number of researchers (Markowitz et al., 2003; Done et al., 1992; Uttley, McHardy & Papadakis, 2002; Czerny et al., 2003), but introduces significant changes. We apply it to the Seyfert 2 galaxy NGC 4945, which is unique among the AGN for which a break in the PDS has been detected in that the mass of the central black hole of $`1.4\times 10^6`$ M is known fairly accurately from mapping the H2O maser emission (Greenhill, Moran & Herrnstein, 1997). A precise knowledge of the break frequency in this source should thus enable us to calibrate the relationship between the break frequency and the mass of the black hole.

NGC 4945 was observed by *RXTE* in late 2002 for a period spanning a total of 50 days. A pointing of 1,400 seconds on average was scheduled approximately every 6 hours, while for 7 days centered within the total observation period the source was monitored intensively. Data from the Proportional Counter Array (PCA) were reduced using standard *RXTE* PCA analysis tools. The average count rate and the time-averaged photon spectrum are consistent with previous observations (Done, Madejski & Smith, 1996; Madejski et al., 2000). The nuclear flux below 8 keV is heavily absorbed at the source, presumably by the same material that is responsible for the megamaser activity. The counts from all three layers in PCU 0 and 2 corresponding to the nominal range of energies from 8 to 30 keV were added to produce a combined light curve (16 s intrinsic binning), resulting in a total of usable data of about 300 ks. The light curve is shown in figure 1.

## 2 Data analysis

The algorithm of Timmer & Koenig (1995) is used to generate simulated light curves based on a user-chosen PDS model that are then sampled to match the observation times of the source light curve. Further processing comprises subtracting the mean, scaling each light curve to the intrinsic variance of the source light curve, and adding random Gaussian numbers multiplied by the error bars of the source light curve. We then calculate the Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982) and re-bin it to 5 points per decade. Using the ensemble of (typically $`>500`$) simulated power density spectra, we fit the distribution of the periodogram power in each frequency bin with a stretched $`\chi ^2`$ distribution: $`p_\nu ^s(\chi ^2,b)=b^1p_\nu (\chi ^2/b)`$ where $`p_\nu ^s`$ is the stretched, $`p_\nu `$ the regular $`\chi ^2`$ distribution, $`\nu `$ is the number of DOF, and $`b`$ is the stretch factor. Finally, we compare the source periodogram to the fitted distributions, using a $`\chi ^2`$-like statistic based on the likelihood ratio as the goodness-of-fit estimator (Baker & Cousins, 1984). The probability that an assumed model with the given set of parameters can be rejected is then the fraction of simulated light curves that have the likelihood ratio $`\chi ^2`$ lower than the source light curve. The process is illustrated as a diagram in Fig. 2, but further explanation is outlined below.

Scaling of simulated light curves to variance of source data (step D): In previous implementations of light curve simulations, the normalization of the PDS model was set to an arbitrary value at first. The best fit normalization was then determined by multiplying all simulated spectra by the same factor until the fit statistic was minimized (Markowitz et al., 2003; Uttley, McHardy & Papadakis, 2002). Due to the stochastic nature of the periodogram, the light curves simulated from a fixed normalization show a spread in variances (Vaughan et al., 2003), and the distribution of the power in each frequency bin will reflect this spread. By comparing the source data to a set of light curves simulated in this way, we are effectively asking the question: *How well do the chosen shape and normalization of the PDS model describe the variability of the source?* In reality, however, we are only interested in characterizing the shape of the periodogram. By allowing the normalization to vary from one light curve to the next, a slightly different question can be investigated: *How likely is it that the source light curve was produced by a process that has the chosen PDS shape?* By scaling each light curve to the intrinsic variance of the source light curve (and adding the Poisson noise level—see below), we ensure that the area under the curve for each simulated spectrum is equal to the one for the source.

Preservation of the Poisson noise level through the normalization of the PDS model (step E): The uncertainties in the source count rate are expected to contribute a constant level of power in the periodogram (usually called the Poisson noise level). Since the raw count rate in the *RXTE* PCA is dominated by the background for this source, the errors on the background-subtracted light-curve are approximately Gaussian. To mimic the effect of measurement errors in the simulated light curves, random Gaussian numbers multiplied by the error bars of the source light curve are added. This is done after the light curves have been scaled to the correct variance to ensure that the Poisson noise level in the simulated light-curves equals the one expected in the source data.

Use of robust fitting statistic to find best-fit parameters of the model PDS (steps H–J): A useful goodness-of-fit measure for non-Gaussian distributions is the rejection probability, calculated as the fraction of simulated light curves that have a value of the test statistic lower than the source data. The choice for test statistic is up to the investigator, and an obvious choice is the $`\chi ^2`$ statistic, which uses the ensemble average and standard deviation of the power in each frequency bin for comparison to the source power. Because of the non-Gaussian nature of the distribution of power in each frequency bin, the applicability of the $`\chi ^2`$ statistic is however severely limited: specifically, even small departures from Gaussian distributions have a marked effect on the standard deviation and thus on the rejection probability. Since the distribution of the power is more symmetric when plotted in the logarithm compared to the linear power, *substantially different goodness-of-fit values can be obtained simply by a change of variable.* The rejection probability, in effect, not only measures how well the model (and its associated parameters) describes the source data, but also the degree of non-Gaussianity. In contrast, by fitting the distributions of the power in each frequency bin with suitably parameterized functions (step H in the analysis), the effect of variable transformation on the fitting statistic is removed: If stretched $`\chi ^2`$ distributions are an adequate description of the distribution of linear power, then under a change of variable the data can equally well be fit with the similarly transformed functions. The probability density for observing the source power as well as the value at the peak of the distribution is unchanged under these transformations. *The likelihood ratio is thus an invariant quantity,* i.e. independent of the chosen form of the variable used to describe the power in the periodogram.

## 3 Results

We investigated two models: the simple (unbroken) power law (power law slope as only parameter), and the knee model of Uttley, McHardy & Papadakis (2002) ($`p(f)=A[1+(f/f_b)^2]^{\alpha /2}`$; parameters: high frequency power law slope $`\alpha `$, break frequency $`f_b`$). Figures 3 and 4 show the rejection probabilities obtained. The unbroken power law results in a minimum rejection probability of 70% at a slope of 1.3, indicating a rather poor fit. The rejection probability is reduced to 40% by introducing a $`10^6`$ Hz break at a high frequency slope of 1.5. It is clear that the exact location of the break is dependent on the particular PDS model used to characterize the source variability.

Uncertainty on rejection probability: To quantify the effect of simulating a finite number of light curves from which the goodness-of-fit statistic is estimated, we simulated more than 800,000 light curves using the unbroken power law model with a slope of 1.3. The process of fitting the distributions of power through calculating the rejection probability (steps H through J in Figure 2) was repeated for sets of light curves of increasing size. The mean of the distribution of rejection probabilities is consistent across the different sizes, while the standard deviation decreases from about 3% to 2.5% by going from 500 light curves to over 16,000. Assuming that the spread of rejection probabilities is only a weak function of the model and its associated parameters, we quote a characteristic uncertainty on our fitting statistic of about $`\pm 3`$%, based on 500 simulated light curves.

Simple power law vs. knee model: The difference in rejection probability between the unbroken power law and the knee model is much larger than the above uncertainty. The flattening of the power towards low frequencies in NGC 4945 has thus been detected at a statistically significant level. The length of our observation is most likely not sufficient to distinguish between models that use different functional forms for the break in the PDS. In particular, the differences in rejection probabilities between the knee model and the broken power law model used in Uttley, McHardy & Papadakis (2002) is not expected to be statistically significant for the currently available light curve of NGC 4945.


text:
# News on Strangeness at Ultrarelativistic Energies - Review of Microscopic Models

## 1 Introduction

The purpose of this overview is to present recent developments in the field of microscopic transport theory for ultrarelativistic heavy ion collisions, that is, for energies at SPS and RHIC. Special emphasis is put on the strangeness production. Moreover, we will particularly focus on the role of the $`\varphi `$-meson as well as on the importance of the color field properties on strange observables. Furthermore, we will also briefly address the topics baryon transport, excitation functions, multi-particle collisions, initial conditions, $`J/\mathrm{\Psi }`$ data, and particle interferometry.

## 2 Color Fields

An example of what is meant by the importance of field effects on the dynamics is given by the quark molecular dynamics (qMD) model . It is used to study microscopically the dynamics of colored quarks. The intention is to get an effective description of the non-perturbative, soft gluonic part of QCD. In particular, the hadronization process itself can be described explicitly. The initial quark distributions are obtained, for example, from the UrQMD model. The non-equilibrium dynamics of hadronization and the loss of correlations among quarks can be studied subsequently by qMD. The semi-classical model has a two-body color potential between the (anti)quarks. In addition, one has to define a dynamical hadronization criterion. The quarks interact via a Cornell potential with color matrices. Thus, they carry color, (flavor, spin, and isospin). The simple model Hamiltonian is

$$H=\underset{i=1}{\overset{N}{}}\sqrt{p_i^2+m_i^2}+\frac{1}{2}\underset{i,j}{}C_{ij}V(\left|𝐫_i𝐫_j\right|)$$

where $`N`$ are the number of quarks in the system. The potential $`V(r)`$ is linearly increasing at large distances (confinement) and has a Coulomb-type behavior at short distances. The Cornell-potential is

$$V(r)=\frac{3}{4}\frac{\alpha _s}{r}+\kappa r.$$

A further approximation is done by allowing only the diagonal colorless gluon exchange, that is, the quarks always keep their color. The color matrix elements $`C_{ij}`$ determine the sign and relative strength of the interaction depending on the color combination of the pair (for details see ).

The rather simple model is surprisingly able to describe a large amount of experimental data, for example, various particle spectra and ratios. Concerning strangeness, we show here the interesting feature of aquiring net strangeness in the quark phase. According to the idea of the strangeness distillation antistrange quarks hadronize earlier in a baryon rich system. Fig. 1 shows the number of net strange quarks (s-$`\overline{\mathrm{s}}`$) as a function of time for the baryon rich system in central Pb+Pb collisions at 30$`A`$GeV. Values up to approximately 10 are reached depending on the string tension. This indicates the possibility to study hyperon matter, hypernuclei, or even strangelets in these high baryon density collisions around the future GSI energies.    Figure 1: qMD calculation for Pb+Pb at 30$`A`$GeV by S. Scherer . The number of net strange quarks is shown as a function of time for different values of the string tension $`\kappa `$.

Another feature of the model is to study the details of the parton-hadron phase transition that is mainly driven by quark rearrangement. Late stage quark-antiquark pair production is a rather rare process in this model because the strong color fields that are needed are screened. Fig. 2 shows the distribution of the mean path length of quarks before hadronizing. Quarks from the same original hadron (solid line) recluster faster than quarks from different original hadrons (dotted line). The hadron formation follows an exponential decay of the quark cluster to hadrons. The reclustering exhibits a different decay length than the rearrangement process (2.2 vs. 4.8 fm). The amount of reclustering ($`50\%`$ for S+Au) demonstrates the survival of correlations in the quark system and possibly incomplete thermalization. This kind of detailed analysis also provides additional insight to the question of the mechanisms of hadronization (recombination vs. fragmentation ).

Figure 2: Quark Molecular Dynamics simulation of the hadronization in S(200 $`A`$GeV)Au collisions. Quarks originating from the same hadron hadronize after a shorter path length $`s`$ than quarks steming from different hadrons. Figure taken from .

## 3 $`\varphi `$-mesons, baryon transport and the sensitivity to field properties

The field properties also enter in the description of heavy ion collisions by means of cascade calculations that are based on hadron, (di)quark, and string degrees of freedom, as for example the UrQMD model . Here the initial particle production is dominated by string excitations and fragmentations, that is, longitudinally excited color flux tubes. In a dense colored medium the color field strength might be considerably enhanced. The particle formation process in hadronic collisions can be viewed as quantum tunneling of quark-antiquark and gluon pairs in the presence of a background color electric field. It is formed between two receding hadrons which are color charged by the exchange of soft gluons while colliding. In nucleus-nucleus collisions the color charges may be considerably greater than in nucleon-nucleon collisions due to the almost simultaneous interaction of several participating nucleons. With increasing energy of the target and projectile the number and density of strings grows, so that they start overlapping, forming clusters, which act as new effective sources for particle production . It has been predicted that the multiplicities of, for example, strange baryons or antibaryons should be strongly enhanced  once the color field strength grows. The abundances of (multiply) strange (anti)baryons in central Pb+Pb collisions at Cern-SPS , for example, can only be explained within the framework of microscopic model calculations  if the elementary production probability of $`s\overline{s}`$ pairs, which is governed in the string models  by the Schwinger mechanism  $`\mathrm{exp}(\pi m_q^2/2\kappa )`$, is considerably enhanced. This corresponds either to a dramatic enhancement of the string tension $`\kappa `$ (from the default $`1`$GeV/fm to $`3`$GeV/fm) or to quark masses $`m_q`$ that are reduced from their constituent quark values to current quark values as motivated by chiral symmetry restoration. In a fully analog way, $`\overline{p}`$ abundances can be explained since an enhanced string tension similarly leads to an increased production probability of antidiquark-diquark pairs which is needed to account for the experimentally observed yields . However, this argument, providing additional motivation, has to be reconsidered if multifusion processes , e.g., $`5\pi \overline{B}B`$, that are neglected in microscopic transport models based on $`2n`$ scatterings, contribute significantly to the baryon pair production. A variation of the string tension from $`\kappa =1`$GeV/fm to $`3`$GeV/fm increases the pair production probability of strange quarks (compared to light quarks) from $`\gamma _s=P(s\overline{s})/P(q\overline{q})=0.37`$ to 0.72. Similarly, the diquark production probability is enhanced from $`\gamma _{qq}=P(qq\overline{q}\overline{q})/P(q\overline{q})=0.093`$ to 0.45. In general, heavier flavors or diquarks ($`Q`$) are suppressed according to the Schwinger formula  by

$$\gamma _Q=\frac{P(Q\overline{Q})}{P(q\overline{q})}=\mathrm{exp}\left(\frac{\pi (m_Q^2m_q^2)}{2\kappa }\right).$$

The effectively enhanced string tension in a densely colored environment also follows from its relation to the Regge slope $`\alpha ^{}`$. In the rotational string picture the string tension $`\kappa `$ is related to the Regge slope $`\alpha ^{}`$ by

$$\kappa =\frac{1}{2\pi \alpha ^{}}.$$

The empirical value of the Regge slope for baryons is $`\alpha ^{}1`$GeV-2  that yields a string tension of approximately $`1`$GeV/fm. However, high-energetic processes dominated by Pomeron exchange, characterizing the multi-gluon exchange processes existent in high-energetic nucleus nucleus collisions, are described by a Regge trajectory (Pomeron) with a smaller slope of $`\alpha _P^{}0.4`$GeV-2 . According to Eq. (3) this translates into a considerably larger (effective) string tension $`\kappa `$.

Strong color fields also have an impact on the formation times of newly produced (di)quarks and thus on collision rates. The formation times are inversely proportional to the string tension

$$t_f1/\kappa .$$

The larger the string tension the shorter and short-living are the strings with a certain total energy.

Figure 3: Rapidity distributions of $`\varphi `$-mesons in central Au+Au collisions at RHIC. Besides the distributions of all $`\varphi `$’s (full circles), the contributions from $`\varphi `$’s produced in (resonant) meson-meson (K$`\overline{\mathrm{K}}`$) collisions (open circles) and from $`\varphi `$’s originating from string decays (open squares) are shown. Each panel corresponds to different values of the string tension ($`\kappa `$) (mass term) and the intrinsic transverse momentum broadening parameter $`\sigma `$.

Moreover, Schwinger’s formula does not only apply to the mass term but to the transverse momentum term, too. While $`\kappa `$ determines the production probability, the intrinsic transverse momentum parameter $`\sigma `$ defines the transverse momentum scale. The vacuum value of $`0.55`$GeV/c (width of a Gaussian) may be considerably enhanced by the strong color fields. Here, we vary both parameters independently to study systematically their individual effects on the observables . Fig. 3 shows the rapidity spectra of $`\varphi `$-mesons at RHIC. The intrinsic $`p_t`$ is varied from 0.55 to 1.5 GeV/c (left to right). The upper and lower panels correspond to the vacuum string tension and the SCF scenario, respectively. Clearly a larger $`\kappa `$ value leads to a strongly enhanced $`\varphi `$ production. Increasing the intrinsic $`p_t`$ decreases the $`\varphi `$-meson yield. The reason for this becomes immediatley clear when looking at the production channels. The production in meson-meson (K$`\overline{\mathrm{K}}`$) collisions gets considerably reduced when $`\sigma `$ is enhanced. The coalescence-like production loses due to a growing phase space. The experimental study of the $`\varphi `$-meson in both decay channels and the observed differences (due to rescattering effects) at RHIC will help to disentangle the different source contributions and thus enlight the role of strong color fields.

Another observable that has been shown to be highly sensitive to SCFs is the netbaryon number transport. Prior to the collision it is non-vanishing only at beam and target rapidities ($`\mathrm{\Delta }y11`$ units for RHIC!). Recent data by the STAR , PHOBOS , BRAHMS , and PHENIX  collaborations, show a considerable number of netprotons ($`p\overline{p}`$) at midrapidity. and antibaryon-to-baryon ratios smaller than unity.

The calculated anti-baryon-to-baryon ratios at midrapidity are shown in Fig. 4 as a function of the strangeness content $`|S|`$, both for $`\kappa =1`$GeV/fm as well as for the strong color field scenario (SCF, $`\kappa =3`$GeV/fm) . The $`\overline{B}/B`$-ratios increase with the strangeness content $`|S|`$ of the baryons and also increase with with impact parameter $`b`$. There is stronger absorption of antibaryons in central collisions (reducing the numerator of the ratio) and more stopping (increasing the denominator). Most important, the $`\overline{B}/B`$-ratios also increase with the color field strength $`\kappa `$ what is basically due to an increased pair production (while more collisions should reduce the ratios). The SCF results are closer to data. Also the observed hyperonization in the netbaryon distributions supports the SCF scenario (the $`(\mathrm{\Lambda }\overline{\mathrm{\Lambda }})/(p\overline{p})`$ ratio is $`0.8`$ for the strong color field (SCF) ($`\kappa =3`$GeV/fm) scenario, corresponding to the experimental value  but smaller for the vacuum calculations.

Figure 4: Antibaryon-to-baryon ratios at midrapidity as a function of the strangeness content $`|S|`$ in Au+Au collisions at RHIC ($`\sqrt{s}_{NN}=200`$GeV). Calculations with a string tension of $`\kappa =1`$GeV/fm are shown on the left and the results with strong color fields ($`\kappa =3`$GeV/fm) are shown on the right .

A complementary or alternative ansatz to describe the baryon transport is provided through the parton cascade model (PCM). Fig. 5 shows the netbaryon distributions as calculated with the VNI/BMS version  of the parton cascade. It is interesting to note that a considerable fraction of the observed netbaryon number at midrapidity is already contained in the parton distribution functions of the original nuclei. Other contributions arise from secondary scatterings of the partons and fragmentation that together yield a netbaryon number compatible with data. It will be interesting to see whether the correct longitudinal dynamics will also imply the strong transverse expansion dynamics as, for example, the relatively large elliptic flow.

Figure 5: Netbaryon rapidity distributions for Au+Au at RHIC. Crosses are calculation with only primary parton scatterings, rhombes include parton rescattering and squares include rescattering and parton fragmentation. The solid lines show the netbaryon content of the partonic distribution functions for gold nuclei, scaled by an average liberation factor of 0.4. Figure taken from .

The production of $`\varphi `$ mesons in Au+Au collisions at RHIC has also been studied with the quark gluon string model (QGSM) . It is found that the inverse slope parameter agree with data but the absolute yield of $`\varphi `$’s is underestimated by a factor 2. Somehow surprising, the fusion of strings does not increase the $`\varphi `$ yield in QGSM. Similar to previous UrQMD results  the early production is governed by string processes. The resonant kaon rescattering production channel dominates in the later stages. Fig. 6 illustrates the freeze-out distribution of the $`\varphi `$ mesons over rapidity and emission time from QGSM.

Figure 6: $`d^2N/dydt`$ distribution of the final-state $`\varphi `$ mesons in the $`(y,t)`$-plane from QGSM

The $`\varphi `$ puzzle has also been addressed with a systematic study of the kaon and dimuon decay channels by the AMPT model . The AMPT model is a multiphase transport model that includes the initial partonic and final hadronic interactions. It starts with minijet partons and strings from the HIJING mode. The partons enter the ZPC parton cascade model. The parton hadron transition is based on the Lund string fragmentation model. Final-state hadronic scatterings are modeled by the ART model. Similar to the other models the resonant K$`\overline{\mathrm{K}}`$ production channel also contributes considerably to the total $`\varphi `$ yield. However, even with in-medium modifications of the kaons and $`\varphi `$’s the observed differences in the yields cannot be fully explained (see figure 7). An additional ad hoc production is needed to come close to the data (of NA50). It is important to remark that the initial $`\varphi `$’s (those that are possibly created through strong color fields) will be primarily visible in the dimuon channel.

Figure 7: Transverse mass spectra of $`\varphi `$’s reconstructed from $`K^+K^{}`$ pairs (solid symbols) and from $`\mu ^+\mu ^{}`$ channel (open symbols) for Pb+Pb at $`158A`$ GeV ($`b3.5`$, $`|y|<1`$) fm in the AMPT model. The results are for without (squares) and with (triangles) in-medium mass modifications. The solid circles are NA49 data  for $`\varphi K^+K^{}`$ decay, and the open circles are NA50 data  for $`\varphi \mu ^+\mu ^{}`$. In the inset is shown as a function of $`m_T`$ the ratio $`R(m_T)`$ for phi mesons decaying to kaon-antikaon pairs that are not scattered to those determined from the dimuon channel. Circles are for an artificially increased number of $`\varphi `$’s by a factor 2 in the HIJING model. Figure from .

## 4 Further Topics

### 4.1 Excitation functions of meson abundancies and multi-meson fusion

New data from the Cern/SPS allow us to study the excitation function of particle production (and spectra, see ) in the interesting energy range of maximum netbaryon densities. Characteristic features in the data as maxima in the K+/$`\pi ^+`$ ratio gave reason to speculate about the observation of the phase transition. How do ’conventional’ hadron transport theories compare to data ? The detailed comparison with two realisations of the cascade model with hadron resonance, (di)quark, and string degrees of freedom (UrQMD and HSD) shows that the kaon yields are described reasonably (see Fig. 8). However, the pion yields and hence the K+/$`\pi ^+`$ ratios are off (see figure 9).

Figure 8: Excitation function of $`K^+`$ and $`K^{}`$ yields from central Au+Au (AGS) or Pb+Pb (SPS) collisions  in comparison to experimental data .

The overestimation of the pion yields in the model calculations could be due to a lack of multiparticle back reactions in the model that also lead to ’shorter’ chemical equilibration times. $`2n`$ collisions (with $`n>2`$) are performed but the inverse channels are not taken care of. Another possible explanation is given through enhanced pion in-medium masses. In fact, the importance of $`32`$ transitions has been explored in an extended HSD transport approach  for antiproton production by meson fusion in $`A+A`$ collisions at the AGS and SPS (see figure 10). In order to achieve a more conclusive answer from transport studies multi-particle interactions deserve further investigation in future generations of transport codes.

Figure 9: The same as in Fig. 8 for the $`K^+/\pi ^+`$ and $`K^{}/\pi ^{}`$ ratios.

Figure 10: Illustration of the flavor rearrangement model for $`B\overline{B}`$ annihilation to 3 mesons and vice versa. The mesons $`M_i`$ may be either pseudo-scalar or vector mesons, respectively. Figure from

Within this approach the abundancies of antiprotons as observed from peripheral to central collisions of $`Pb+Pb`$ at the SPS and $`Au+Au`$ at the AGS can approximately be described.

Fig. 11 shows the reaction rate $`B+\overline{B}mesons`$ vs. the backward reaction rate for central Pb+Pb collisions at 160 $`A`$GeV. Both rates turn out to be comparable within the statistics demonstrating an approximate local chemical equilibrium. The final $`\overline{p}`$’s origin dominantly ($``$ 88%) from the 3 meson fusion reactions. These studies indicate that these back reactions may have dramatic effects and therefore ask for further investigations.

Figure 11: The annihilation rate $`B\overline{B}`$ mesons (dashed histogram) for central $`Pb+Pb`$ collisions at 160 $`A`$GeV as a function of time in comparison to the backward reaction rate (solid histogram) within the HSD transport approach. Figure from .

### 4.2 J/$`\psi `$

Another possible QGP signal is charmonium production and its anomalous suppression. It is worth to mention that the new data from NA50 , however, do not deviate significantly from the predictions of standard-type cascade calculations (UrQMD or HSD)  (Fig. 12). Both the $`E_\mathrm{T}`$ dependence and the $`E_{\mathrm{ZDC}}`$ dependence do no longer support the observation of a strong anomalous suppression.

Figure 12: The ratio of $`J/\mathrm{\Psi }`$ to Drell-Yan production as a function of forward energy $`E_{\mathrm{ZDC}}`$ (left) and transverse energy $`E_\mathrm{T}`$ (right) for Pb+Pb at 158 AGeV. The preliminary data from the year 2000 NA50 Pb+Pb run  are compared to UrQMD results . Figure from

### 4.3 Initial conditions

The newly developed NEXUS3 model addresses the initial conditions of pp or AA collisions . It provides a self-consistent description of the energy sharing among the participating partons in the collision. The idea is called parton based Gribov-Regge theory, that is, soft and hard processes are described by (multi-)pomeron exchange. As one consequence strings are not necessarily attached to the valence quarks but may be connected to sea quarks. An observable supporting this picture may be given by the $`\overline{\mathrm{\Omega }}/\mathrm{\Omega }`$ ratio in pp collisions which is smaller or close to unity for NEXUS3 but larger for standard string models (or the older versions of NEXUS) .

### 4.4 Particle interferometry

In brief, progress has been made at several frontiers. More realistic models were developed and used to calculate the multidimensional correlation functions. Hybrid-type models, for example, that allow one to study explicitly the sensitivities to the QCD equation of state by simultaneously describing the freeze-out in a realistic way were used to discuss the so-called RHIC HBT-puzzle. The impact of in-medium modifications on HBT results were discussed for the first time. For a discussion and references, see for example .

## 5 Outlook

Topics not discussed here are, for example, the necessary tests of the microscopic models by simulating infinite matter. The upcoming RHIC pA data will put severe constraints for the models. An exact understanding of the elliptic flow $`v_2`$ is needed (collision numbers, formation times). How is the (partially) hydrodynamic behavior realized on a microscopic level? Also, there has been quite some progress in the charm sector, for example, the prediction of a substantial production of $`J/\mathrm{\Psi }`$’s in D$`\overline{\mathrm{D}}`$ collisions , in analogy to the $`\varphi `$ production via K$`\overline{\mathrm{K}}`$ .

## 6 Summary

Here, we focused on the role of strong color fields that lead to (i) enhanced heavy flavor and diquark production probabilities, (ii) modified, shorter formation times, (iii) stronger intrinsic transverse momenta. The effects on the $`\varphi `$ meson include (a) enhanced yields, (b) a change of the composition of the spectra, (c) a hardening of the spectra, (d) modified baryon dynamics (hyperonization). Moreover, we briefly discussed the competition of late stage vs. early production and the rescattering of the decay daughters. The differences between the kaonic and muonic $`\varphi `$-mesons are desired!

The quark molecular dynamics model provides a (simple) dynamical description of the hadronization process in the soft regime and needs to be further tested.

Parton cascades are now technically established and provide a valuable tool to learn about the validity range of pQCD or necessary modifications.

The study of excitation functions was shown to be essential to detect deficiencies in the models. Additional new developments include the investigation of multi-particle collisions, of charm in microscopic models, and the formulation of the initial conditions in a field theoretical language (parton ladders).

This work was supported by the Alexander von Humboldt-Foundation through a Feodor Lynen Fellowship, BMBF, DFG, GSI, and NSF grant PHY-03-11859. We thank S. Kesavan, V. Koch, J. Randrup, S. Scherer, H. Stöcker, and N. Xu, for many helpful comments and/or for providing figures.

## References


text:
# A possible hydrodynamical equivalence for astrophysical jets

## 1 Introduction

When astrophysical jets are formed, they begin to expand into their surrounding medium. A jet is capable of this expansion because the outflow, produced by the central engine, generates a ram pressure capable of elongating it. At the head of a particular jet, the kinetic energy of the particles is thermalised through a shock wave or *hot spot*. After a particular fluid particle crosses the hot spot, it is recycled back to the neighbourhood of the jet itself. This generates a cavity or *cocoon* that, together with the jet, expands through the environment. The cocoon keeps the jet collimated so that its transverse length does not grow too much as it expands. This standard picture for the expansion of an extragalactic (hydrodynamical) radio jet was first proposed by Scheuer (1974).

Today, we know that jets are not only associated with quasars, but that they occur at a range of different astrophysical scales. Some jets are observed in recently formed star systems. These are associated with Herbig–Haro objects and are commonly named Herbig–Haro jets (Reipurth & Bally, 2001). Jets can also form in galactic stellar binary systems (Mirabel & Rodriguez, 1994), and are so similar to their quasar extragalactic counterparts that they were termed $`\mu `$–quasar jets. Finally, observations of gamma ray burst afterglows (of the long variety) in distant galaxies, have led recently to the idea that jets are associated with these sources (see Mészáros (2002) for a review).

Aside from a difference of $`67`$ orders of magnitude covered by the maximum size of jets, their kinetic power output varies depending on the nature of the jet itself. Observations show that some extragalactic, and all Herbig–Haro jets, expand at non–relativistic velocities. On the other hand, the fluid particles inside quasar, $`\mu `$–quasar and $`\gamma `$–ray burst jets expand at relativistic velocities. Scheuer’s (1974) simple hydrodynamical model has been used as the basic model for all kind of jets that we observe in the universe.

Using dimensional analysis arguments, we propose in this article that there is a simple way to relate the maximum size of an astrophysical jet with the properties of its surrounding environment. The result appears to be valid both for low–velocity jets in the Newtonian regime and for highly relativistic outflows such as those found in $`\gamma `$–ray bursts.

Dimensional analysis techniques have been applied to the study of powerful jets in radio galaxies interacting with their external environment by Falle (1991) and Kaiser & Alexander (1997). However, the nature of their models show a self–similar indefinite growth of cocoon and jet, possibly due to the omission of gravity from the forces acting on the growing cocoon, since they concentrated their analysis on the initial growth phase.

## 2 Non–relativistic hydrodynamical jets

Let us first consider a non–relativistic jet, for example an AGN or a Herbig–Haro object jet, that has a characteristic length $`r`$. The physical processes and the dimensionless combination of important parameters that determine the dynamics of the expansion of the jet are many and complicated (Begelman et al., 1984; Blandford, 1990). For example, internal shock waves that develop inside jets, radiative cooling, kinetic power, magnetic fields, mass of the central engine and accretion power are some of the important physical quantities that enter into the description of an expanding jet. For the sake of simplicity, suppose that we do not include magnetic fields on the description of the jet’s dynamics and consider a purely hydrodynamical jet. Of all the possible dimensionless combinations of physical parameters that describe the problem, one was constructed by Mendoza & Longair (2001). This parameter measures the strength of the kinetic energy in the jet as compared to the gravitational energy of the gas cloud where it is embedded. In order to build such parameter, we proceed as follows: the gravity of the surrounding cloud is taken into account by introducing Newton’s gravitational constant $`G`$, and the physical properties of the cloud are determined by its average density $`\rho `$. In general terms, this density does not only represent the gas density of the surrounding material. It represents the *total* average density evaluated at a certain distance from the progenitor jet source. This means that the average density contains the gas density plus the stellar mass density (for galactic jets) and dark matter mass density (for extragalactic jets). Finally, since the expanding jet is roughly in pressure equilibrium with the surrounding cloud, the jet is roughly characterised by its bulk velocity $`v`$.

If dissipative processes inside the jet (such as radiative cooling and drag acting on the flow) are negligible then the four fundamental quantities ($`r,G,\rho \text{and}v`$) are a good estimate of the dynamical behaviour of the expanding jet in the cloud. This assumes that the expansion of the jet is close to ballistic and is of course a first order approximation, dealing only with the purely hydrodynamical aspect of the problem through dimensional analysis. One of the many simplifications being made is to consider a constant density medium, which is in all likelihood not true. However, under the reasonable assumption of a jet expanding from the centre of a spherically symmetric region, or from an overdensity within a cloud, the average density is well defined at each radius.

For any mechanical problem there are three independent dimensions, namely the dimensions of time $`t`$, length $`l`$, and mass $`m`$. Having restricted the problem to the above physics, Buckingham’s theorem of dimensional analysis (Sedov, 1993) demands the existence of a unique non–dimensional parameter $`\mathrm{\Lambda }`$ that will fully determine the solution of the problem. In the case of jets, this parameter is given by

$$\mathrm{\Lambda }G\frac{\rho r^2}{M^2a^2},$$
(1)
where $`Mv/a`$. Formally this number is not the Mach number of the jet flow since it is defined as the ratio of the jet velocity to the sound speed $`a`$ in the cloud (Mendoza & Longair, 2001). These authors showed that, for the particular problem they were discussing, the dimensionless parameter $`\mathrm{\Lambda }`$ measured in units of $`M^2`$ had the same value for jets formed in giant molecular clouds and for the gaseous haloes of galaxies (the length $`r`$ in their analysis was a characteristic length for which deflection of jets by pressure gradients might occur). This equivalence hints at an underlying mechanism that makes jets look the same at such different scales.

As noted by Mendoza & Longair (2001) the parameter $`\mathrm{\Lambda }`$ can also be written as

$$\mathrm{\Lambda }=\frac{3}{4\pi }\frac{G𝖬(r)/r}{v^2}1.2\times 10^{17}\left(\frac{𝖬(r)}{M_{}}\right)\left(\frac{c}{v}\right)^2\left(\frac{r}{\text{kpc}}\right)^1,$$
(2)
where $`𝖬(r)`$ is the total mass within a sphere of radius $`r`$ and $`c`$ is the speed of light. The right hand side of equation (2) is roughly the ratio of the gravitational potential energy to the kinetic energy of a fluid element of the jet evaluated at position $`r`$. Formally, the total gravitational potential energy acting on a fluid element on the jet needs to include not only the gas, but also the central mass and the dark matter (for the case of extragalactic jets) within radius $`r`$. Let us consider from now on that $`r`$ represents the maximum length that a specific class of jets can have. This maximum length is constrained by observations.

If a particular jet is such that $`\mathrm{\Lambda }1`$, then the total gravity acting on the jet is so strong that the jet cannot expand away from the central engine that is producing it. On the other hand, when $`\mathrm{\Lambda }1`$ the jet can freely expand away from the progenitor source, drilling a hole through the parent cloud. Since the numerical factor appearing in equation (2) is very small, it follows that for any real jet, the parameter $`\mathrm{\Lambda }1`$. For example, if we consider a Herbig–Haro jet, for which the total mass within a radius $`r`$ is a few $`M_{}`$ and the expansion velocity $`v10^3c`$, where $`c`$ is the speed of light, in which typical lengths $`r1\text{pc}`$, we then obtain a value of $`\mathrm{\Lambda }10^7`$. For the case of powerfull extragalactic jets, as will be discussed later, $`M(r)10^{15}`$, $`r1\text{Mpc}`$ and $`vc`$ which implies that $`\mathrm{\Lambda }10^5`$. As will be shown below, in the case of a relativistic jet we would replace $`\beta `$ by $`\gamma \beta `$, where $`\beta =v/c`$ and $`\gamma `$ is the Lorentz factor of the flow. This replacement reduces even more the value of $`\mathrm{\Lambda }`$.

Since $`\mathrm{\Lambda }1`$ for all types of jets (including $`\mu `$–quasar and $`\gamma `$–ray burst jets), it would appear that from merely considering gravitational potential and kinetic energies, jets should still be capable of growing several orders of magnitude in size beyond what is observed. Typical sizes of jets are much smaller than what the energy balance of $`\mathrm{\Lambda }1`$ suggests. This indicates that the mechanism determining the maximum extent of a jet must lie elsewhere.

In view of this fact, let us return to equation (1) in its original form, keeping the ratio $`v/a`$ explicit, and performing a re–scaling. The maximum length that a Herbig–Haro jet can have is $`r1\text{pc}`$ (Reipurth & Bally, 2001) and the surrounding medium where those jets expand are cold molecular clouds for which the particle number density $`n_\mathrm{H}10^2\text{cm}^3`$ and the temperature $`T10\text{K}`$ (Spitzer, 1998; Hartmann, 1998). We can thus write equation (2) as

$$\mathrm{\Lambda }\frac{10^1}{M^2}\left(\frac{M}{10\text{M}_{}}\right)\left(\frac{r}{1\mathrm{pc}}\right)^1\left(\frac{T}{10\mathrm{K}}\right)^1.$$
(3)
The total (stellar plus gas) mass $`𝖬(r)`$ within a radius of $`1\text{pc}`$ is such that $`2𝖬(r)10\text{M}_{}`$, so that the upper limit of equation (3) is $`\mathrm{\Lambda }10^1/M^2`$.

The largest extragalactic jets have typical lengths $`r2\text{Mpc}`$ (Ferrari, 1998) and their surrounding intergalactic medium is such that $`n_\mathrm{H}10^4\text{cm}^3,`$ and $`T10^8\text{K}`$ (Longair, 1992, 1998). The total gas mass within this radius is $`𝖬_{\text{gas}}10^{13}\text{M}_{}`$. Clusters of galaxies are such that the fraction of gas mass to dark matter is $`0.1`$ (Allen et al., 2002, 2004, and references within). This means that the total mass within a radius of $`2\text{Mpc}`$ is $`𝖬(r)10^{14}\text{M}_{}`$. With these values equation (3) implies that $`\mathrm{\Lambda }10^1/M^2`$, the same as the one obtained for Herbig–Haro jets.

The fact that the numerical value of $`\mathrm{\Lambda }`$ in units of $`M^2`$ is the same for Herbig–Haro jets and for non–relativistic extragalactic jets appears to indicate that a certain underlying physical mechanism is involved, fixing a maximum jet size. In order to investigate this, let us now rewrite equation (1) as

$$\mathrm{\Lambda }\frac{1}{M^2}\frac{r^2}{\lambda _\text{J}^2},$$
(4)
where $`\lambda _\text{J}a/\sqrt{G\rho }`$ is the cloud’s Jeans length. Note that $`\lambda _\text{J}`$ is the scale–length for gravitational instability in a fluid having sound speed $`a`$ and subject to the gravity produced by a total density $`\rho `$, which includes all matter present. The above comes from setting $`\rho _{\text{gas}}^1p_{\text{gas}}=\varphi `$, where $`\varphi `$ is the gravitational potential.

Comparing equation (4) with (3) we find that if the relation $`\mathrm{\Lambda }10^1/M^2`$ can be treated as universal, then

$$\lambda _\text{J}\sqrt{10}r.$$
(5)
One way of interpreting this general result is the following. It appears that jets can grow inside gaseous clouds, up to the point where they can no longer do so because a gravitational instability of the surrounding material sets in. In fact, the numerical factor of $`\sqrt{10}`$ that appears in equation (5) can be thought of as a form factor for the problem. The Jeans stability criterion can be applied to the gas inside the cocoon plus the swept up gas contained within the outer shell surrounding it in the presence of an extra central gravitational potential including dark matter or central star. The relevant pressure $`p_\text{c}`$ inside the cocoon decreases as $`p_\text{c}t^\alpha `$ with the expansion of the source in time $`t`$. Here $`\alpha 4/5\text{}\mathrm{\hspace{0.17em}2}`$ depending on the density profile of the external gas (Kaiser & Alexander, 1997). Since the gravitational force produced by the shocked gas inside the shell at the boundaries of the cocoon remains the same and because the pressure inside the cocoon tends to zero with the expansion of the source, there might be a time when a gravitational instability on the “ellipsoidal” cocoon will develop, leading to a collapse of the shell. In fact, the regions closest to the head of the jet will not collapse because $`\mathrm{\Lambda }1`$. However, regions of the shell closest to the central engine will eventually collapse for sufficiently large cocoon sizes, leading to the inflow of the external gas, which perhaps destroys the jet.

There is however another physical interpretation that can be given to equation (5). Alexander (2002) has proposed that the important instability which affects the late evolution in powerful FR–II radiogalaxies is the Rayleigh–Taylor instability (Frieman, 1954; Chandrasekhar, 1961). Due to the fact that Rayleigh–Taylor instabilities grow proportional to $`\mathrm{exp}\left(4\pi G\rho t\right)`$ in time $`t`$, they scale in exactly the same manner as the Jeans length when the wavelength $`\lambda `$ of the oscillations is such that $`\lambda \lambda _\text{J}`$. Once the gravitational force of the swept up gas is of the same order of the cocoon pressure, the interface between the swept up gas and the cocoon becomes Rayleigh–Taylor unstable leading to the break in of the surrounding medium, most likely inhibiting the growth of the jet.

The mean size of powerful radio sources is $`300\text{Kpc}`$ with maximum sizes reaching up to $`2\text{Mpc}`$. In the scenario presented here, these mean sized objects are understood as transient, being in the process of expansion. They will continue to grow until they reach the limit given by equation (5) provided the central engine remains active long enough. That no jets larger than this limit have been found, and the otherwise extraordinary coincidence of equation (5) at stellar and extragalactic scales suggest that the physics behind these objects is the same. Of course, other dimensionless parameters can be constructed, related to the many other physical aspects of the problem, such as magnetic fields or nature of the central engine. However, the numerical results of the previous section suggest that the particular feature of the problem which we are examining, i.e. maximum jet length, relates directly, or at least primarily, to $`\mathrm{\Lambda }`$ through equation (5).

## 3 Relativistic Jets–Generalising Further

When a fully relativistic jet is introduced in the discussion, the dimensional analysis is slightly more complex. This is because the speed of light $`c`$ has to be introduced as a fundamental parameter. In this case, the quantity $`M`$ is given by the ratio of the proper velocity of the jet $`\gamma v`$ to the proper velocity of sound $`\gamma _aa`$ in the medium (Chiu, 1973; Mendoza & Longair, 2002),

$$M=\frac{\gamma v}{\gamma _aa},$$
(6)
where $`\gamma =\left(1v^2/c^2\right)^{1/2}`$ and $`\gamma _a=\left(1a^2/c^2\right)^{1/2}`$ represent the Lorentz factors of the jet material and the sound speed respectively. Since all quantities, except the quantity $`M`$, that appear in equation (1) are measured in the proper frame of the gas, and because the sound speed of the gas in the cloud is non–relativistic, it is clear that the dimensionless parameter $`\mathrm{\Lambda }`$ given by equation (1) also has the same form in the relativistic case. The only difference is that the number $`M`$ is now defined by equation (6). Following the same procedure as before, it follows that equation (4) should be valid in this situation as well.

For the case of relativistic extragalactic jets, i.e. quasar jets, that expand through an intergalactic medium similar to their non–relativistic counterparts, it follows that $`\mathrm{\Lambda }10^1/M^2`$, so for this particular case equation (5) is also valid.

We can now repeat the calculations for the case of $`\mu `$–quasars to test the generality of the result. Most $`\mu `$–quasars show jets with typical lengths $`r1\text{pc}`$ (see Corbel et al., 2002, and references therein). The total (stellar plus gas) mass within $`1\text{pc}`$ is such that $`𝖬(r)10\text{M}_{}`$. Since the temperature of the interstellar medium is $`T10\text{K}`$ it follows that $`\mu `$–quasar jets satisfy the condition given by equation (6).

As a final application in the extreme, consider cosmological $`\gamma `$–ray bursts, GRBs (Fishman & Meegan, 1995; Mészáros, 2002). The relevant parameters are less well known than for the previous examples, but a comparison is nevertheless possible. These events typically release $`10^{51}`$ erg for a few seconds in a highly relativistic outflow, with Lorentz factors $`\gamma 1001000`$. The long variety, for which X–ray, optical and radio counterparts have been observed (van Paradijs et al., 2000), is probably related to the collapse of massive stars (Woosley, 1993; Stanek & et al., 2003; Hjorth & et al., 2003). The outflow itself is apparently highly collimated in many cases, manifesting itself as a relativistic jet, and the physical size of the source has actually been measured through radio scintillation, and is typically $`r10^21\text{pc}`$ (Frail et al., 1997, 2000). If the progenitors of GRBs are massive stars, it is reasonable to assume that they will be found in star–forming regions within the ISM. The density of the circumburst environment has also been determined from broadband observations in several cases, and ranges from $`0.1\text{cm}^3`$ to $`100\text{cm}^3`$, with a canonical value of $`10\text{cm}^3`$ (Panaitescu & Kumar, 2002). If the jet is directly interacting with the external ISM, the temperature could be as low as $`10\text{K}`$. This gives a value of $`\mathrm{\Lambda }`$ that does not violate our constraints, and is lower than in the previous cases. If however, the jet is within the wind–fed bubble produced by the massive progenitor star prior to the GRB, the temperature will be much higher, perhaps $`10^6\text{K}`$ (see e.g. Chevalier et al., 2004). Inspection of equation (3) shows that this will decrease $`\mathrm{\Lambda }`$ even further, thus also satisfying our proposed constraint that the right–hand side of equation (3) is an upper stability limit for $`\mathrm{\Lambda }`$. One additional fact makes GRB jets different from the other types considered here: the central engine typically lasts for only a few tens of seconds for long–duration bursts. It is most likely the case that this is what determines the extent of the jet, but it nevertheless does so within the limits imposed by equation (3) under the above interpretation.

## 4 Discussion

The results presented in this article suggest that astrophysical jets at all scales seem to obey the same physical law given by equation (5), regardless of their surrounding environment or the nature of the jet itself. It is remarkable that a class of objects where the physics is super relativistic and where the overall phenomena have a duration of milli–seconds, over $`15`$ orders of magnitude below that of large radio galaxies, appear to satisfy (or at any rate not to deviate significantly from) the predictions of equation (5).

The idea of unification of jets at all physical scales that expand in very different environments has of course been considered before. Much work has been made regarding the mechanism by which jets are *produced*, from $`\mu `$–quasar to AGN, with the magnetic field usually playing an important role (Blandford & Payne, 1982; Meier et al., 2001; Ghisellini & Celotti, 2002; Price et al., 2003). Using the observed time variability, the scaling is in terms of the mass of the central object, be it a black hole or a young star (e.g. Mirabel & Rodriguez (1994)). In this article we have proposed a unification scheme based on the physical properties between the expanding jet and its surroundings, described pictorically in Fig. 1. The straight line is given by equation (1), and we have argued here that real jets must lie above it.

Other physical causes of a jet’s maximum size are possible, most obviously the cut of an injection of particles and energy generated by the central engine. Also, for a purely hydrodynamical jet with a fixed opening angle $`\theta `$ that expands through an external gaseous environment at constant pressure, lateral expansion of the jet becomes important (Begelman et al., 1984). As a constant momentum is distributed over a growing area, ram pressure decreases, and a natural maximum extent might be reached through pressure equilibrium. Jet opening angles have been studied and observed on many astrophysical sources, mainly on extragalactic radio sources. These observations show globally well collimated jets over their large scale structure. Indeed, magnetic fields and internal shock waves (Begelman et al., 1984; Blandford, 1990) are probably important physical mechanisms that recollimate the jet flow at large scales. We hence believe that this explanation does not account properly for the maximum length of jets.

Also, well collimated self–similar jets and cocoons have been modelled very successfully (Falle, 1991; Kaiser & Alexander, 1997), albeit neglecting self-gravity of the medium into which the jet expands, other than in fixing the initial density profile of the surrounding medium. This very probably explains why no characteristic or limiting size for jets appears in their analysis (their jets actually continue to expand indefinitely). The analysis of the previous sections would indicate a limiting size which is perhaps due to gravitational instabilities in the growing cocoon, which collapses onto the jet. Most probably this instability is a Rayleigh–Taylor instability arising from the gravitational field produced by all matter within the cocoon interacting with the swept up gas as the source expands (Alexander, 2002).

To conclude, we propose in this paper that purely hydrodynamical jets of characteristic sizes are stable because there is a fundamental combination of definite physical quantities that limit their growth. Once the surrounding cocoon of jets reaches a certain limit, it probably becomes Rayleigh–Taylor unstable and might frustrate the growth of the jet. Evidently, additional physical phenomena (e.g. magnetic fields) play a role in the formation and evolution of astrophysical jets, and should be considered in a more detailed description.

## 5 Acknowledgements

We would like to thank Tigran Arshakian, Jorge Cantó, Manoj Choudhury and Enrico Ramirez–Ruiz for many fruitful discussions on $`\mu `$–quasars, $`\gamma `$–ray bursts and jets in general while preparing this letter. We very much thank the comments made by an anonymous referee, which gave a better physical understanding to the problem. S. Mendoza gratefully acknowledges financial support from CONACyT (41443) and DGAPA–UNAM (IN119203).