text:
# I. Introduction and Statements of the Results

## I. Introduction and Statements of the Results

Effects play very important role in the quantum theory of measurement (see, for example, ). In the Hilbert space formalism of the theory, the so-called Hilbert space effects are the positive bounded linear operators on a Hilbert space $`H`$ which are bounded by the identity $`I`$. The set of all effects on $`H`$ is denoted by $`E(H)`$.

There are several operations and relations defined on $`E(H)`$ which are important from different aspects of the theory. What concerns the present paper, here we are interested in the ortho-order structure and in the sequential structure on $`E(H)`$. The first one is obtained as follows. The usual ordering $``$ among self-adjoint bounded linear operators gives rise to a partial order on $`E(H)`$ and the operation $`{}_{}{}^{}:AIA`$ defines a kind of orthocomplementation on $`E(H)`$. This relation and operation together give the ortho-order structure on $`E(H)`$ . As for the second structure, it comes from the sequential product which is defined as follows. If $`A,BE(H)`$, then their sequential product is $`AB=\sqrt{A}B\sqrt{A}`$ .

Supposing $`dimH3`$, the automorphisms of $`E(H)`$ with respect to the ortho-order structure (called ortho-order automorphisms) as well as the ones with respect to the sequential product (called sequential automorphisms) are known to be implemented by unitary or antiunitary operators. This means that all those automorphisms are of the form

(1)
$$\varphi (A)=UAU^{}(AE(H))$$

where $`U`$ is a unitary or antiunitary operator. The result concerning ortho-order automorphisms was obtained by Ludwig in \[10, Section V.5\] (the proof was later clarified in ) while the corresponding result on sequential automorphisms was given by Gudder and Greechie in . In the paper  we showed that Ludwig’s theorem holds also in the 2-dimensional case.

In our recent paper  we presented some characterizations of the ortho-order automorphisms of $`E(H)`$ by means of their preserver properties. This investigation was motivated by the extensive study of preserver problems in matrix theory and in operator theory. Roughly speaking, preserver problems are concerned with the description of all transformations (called preservers) on a given set (preferably equipped with an algebraic structure) which preserve a certain relation between the elements, or a given subset of elements, or a quantity attached to the elements of the underlying set. In the present paper we continue the investigation started in  by studying the bijective maps $`\varphi :E(H)E(H)`$ which preserve the order and zero product in both directions, i.e., which satisfy

(a)
$$AB\varphi (A)\varphi (B)(A,BE(H))$$

and

(b)
$$AB=0\varphi (A)\varphi (B)=0(A,BE(H)).$$

In what follows we present a complete description of the structure of those maps and give several applications of the corresponding result. Among others, we obtain the form of all bijective transformations of $`E(H)`$ which

* preserve the order in both directions and
* map one single nontrivial scalar operator $`\lambda I`$ $`(\lambda 0,1)`$ to an operator of the same type.

At the first glance, it might be rather surprising that such maps are of a nice form but it turns out that they in fact satisfy (a) and (b). Next, we easily recover one of the main results in  on maps preserving the order and coexistency in both directions (the definition of coexistency is given below). What is probably more important, we also show that the main result of the present paper readily implies the former results on the structure of ortho-order automorphisms and sequential automorphisms of $`E(H)`$ given in  and in , respectively. In fact, we prove that those automorphisms also satisfy (a) and (b). It should be emphasized that, due to the proof of our main result, all statements presented in the paper are valid in 2-dimensional case as well. In particular, this holds for the result on the form of sequential automorphisms of $`E(H)`$ which is a new result.

Finally, in the last paragraph of the paper we demonstrate that even on the set $`E(𝒜)`$ of all effects belonging to a general von Neumann algebra, the ortho-order automorphisms and the sequential automorphisms (the definitions should be self-explanatory) belong to our new class of preservers, i.e., they preserve the order and zero product in both directions. This observation is worth mentioning since in that generality it is quite hard to see any connection between those two kinds of automorphisms. Therefore, we believe that our preservers deserve attention and it has sense to study them in other contexts as well.

Now we turn to the precise formulations of the results of the paper. *In what follows we assume that $`H`$ is a (complex) Hilbert space with $`dimH2`$.* Our main result reads as follows.

###### Theorem 1.

Let $`\varphi :E(H)E(H)`$ be a bijective map which preserves the order and zero product in both directions. Then there is an either unitary or antiunitary operator $`U`$ on $`H`$ and a real number $`p<1`$ such that with the function $`f_p(x)=\frac{x}{xp+(1p)}`$ $`(x[0,1])`$ we have

$$\varphi (A)=Uf_p(A)U^{}(AE(H)).$$

Here, $`f_p(A)`$ denotes the image of the function $`f_p`$ under the continuous function calculus belonging to the operator $`A`$.

To get the form of all bijective maps on $`E(H)`$ with the properties (i), (ii), we need the following proposition which might be interesting on its own right.

###### Proposition 2.

Let $`\varphi :E(H)E(H)`$ be a bijective map which preserves the order in both directions and suppose that there is a single pair of scalars $`\lambda ,\mu ]0,1[`$ such that $`\varphi (\lambda I)=\mu I`$. Then $`\varphi `$ preserves zero product in both directions.

The above results have the following immediate consequences. To the second statement in the corollary below observe that scalar operators in $`E(H)`$ can be characterized as effects which commute with every other effect.

###### Corollary 3.

Let $`\varphi :E(H)E(H)`$ be a bijective map which preserves the order in both directions and suppose that there is a single pair of scalars $`\lambda ,\mu ]0,1[`$ such that $`\varphi (\lambda I)=\mu I`$. Then there exist an either unitary or antiunitary operator $`U`$ on $`H`$ and a real number $`p<1`$ such that with the function $`f_p(x)=\frac{x}{xp+(1p)}`$ $`(x[0,1])`$ we have

$$\varphi (A)=Uf_p(A)U^{}(AE(H)).$$

In particular, we obtain the same form for any bijection of $`E(H)`$ which preserves the order and commutativity in both directions.

It is easy to see that if a function $`f_p`$ appearing above satisfies $`f_p(\lambda )=\lambda `$ for some $`\lambda ]0,1[`$, then $`p=0`$ and hence $`f_p`$ is the identity. This gives us the following corollary stating that if an order preserving bijection $`\varphi `$ of $`E(H)`$ fixes one single nontrivial scalar operator, then $`\varphi `$ is implemented by an either unitary or antiunitary operator.

###### Corollary 4.

Let $`\varphi :E(H)E(H)`$ be a bijective map which preserves the order in both directions and suppose that there is a scalar $`\lambda ]0,1[`$ such that $`\varphi (\lambda I)=\lambda I`$. Then there exists an either unitary or antiunitary operator $`U`$ on $`H`$ such that

$$\varphi (A)=UAU^{}(AE(H)).$$

The following corollary in the case when $`dimH3`$ appears in  as Theorem 1. Here we shall present a short proof of it based on our main result which works also in the 2-dimensional case. Recall that two effects $`A,BE(H)`$ are called coexistent if they are in the range of a positive operator valued measure or, equivalently, if there are effects $`E,F,GE(H)`$ such that

$$A=E+G,B=F+G,E+F+GE(H).$$

###### Corollary 5.

Let $`\varphi :E(H)E(H)`$ be a bijective map which preserves the order and coexistency in both directions. Then there is an either unitary or antiunitary operator $`U`$ on $`H`$ such that

$$\varphi (A)=UAU^{}(AE(H)).$$

As for the last two corollaries below, we shall see that the structure of the ortho-order automorphisms and that of the sequential automorphisms of $`E(H)`$ can be deduced from our main result even in the 2-dimensional case. The proofs are based on the observation that both kinds of automorphisms preserve the order and zero product in both directions. (For a similar observation concerning effects in general von Neumann algebras, see the last paragraph of the paper.) Recall that a bijective map $`\varphi :E(H)E(H)`$ is called an ortho-order automorphism if for any $`A,BE(H)`$ we have

$$AB\varphi (A)\varphi (B),$$

(i.e., $`\varphi `$ preserves the order in both directions) and

$$\varphi (A^{})=\varphi (A)^{}$$

(i.e., $`\varphi `$ preserves the orthocomplements). Moreover, a bijective map $`\varphi :E(H)E(H)`$ is called a sequential automorphism if

$$\varphi (AB)=\varphi (A)\varphi (B)$$

holds for every $`A,BE(H)`$.

###### Corollary 6.

Let $`\varphi :E(H)E(H)`$ be an ortho-order automorphism. Then $`\varphi `$ preserves the order and zero product in both directions. In fact, there is an either unitary or antiunitary operator $`U`$ on $`H`$ such that

$$\varphi (A)=UAU^{}(AE(H)).$$

###### Corollary 7.

Let $`\varphi :E(H)E(H)`$ be a sequential automorphism. Then $`\varphi `$ preserves the order and zero product in both directions. In fact, there is an either unitary or antiunitary operator $`U`$ on $`H`$ such that

$$\varphi (A)=UAU^{}(AE(H)).$$

Finally, we remark that the converse statements in all of the above results excluding Proposition 2 are also true. We mean that if a transformation is of the form which appears in the formulation of the corresponding result, then it has the properties which were required there. Since the verification of this observation needs only elementary computations, hence we omit them.

## II. Proofs

First we emphasize that the proof of our main result has many common points with the proof of \[12, Theorem\]. In fact, the argument to be presented below can be considered as an adaptation of the proof given in  for another situation. So, it would have been a possibility just to point out where and how the proof in  should be modified to get the statement of the main result of this paper but we think that such a proof is quite hard to follow and hence it does not meet the most elementary requirements. Hence, we decided to present a complete, self-contained proof.

We begin with some notation and useful facts that we shall apply in our arguments.

First, we note that the concept of the strength of an effect along a ray plays very important role in what follows. This notion was introduced by Busch and Gudder in . If $`A`$ is an effect on $`H`$, $`\phi `$ is a unit vector in $`H`$ and $`P_\phi `$ is the rank-one projection onto the subspace generated by $`\phi `$, then the quantity

$$\lambda (A,P_\phi )=sup\{\lambda [0,1]:\lambda P_\phi A\}$$

is called the strength of $`A`$ along the ray represented by $`\phi `$. Due to \[2, Theorem 4\] there is a very useful formula to compute the strength. In fact, we have

(2)
$$\lambda (A,P_\phi )=\{\begin{array}{cc}A^{1/2}\phi ^2,\hfill & \text{if }\text{ }\phi \mathrm{rng}(A^{1/2})\text{;}\hfill \\ 0,\hfill & \text{else}\text{.}\hfill \end{array}$$

(Here, $`\mathrm{rng}`$ denotes the range of an operator and $`A^{1/2}`$ denotes the inverse of $`A^{1/2}`$ on its range.)

Let $`\varphi :E(H)E(H)`$ be a bijective map which preserves the order in both directions. It was proved by Ludwig in \[10, Theorem 5.8., p. 219\] that $`\varphi `$ necessarily preserves the projections in both directions. It is then trivial to see that $`\varphi `$ also preserves the rank of the projections (cf. the proof of \[11, Theorem 1\]).

An easy fact follows what we shall use several times. Namely, if $`A,B`$ are effects, $`B`$ is of rank one and $`AB`$, then $`A`$ is a scalar multiple of $`B`$. This observation and the previous one have, among others, the following corollary. Let $`\varphi `$ be as above, i.e., suppose that it is a bijection of $`E(H)`$ which preserves the order in both directions. Then for every rank-one projection $`P`$, there is a function $`f_P:[0,1][0,1]`$ such that

$$\varphi (tP)=f_P(t)\varphi (P)(t[0,1]).$$

By the order preserving property of $`\varphi `$ and $`\varphi ^1`$ we see that $`f_P`$ is strictly increasing and bijective. In fact, we have

(3)
$$\varphi ^1(t\varphi (P))=f_P^1(t)P(t[0,1]).$$

Now, we turn to the proofs. In the proof of our main result Theorem 1 we need the following proposition which presents the solution of a functional equation.

###### Proposition 8.

Let $`f,g:]0,1[]0,1[`$ be functions and suppose that $`f`$ is a strictly monotone increasing bijection. Let

(4)
$$f\left(\frac{x}{x+(1x)y}\right)=\frac{f(x)}{f(x)+(1f(x))g(y)}(x,y]0,1[).$$

Then there are positive real numbers $`a,b,c`$ such that

$$f(x)=\frac{x^c}{x^c+a(1x)^c}(x]0,1[)$$

and

$$g(y)=by^c(y]0,1[).$$

###### Proof.

First note that since the function $`f`$ is continuous, equation (4) implies the continuity of $`g`$.

Next observe that with the notation

$$\begin{array}{cccc}\hfill \alpha (t)& =& \frac{1}{1+e^t}\hfill & (t),\hfill \\ \hfill \beta (x)& =& \mathrm{ln}\frac{1x}{x}\hfill & (x]0,1[),\hfill \\ \hfill \gamma (y)& =& \mathrm{ln}y\hfill & (y]0,1[),\hfill \end{array}$$

we have the identity

$$\frac{x}{x+(1x)y}=\frac{1}{1+\mathrm{exp}\left(\mathrm{ln}\frac{1x}{x}+\mathrm{ln}y\right)}=\alpha \left(\beta (x)+\gamma (y)\right)$$

for all $`x,y]0,1[`$. Therefore, equation (4) can be rewritten as

(5)
$$f\alpha (\beta (x)+\gamma (y))=\alpha (\beta f(x)+\gamma g(y))(x,y]0,1[).$$

Substituting $`x=\beta ^1(u)`$ and $`y=\gamma ^1(v)`$ into (5) and applying the inverse function of $`\alpha `$ to both sides of (5), we get

(6)
$$\alpha ^1f\alpha (u+v)=\beta f\beta ^1(u)+\gamma g\gamma ^1(v)$$

for all $`u`$ and $`v]\mathrm{},0[`$. This means that the functions

$$F=\alpha ^1f\alpha ,G=\beta f\beta ^1,\text{and}H=\gamma g\gamma ^1$$

satisfy the following so-called Pexider equation

$$F(u+v)=G(u)+H(v)(u,v]\mathrm{},0[).$$

Then, by known results of the theory of functional equations (cf. , or ) and by the continuity of $`F`$, $`G`$, $`H`$, it follows that there exist constants $`a,b,c`$ such that

$`F(w)`$ $`=`$ $`cw+a+b(w),`$
(7) $`G(u)`$ $`=`$ $`cu+a(u),`$
(8) $`H(v)`$ $`=`$ $`cv+b(v]\mathrm{},0[).`$
Using (7) and the definition of $`G`$, we get that $`\beta f(x)=c\beta (x)+a`$. Easy computation yields that

$$f(x)=\frac{x^c}{x^c+e^a(1x)^c}(x]0,1[).$$

Similarly, the definition of $`H`$ and $`\gamma `$, and equation (8) give

$$g(y)=e^by^c(y]0,1[).$$

The function $`f`$ being strictly increasing, $`G`$ is also increasing and hence we get $`c>0`$. ∎

Now, we are in a position to prove our main result.

###### Proof of Theorem 1.

The clue of the proof is to show that the functions $`f_P`$ (see the first part of this section) do not depend on the rank-one projections $`P`$. This will be done in what follows.

Let $`P,Q`$ be arbitrary mutually orthogonal rank-one projections. By the order, rank and orthogonality preserving properties of $`\varphi `$ on set of all projections we clearly have

$$\varphi (P+Q)=\varphi (P)+\varphi (Q).$$

Let $`\lambda [0,1]`$. From the inequality

$$\varphi (Q)\varphi (\lambda P+Q)\varphi (P+Q)=\varphi (P)+\varphi (Q)$$

we obtain

$$0\varphi (\lambda P+Q)\varphi (Q)\varphi (P).$$

As $`\varphi (P)`$ is of rank-one, according to the introduction of the present section, this implies that there is a scalar $`h_P(\lambda )[0,1]`$ such that

$$\varphi (\lambda P+Q)\varphi (Q)=h_P(\lambda )\varphi (P)$$

or, equivalently, that

$$\varphi (\lambda P+Q)=h_P(\lambda )\varphi (P)+\varphi (Q).$$

Since $`\varphi ^1`$ has the same properties as $`\varphi `$, it can be seen that the function $`h_P:[0,1][0,1]`$ is a strictly monotone increasing bijection. In fact, we have

(9)
$$\varphi ^1(\lambda \varphi (P)+\varphi (Q))=h_P^1(\lambda )P+Q.$$

We assert that $`h_P=f_P`$. Indeed, since

$$f_P(\lambda )\varphi (P)=\varphi (\lambda P)\varphi (\lambda P+Q)=h_P(\lambda )\varphi (P)+\varphi (Q),$$

it follows that $`f_Ph_P`$. By (3) and (9), if one considers $`\varphi ^1`$, it follows that $`f_P^1h_P^1`$. Since the functions $`f_P,h_P:[0,1][0,1]`$ are monotone increasing we then conclude that $`f_P=h_P`$. From the inequality

$$\begin{array}{c}f_P(\lambda )\varphi (P)=\varphi (\lambda P)\varphi (\lambda (P+Q))\\ \varphi (\lambda P+Q)=h_P(\lambda )\varphi (P)+\varphi (Q)=f_P(\lambda )\varphi (P)+\varphi (Q)\end{array}$$

we infer that

$$0\varphi (\lambda (P+Q))f_P(\lambda )\varphi (P)\varphi (Q).$$

As $`\varphi (Q)`$ is of rank one, this implies that

(10)
$$\varphi (\lambda (P+Q))=f_P(\lambda )\varphi (P)+k_P(\lambda )\varphi (Q)$$

holds for some scalar $`k_P(\lambda )[0,1]`$.

With the notation $`F=P+Q`$ it follows from the equality (10) that the operator $`\varphi (\lambda F)`$ is diagonizable with respect to any orthonormal basis in the range of the projection $`\varphi (F)=\varphi (P)+\varphi (Q)`$. We obtain that $`\varphi (\lambda F)`$ is a constant multiple of $`\varphi (F)`$. This gives us that we have $`f_P=k_P`$ and hence

$$\varphi (\lambda (P+Q))=f_P(\lambda )\varphi (P)+f_P(\lambda )\varphi (Q).$$

Let $`R`$ be a rank-one projection whose range is included in the subspace generated by $`\mathrm{rng}P`$ and $`\mathrm{rng}Q`$. Then we have

$$f_R(\lambda )\varphi (R)=\varphi (\lambda R)\varphi (\lambda F)=f_P(\lambda )(\varphi (P)+\varphi (Q))=f_P(\lambda )\varphi (F).$$

This gives us that

(11)
$$f_Rf_P$$

whenever $`P,R`$ are rank-one projections. It is then trivial that there is in fact equality in (11) and this proves that $`f_P`$ does not depend on $`P`$. Denote by $`f`$ this common function.

Our next claim is to show that $`f`$ satisfies a functional equation of the form (4). Fix mutually orthogonal rank-one projections $`P,Q`$ on $`H`$. Pick $`\mu ]0,1[`$ and let $`E=\mu P+Q`$. Take any rank-one projection $`R`$ on $`H`$ whose range is contained in the subspace generated by the ranges of $`P`$ and $`Q`$ and which is neither equal nor orthogonal to $`P`$. Using the formula (2), one can easily verify that

$$\lambda (E,R)=\frac{\mu }{\mu +(1\mu )\mathrm{tr}PR}.$$

By the definition of $`\lambda (E,R)`$ and the order preserving property of $`\varphi `$ it is clear that

$$\begin{array}{c}f(\lambda (E,R))=sup\{f(\lambda ):\lambda RE\}=sup\{f(\lambda ):\varphi (\lambda R)\varphi (E)\}=\\ sup\{f(\lambda ):f(\lambda )\varphi (R)\varphi (E)\}=\lambda (\varphi (E),\varphi (R)).\end{array}$$

Since $`\varphi (E)=\varphi (\mu P+Q)=f(\mu )\varphi (P)+\varphi (Q)`$, we obtain the equality

$$f\left(\frac{\mu }{\mu +(1\mu )\mathrm{tr}PR}\right)=\frac{f(\mu )}{f(\mu )+(1f(\mu ))\mathrm{tr}\varphi (P)\varphi (R)}.$$

As the quantities $`\mathrm{tr}PR`$ and $`\mathrm{tr}\varphi (P)\varphi (R)`$ do not depend on $`\mu `$, it follows from this equality that $`\mathrm{tr}\varphi (P)\varphi (R)`$ can be uniquely expressed as a function of $`\mathrm{tr}PR`$. Denoting $`g(\mathrm{tr}PR)=\mathrm{tr}\varphi (P)\varphi (R)`$, we get a bijective function $`g:]0,1[]0,1[`$ for which

$$f\left(\frac{\mu }{\mu +(1\mu )\nu }\right)=\frac{f(\mu )}{f(\mu )+(1f(\mu ))g(\nu )}(\mu ,\nu ]0,1[).$$

This gives us the desired functional equation for $`f`$ and $`g`$. By Proposition 8 we obtain that there are positive real numbers $`a,b,c`$ such that

$$f(x)=\frac{x^c}{x^c+a(1x)^c}(x]0,1[)$$

and

$$g(y)=by^c(y]0,1[).$$

But our function $`g`$ has the additional property that $`g(1x)=1g(x)`$ $`(x]0,1[)`$. In fact, this follows from the equality

$$\begin{array}{c}g(\mathrm{tr}PR)+g(1\mathrm{tr}PR)=g(\mathrm{tr}PR)+g(\mathrm{tr}QR)=\\ \mathrm{tr}\varphi (P)\varphi (R)+\mathrm{tr}\varphi (Q)\varphi (R)=1.\end{array}$$

One can easily deduce that we necessarily have $`b=1,c=1`$, i.e., $`g`$ is the identity on $`]0,1[`$. This further implies that our function $`f`$ is of the form

$$f(x)=\frac{x}{x+a(1x)}=\frac{x}{x(1a)+a}(x]0,1[).$$

Because of continuity, the above equality holds also on the whole interval $`[0,1]`$. Hence, we have that $`f`$ is of the form $`f=f_p`$ where $`p=1a`$.

Since the function $`g`$ above is the identity, we have

$$\mathrm{tr}PQ=\mathrm{tr}\varphi (P)\varphi (Q)$$

for all rank-one projections $`P,Q`$ on $`H`$. Hence, using Wigner’s celebrated theorem on symmetry transformations (sometimes called unitary-antiunitary theorem) we obtain that there exists an either unitary or antiunitary operator $`U`$ on $`H`$ such that

$$\varphi (P)=UPU^{}$$

holds for every rank-one projection $`P`$ on $`H`$. Consider the transformation

$$\psi :Af^1(U^{}\varphi (A)U)$$

on $`E(H)`$. It is not hard to see that this map is a bijection of $`E(H)`$ which preserves the order (as well as zero product) in both directions and it has the additional property that it fixes the so-called weak atoms, that is, the effects of the form $`\lambda P`$ where $`\lambda [0,1]`$ and $`P`$ is a rank-one projection. As, according to \[2, Corollary 3\], every effect is equal to the supremum of the set of all weak atoms it majorizes, we have that $`\psi `$ is the identity on $`E(H)`$. Transforming back, we see that

$$\varphi (A)=Uf(A)U^{}(AE(H)).$$

The proof is complete. ∎

###### Proof of Proposition 2.

Let $`\varphi :E(H)E(H)`$ be an order preserving bijection and $`\lambda ,\mu `$ be a pair of nontrivial scalars for which we have $`\varphi (\lambda I)=\mu I`$. Keeping the notation introduced in the first part of this section, we claim that

(12)
$$f_P(\lambda )=\mu ,$$

i.e., that $`\varphi (\lambda P)=\mu \varphi (P)`$ holds for every rank-one projection $`P`$. Indeed, from

$$f_P(\lambda )\varphi (P)=\varphi (\lambda P)\varphi (\lambda I)=\mu I$$

we deduce that $`f_P(\lambda )\mu `$. Now, considering $`\varphi ^1`$ and $`\varphi (P)`$ in the place of $`\varphi `$ and $`P`$, respectively, we also have $`f_P^1(\mu )\lambda `$. Since $`f_P`$ is increasing, this implies $`\mu f_P(\lambda )`$ and hence we get (12).

We next assert that $`\varphi `$ preserves the orthogonality between rank-one projections. To see this, let $`P,Q`$ be mutually orthogonal rank-one projections. Denote by $`P^{}=IP`$ the orthogonal complement of $`P`$. Consider the effect $`E=\lambda P+P^{}`$. Clearly, we have $`\lambda IEI`$, the strength of $`E`$ along $`P`$ is $`\lambda `$ and along $`Q`$ (which is a subprojection of $`P^{}`$) is $`1`$. It follows from the order preserving property of $`\varphi `$ and from (12) that

* $`\mu I\varphi (E)I`$,
* the strength of $`\varphi (E)`$ along $`\varphi (P)`$ is $`\mu `$,
* the strength of $`\varphi (E)`$ along $`\varphi (Q)`$ is $`1`$.

Now, Lemma 3 in  tells us that in this case the ranges of $`\varphi (P)`$ and $`\varphi (Q)`$ are subspaces of the eigenspaces of $`\varphi (E)`$ corresponding to the eigenvalues $`\mu `$ and $`1`$, respectively. This yields that the ranges of $`\varphi (P)`$ and $`\varphi (Q)`$ are orthogonal to each other.

As every projection is equal to the supremum of the set of all rank-one projections it majorizes, it follows that $`\varphi `$ preserves the orthogonality of projections of any rank. It is also easy to verify that $`\varphi `$ preserves the range projections of effects. This means that if $`R`$ is the range projection of $`A`$ (i.e., the projection onto $`\overline{\mathrm{rng}A}`$), then $`\varphi (R)`$ is the range projection of $`\varphi (A)`$. Indeed, this preserver property follows from the simple fact that the range projection of the effect $`A`$ is equal to the infimum of the set of all projections which are greater than or equal to $`A`$. It is clear that for any $`A,BE(H)`$, we have $`AB=0`$ if and only if the range projections of $`A`$ and $`B`$ are orthogonal. Using these observations we can infer that

$$AB=0\varphi (A)\varphi (B)=0.$$

This completes the proof. ∎

###### Proof of Corollary 5.

By \[11, Lemma 2\] an effect is coexistent with every other effect if and only if it is a scalar multiple of the identity. This implies that our transformation $`\varphi `$ maps scalar operators to scalar operators. Hence, by Corollary 3 we infer that up to unitary-antiunitary equivalence, $`\varphi `$ is of the form

$$\varphi (A)=f_p(A)(AE(H))$$

for some $`p<1`$. We claim that $`p=0`$, i.e., $`f_p`$ is the identity. To see this, let $`P,Q`$ be different rank-one projections. It was proved in \[11, Lemma 2\] that two rank-one effects with different ranges are coexistent if and only if their sum is an effect. Since $`\lambda P`$ and $`(1\lambda )Q`$ are always coexistent (indeed, their sum is an effect), we obtain that $`\varphi (\lambda P)`$ and $`\varphi ((1\lambda )Q)`$ must be also coexistent. By the just mentioned characterization we obtain that for any $`\lambda ]0,1[`$ we have

$$f_p(\lambda )P+f_p(1\lambda )Q=\varphi (\lambda P)+\varphi ((1\lambda )Q)I.$$

If we let $`Q`$ tend to $`P`$, we infer from this inequality that

(13)
$$f_p(\lambda )+f_p(1\lambda )1(\lambda [0,1]).$$

Since $`\varphi ^1`$ has the same properties as $`\varphi `$, it also follows that

$$f_p^1(f_p(\lambda ))+f_p^1(1f_p(\lambda ))1.$$

This further implies

$$f_p^1(1f_p(\lambda ))1\lambda $$

and by the monotonicity of $`f_p`$ we obtain

$$1f_p(\lambda )f_p(1\lambda ).$$

Comparing this with (13), we see that

$$f_p(1\lambda )=1f_p(\lambda )$$

holds for every $`\lambda [0,1]`$. It is easy to show that this implies $`p=0`$ which completes the proof. ∎

###### Proof of Corollary 6.

We first show that $`\varphi `$ preserves zero product in both directions. As $`\varphi `$ preserves the order in both directions, we know that $`\varphi `$ preserves the projections in both direction. Since $`\varphi `$ is an ortho-order automorphism, we obtain that it also preserves the orthogonality between projections. Now, one can use the same argument as in the last paragraph of the proof of Proposition 2 to verify that $`\varphi `$ preserves zero product in both directions. Since $`\varphi `$ preserves the orthocomplements, we can compute

$$\varphi \left(\frac{1}{2}I\right)=\varphi \left(\left(\frac{1}{2}I\right)^{}\right)=\varphi \left(\frac{1}{2}I\right)^{}$$

which implies that $`\varphi (\frac{1}{2}I)=\frac{1}{2}I`$. One can apply Corollary 4 to complete the proof. ∎

###### Proof of Corollary 7.

It was proved in  that the order on $`E(H)`$ is completely determined by the sequential product. More precisely, \[7, Theorem 5.1\] tells us that for any $`A,BE(H)`$ we have

$$ABCE(H):A=BC.$$

As $`\varphi `$ is a sequential automorphism, using this characterization it follows that $`\varphi `$ preserves the order in both directions. It is easy to see that

$$AB=(\sqrt{B}\sqrt{A})^{}\sqrt{B}\sqrt{A}=0AB=0.$$

Therefore, $`\varphi `$ also preserves zero product in both directions. By Theorem 1 we obtain that up to unitary-antiunitary equivalence, $`\varphi `$ is of the form

(14)
$$\varphi (A)=f_p(A)(AE(H))$$

for some $`p<1`$. We assert that $`p=0`$, i.e., $`f_p`$ is the identity. In fact, as $`\varphi `$ is a sequential automorphism, we obtain from (14) that $`f_p`$ is a multiplicative function on the unit interval. It is then obvious that $`p=0`$ and the proof is complete. ∎

In conclusion, according to our promise given in the introduction, we make some remarks on the relation between the class of our preservers and the collections of ortho-order automorphisms, resp. sequential automorphisms in the general setting of von Neumann algebras. So, let $`𝒜`$ be a von Neumann algebra of operators acting on the Hilbert space $`H`$. Denote by $`E(𝒜)`$ the set of all effects which belong to $`𝒜`$, i.e., let $`E(𝒜)=E(H)𝒜`$. The definitions of the order, orthocomplementation, and sequential product are straightforward and so are the definitions of ortho-order automorphisms and sequential automorphisms.

First, let $`\varphi `$ be an ortho-order automorphism of $`E(𝒜)`$. It is easy to see that the sharp elements in $`E(𝒜)`$ (i.e., the elements $`A`$ for which the infimum of $`A`$ and $`A^{}=IA`$ is $`0`$) are exactly the projections. Hence we obtain that $`\varphi `$ preserves the projections in $`E(𝒜)`$ as well as their orthogonality in both directions. Since, as it is well-known, the range projection of any element of a von Neumann algebra also belongs to the algebra, we see that $`\varphi `$ preserves the range projections of the elements of $`E(𝒜)`$ in the same sense as it was mentioned in the proof of Proposition 2. Then one can argue as in that proof to verify that $`\varphi `$ preserves zero product in both directions. So, we obtain that every ortho-order automorphism of $`E(𝒜)`$ belongs to our class of preservers, that is, those automorphisms preserve the order and zero product in both directions.

Now, let $`\varphi :E(𝒜)E(𝒜)`$ be a sequential automorphism. It is not hard to prove that the above mentioned characterization of the order by means of the sequential product due to Gudder and Greechie holds true also in the setting of von Neumann algebras. This means that for any $`A,BE(𝒜)`$ we have

(15)
$$ABCE(𝒜):A=BC.$$

In fact, if $`A,BE(𝒜)`$ and $`AB`$, then by \[7, Theorem 5.1\] there is an operator $`CE(H)`$ such that $`A=BC`$. Following the proof in , one can see that the existence of this operator $`C`$ is a consequence of a well-known result of Douglas . In fact, in the corresponding part of the proof of Douglas’ result this $`C`$ was constructed. Examining the construction, it is not hard to verify that $`C`$ belongs to the von Neumann algebra $`𝒜`$, i.e., we have $`CE(𝒜)`$. This gives one implication from the asserted equivalence in (15). The other implication is trivial. Then, just as in the proof of Corollary 7, one can show that $`\varphi `$ preserves the order and zero product in both directions.

To sum up, it has turned out that the ortho-order automorphisms and the sequential automorphisms all belong to our new class of preservers even in the setting of von Neumann algebras. This seems to be a worthwhile observation as in that generality it is quite hard to see any connection between those two kinds of automorphisms.

## Acknowledgements

This paper was written when the author held a Humboldt Research Fellowship. He is very grateful to the Alexander von Humboldt Foundation for providing ideal conditions for research, and to his host Werner Timmermann (TU Dresden, Germany) for very warm hospitality and, what concerns the present paper, for his kind advises that helped to improve the presentation of the manuscript. The author also acknowledges support from the Hungarian National Foundation for Scientific Research (OTKA), Grant No. T030082, T031995, and from the Ministry of Education, Hungary, Grant No. FKFP 0349/2000


format: text

text:
# Statistical Hadronization Probed by Resonances

## I Introduction

The Fermi statistical model of particle production Fer50 ; Pom51 ; Lan53 ; Hag65  has been used extensively in the field of relativistic heavy ion collisions jansbook . Particle abundancies and spectra both at Super Proton Synchrotron (SPS) PBM99 ; NA57 ; becattini ; van-leeuwen ; ourspspaper ; sqm2001 ; PBM01  and Relativistic Heavy Ion Collider (RHIC) burward-hoy ; castillo ; PBMRHIC ; bugaev\_freeze ; florkowski ; rafelski2002  energies have been analyzed in this way. The quality of fits to experimental results was such that it became possible to discuss hadronization conditions quantitatively, but the conclusions of the groups differ. For example the values of temperature range from as low as 110 MeV burward-hoy ; castillo ; van-leeuwen ; bugaev\_freeze  to 140 MeV ourspspaper ; rafelski2002  to as high as 160, 170, and 180 MeV florkowski ; PBM99 ; PBM01 ; PBMRHIC ; castillo ; becattini  at both SPS and RHIC energies.

These differences are on a closer inspection not very surprising, since the (tacit) assumptions made about hadronization mechanisms differ. However, this means that before we can say that the freeze-out temperature has been determined, we must understand precisely the origins of these differences, and proceed to ascertain which model is applicable. We shall suggest experimental observables which will be particularly sensitive to the differences between the hadronization scenarios, in the hope that further experimental study will allow to understand the statistical hadronization mechanism.

We begin with an overview of the differences between hadronization scenarios and their relation to the physical assumptions used. Every model discussed here has been extensively studied before, and has gained acceptance of some part of the heavy ion community. The theoretical principles that we invoke are well understood, and the methods we use can be found scattered in literature. We shall concentrate here on an analysis of resonances produced in a heavy ion collision. Direct detection of hadronic short-lived resonances has become possible through invariant mass reconstruction fachini ; vanburen ; markert ; friese . Resonances have already been proposed as candidates for differentiating between freeze-out models resonances . Resonances are a sensitive probe of the freeze-out temperature, since the ratio of yields of particles with the same quark composition is insensitive to both fugacities and phase space occupancies, and mass differences are greater than the hadronization temperature considered.

Here we shall develop this reasoning one step further: If particles were to be emitted from a static thermal source, and feed down corrections were performed, the ratio of resonance to daughter particle would be independent of $`m_{}`$. That this ratio is in general somewhat $`m_{}`$ dependent is in this situation due to dynamical effects, such as hadronizing matter flow and freeze-out geometry. For this reason, it can be expected that the $`m_{}`$ dependence of this ratio can help isolate these effects and thus remove the ambiguities in the present freeze-out models.

In Sec. II we review hadronization models and discuss their ambiguities. We then show in Sec. III how resonance $`m_{}`$ ratios can be used to distinguish between particle hadronization models. We close with a short discussion of open issues.

## II Statistical hadronization

### II.1 General remarks

Nearly all hadronic spectra comprise a significant O(50%) component from resonance decays. Fits to data, which are not allowing for the decay contributions have a very limited usefulness. Particle spectra and thus yields are in general, controlled by the properties of statistical hadronization models. However, some recent work fits the particle slopes only burward-hoy ; van-leeuwen ; NA57 ; castillo , treating the normalization of each particle as a free parameter. This approach can be argued for assuming a long-lived posthadronization “interacting hadron gas phase” in which individual hadron abundances subject to inelastic interactions evolve away from chemical equilibrium. This particular reaction picture clashes with e.g. the fact that short-lived resonance ratios can be described within the statistical hadronization model using the chemical (statistical hadronization) freeze-out temperature obtained in stable particle studies PBMRHIC . This implies that in principle the relative normalization of the particle spectra should be derived from a hadronization scenario involving flavor chemical potentials. In fact a study of RHIC spectra finds that the normalization can be accounted for florkowski , and that the chemical equilibration temperature also describes particle spectra well. This is suggesting that any posthadronization reinteraction phase is short and has minor influence on the particle yields.

The problem is that the different ways to derive hadronization particle distributions have a profound effect on the resulting fitted temperature. Temperature affects the absolute number of particles through several mechanisms and anticorrelates with the phase space occupancy parameters $`\gamma _i,i=u,d,s`$ ourspspaper ; rafelski2002 . It has been found that the introduction of these parameters, motivated by the need to conserve entropy at hadronization rafelski2002  decrease the $`\chi ^2/`$ per degree of freedom considerably and lowers the freeze-out temperature by 30 MeV rafelski2002 . Other workers assume the light flavors are in chemical equilibrium PBM99 ; becattini ; PBM01 ; PBMRHIC .

Considering ratios of resonances to ground state particles eliminates the fitted temperature’s sensitivity to chemical equilibration, since the numerator and denominator have the same quark composition. In every hadronization model considered here the chemical parameters cancel out and only temperature and fireball freeze-out geometry and dynamics influence the observed ratios.

When fitting the particle spectra, the system’s spatial shape and the way the freeze-out progresses in time have a considerable effect on the form of particle distributions, and hence on the fitted temperature and matter flow. The impact of freeze-out geometry and dynamics on particle spectra were examined well before RHIC data became available fireballspectra ; fireballspectra2  and it was realized that an understanding of freeze-out is essential for the statistical analysis of the fireball heinzkolb . Even though this matter has been clearly recognized, a systematic analysis of how freeze-out geometry affects particle distributions is for the first time attempted here. In fact, each of the models used in the study of particle spectra van-leeuwen ; ourspspaper ; castillo ; burward-hoy ; florkowski  employs a different choice of freeze-out geometry, based on different, often tacitly assumed, hadronization scenarios. Thus an understanding for the influence of hadronization mechanism is impossible to deduce from this diversity.

However, every study of firball hadronization we are aware of uses the Cooper-Frye formula cooperfrye :

$$E\frac{dN}{d^3p}=p^\mu d^3\mathrm{\Sigma }_\mu f(p^\mu u_\mu ,T,\lambda )\theta (p^\mu d^3\mathrm{\Sigma }_\mu ),$$
(1)
where $`p^\mu `$ is the particle’s four-momentum, $`u^\mu `$ is the systems velocity profile, T is the temperature, $`\lambda `$ is a chemical potential, $`f(E,T,\lambda )`$ is the statistical distribution of the emitted particles in terms of energy and conserved quantum numbers and $`\mathrm{\Sigma }^\mu `$ describes the hadronization geometry. It is the covariant generalization of a volume element of the fireball, i.e., a “3D” surface in spacetime from which particles are emitted. $`\theta (p^\mu d^3\mathrm{\Sigma }_\mu )`$ is the step function, which eliminates the possible inward emission cf\_trunc1 ; cf\_trunc2 .

The Cooper-Frye formula, Eq. (1), is believed to be the most general way to implement statistical hadronization emission of particles. For it to represent a physical description of the system, the following two conditions have to be met:
(i) Statistical hadronization must apply. The particles emitted from a volume element (in it’s co-moving frame), will be distributed according to the Bose-Einstein or Fermi-Dirac distributions $`f(E,T,\lambda )`$ for some temperature $`T`$ and fugacity $`\lambda `$.
(ii) A “small” volume element hadronizes rapidly in it’s rest frame, that is, no long lived mixed Quark gluon plasma (QGP) -hadronic confined phase exists.

If this second condition is satisfied, it becomes possible to define a hadronization hypersurface $`\mathrm{\Sigma }^\mu =(t_f(x,y,z),x,y,z))`$ which specifies at which time $`t_f`$ hadrons are emitted from the point $`(x,y,z)`$. In this fast hadronization case differing $`\mathrm{\Sigma }^\mu `$ can be considered for physically differing models. The different choices of $`\mathrm{\Sigma }^\mu `$ correspond to physically different scenarios, and it becomes possible, in principle, to distinguish them experimentally. However, if a long-lived mixed phase does exist, it might well be that the Cooper-Frye formula can be used as an approximation technique to transform a hydrodynamically evolving system into hadrons and authors who worked with a long mixed phase have chosen this approach, see e.g. Ref. shuryak .

### II.2 Freeze-out geometry

The high baryon stopping power observed at SPS energies allNA49 ; piNA49 ; lamNA49  has prompted some authors to use a spherical expansion and freeze-out as an ansatz ourspspaper .

However, at RHIC collision energies the measured $`dN/d\eta `$ phobos ; brahms  indicates that around mid-rapidity the system conditions can be approximated by the Bjorken picture bjorken .

To describe particle spectra measured around midrapidity, therefore, boost invariance becomes the dominant symmetry on which freeze-out geometry should be based. To construct such a hadronization scenario, we consider that the most general cylindrically symmetric flow profile

$$u^\mu =\left(\begin{array}{c}\mathrm{cosh}(y_L)\mathrm{cosh}(y_{})\hfill \\ \mathrm{sinh}(y_{})\mathrm{cos}(\theta )\hfill \\ \mathrm{sinh}(y_{})\mathrm{sin}(\theta )\hfill \\ \mathrm{sinh}(y_L)\mathrm{cosh}(y_{})\hfill \end{array}\right),p^\mu =\left(\begin{array}{c}m_{}\mathrm{cosh}(y)\\ p_{}\mathrm{cos}(\varphi )\\ p_{}\mathrm{sin}(\varphi )\\ m_{}\mathrm{sinh}(y)\end{array}\right)$$
(2)
(the last, longitudinal coordinate is defined along the beam direction) leads to the following rest energy rafrap :

$$p_\mu u^\mu =m_{}\mathrm{cosh}(y_{})\mathrm{cosh}(yy_L)\mathrm{sinh}(y_{})\mathrm{cos}(\theta \varphi )p_{}.$$
(3)
The requirement for the Bjorken picture is that the emission volume element has the same $`y_L`$ dependence:

$$p_\mu d^3\mathrm{\Sigma }^\mu A\mathrm{cosh}(yy_L)+B.$$
(4)
This constrains the freeze-out hypersurface to be of the form

$$\mathrm{\Sigma }^\mu =(t_f\mathrm{cosh}(y_L),x,y,t_f\mathrm{sinh}(y_L)).$$
(5)
Here $`t_f`$ is a parameter invariant under boosts in the z direction, whose physical significance depends on the model considered.

For central collisions, a further simplifying constraint is provided by the cylindrical symmetry, which forces $`t_f`$, as well as $`y_L`$ and $`y_{}`$, to be independent of the angles $`\theta `$ and $`\varphi `$. The freeze-out hypersurface can be parametrized, in this case, as

$$\mathrm{\Sigma }^\mu =(t_f(r)\mathrm{cosh}(y_L),r\mathrm{sin}(\theta ),r\mathrm{cos}(\theta ),t_f(r)\mathrm{sinh}(y_L)),$$
(6)
$`d^3\mathrm{\Sigma }^\mu =t_frdrd\theta dy_L`$
$`(\mathrm{cosh}(y_L){\displaystyle \frac{t_f}{r}}\mathrm{cos}(\theta ),{\displaystyle \frac{t_f}{r}}\mathrm{sin}(\theta ),\mathrm{sinh}(y_L))`$ (7)
And the emission element takes the form

$`p^\mu d^3\mathrm{\Sigma }_\mu =`$
$`\left[m_{}\mathrm{cosh}(yy_L)p_{}{\displaystyle \frac{t_f}{r}}\mathrm{cos}(\theta \varphi )\right]t_frdrd\theta dy_L,`$ (8)
with the same dependence on the angle as Eq. (3). Equation (1) can then be integrated over all the possible values of $`y_L`$ and $`\theta \varphi `$ to give a particle spectrum depending purely on the transverse mass, temperature, and $`y_{}`$. The fits in burward-hoy ; van-leeuwen ; castillo ; florkowski  are based on such an ansatz.

What distinguishes the models currently considered is the time component of the freeze-out surface. The most general freeze-out hypersurface compatible with cylindrical symmetry is provided by Eq.( 6). Generally, $`t_f`$ (a generic function of r) represents the time, in a frame comoving with the longitudinal flow, at which the surface at distance r freezes out.

The fits in Refs. burward-hoy ; van-leeuwen ; castillo  are based on a particular case of such a freeze-out surface, in which $`t_f`$ is completely independent of $`r`$ ($`t_f/r=0`$). Such a picture’s physical reasonableness can be questioned, e.g., why should spatially distant volume elements, presumably with different densities and moving at different transverse velocities, all freeze out simultaneously in a longitudinally comoving frame?. However, such a simple model can perhaps serve as an approximation.

More generally, the “burning log” model fireballspectra ; blastwave  (sometimes referred to as “blast wave”; This term, however, is also used to refer to the $`t_f/r=0`$ model described in the preceding paragraph) assumes that the emission occurs through a three-dimensional hadronization surface which is moving at a constant “velocity” ($`1/v_f=t_f/r`$ throughout the fireball. Both boost-invariant and spherically symmetric versions of burning log model were considered. Even if the hadronization velocity encompasses an extra parameter $`v_f`$, it is worth considering since it is based on a physically motivated hadronization picture. Moreover, the burning log picture is a suitable framework in the study of sudden hadronization. Sudden hadronization occurs when the fireball encounters a mechanical instability sudden , which combined with the fireball’s high transverse flow ensures that the emission surface spreads to the interior of the fireball with $`v_fc`$. All of the indications suggested for such a picture seem to be borne out by both SPS and RHIC data rafelski2002 ; sudden ; pratt .

An approach based on the hypothesis of initial state “synchronization” by the primary instant of collision and the following independent but equivalent evolution of all volume elements assumes that each element of the system undergoes freeze-out at the same proper time $`\tau `$. In this framework each fireball element expands and cools down independently, hadronizing when it’s temperature and density reach the critical value. This model was successfully used to describe RHIC $`m_{}`$-spectra florkowski . In this approach $`t_f`$ in Eq. (6) is equal to $`\tau \mathrm{cosh}(y_{})`$ and the hadronization hypersurface in Eq. (II.2) becomes proportional to the flow vector:

$`\mathrm{\Sigma }^\mu `$ $`=`$ $`\tau u^\mu `$ (9)
$`d^3\mathrm{\Sigma }^\mu `$ $`=`$ $`\tau rdrd\theta dy_Lu^\mu =dVu^\mu `$ (10)
$`r`$ $`=`$ $`\tau \mathrm{sinh}(y_{}).`$ (11)
In this hadronization model the heavy ion fireball behaves similarly to the expanding Hubble universe. In the ‘Hubble’ scenario, the Cooper-Frye formula reduces to the Touscheck Covariant Boltzmann distribution rafrap ; jansbook ; tous1 ; tous2 .

$`{\displaystyle \frac{V_0d^3p}{(2\pi )^3}}e^{E/T}`$ $``$ $`{\displaystyle \frac{V_\mu p^\mu }{(2\pi )^3}}d^4p\mathrm{\hspace{0.17em}2}\delta _0(p^2m^2)e^{p_\mu u^\mu /T}`$ (12)
$`V^\mu `$ $`=`$ $`V_0u^\mu `$ (13)
(Where V is the co moving fireball’s volume element in the local rest frame.)

To summarize and illustrate the diversity of distinct hadronization geometries we present in Table 1 and Fig. 1 the freeze-out scenarios examined here. As we shall see the choice of freeze-out geometry produces in a fit of experimental data a non trivial effect capable of altering significantly the understanding of statistical hadronization parameters.

### II.3 Flow profile

Hydrodynamical expansion of the fireball implies in general that each volume element will have a different density and transverse expansion rate. For this reason, the integral over $`d^3\mathrm{\Sigma }`$ can span a range of flows, weighted by density. In first approximation one can fit data using just an “average” flow velocity throughout the entire fireball NA57 ; ourspspaper :

$`E{\displaystyle \frac{dN}{d^3p}}`$ $`=`$ $`{\displaystyle r𝑑r(Ep_{}\frac{dt_f}{dr})f(T,y_{}(r),\lambda )}`$ (14)
$``$ $`(Ep_{}{\displaystyle \frac{dt_f}{dr}})f(T,y_{},\lambda ).`$
However, if one wants to properly identify $`\frac{dt_f}{dr}`$, the flow profile should be taken into account comingpaper . Hydrodynamic simulations shuryak  accompanied by assumption that freeze-out happens when a volume element reaches a critical energy density indicate that the transverse rapidity will depend linearly with the radius i.e. $`v_{}\mathrm{tanh}(r)`$. This condition, however, is appropriate for a static freeze-out and will not in general hold if the freeze-out is sudden. Other flow profiles have been tried in the literature, arising from dynamical hypothesis. For example, the assumption that the freeze-out occurs at the same time $`t_f`$ results in a quadratic ($`vr^2`$) flow profile hydroheinz , which has also been used recently in fits to data castillo . In the Hubble fireball florkowski  the freeze-out conditions will also result in a distinctive flow profile. Specifically with $`\mathrm{\Sigma }^\mu u^\mu `$, we have $`\gamma vr`$.

Density profiles also depend on the assumed initial condition and the equation of state of the expanding QGP. It has been shown wa98  that different density choices have a considerable effect on both the temperature and flow fits at SPS energies.

Figure. 2 shows how the choice in hadronization dynamics and flow profiles at same given freeze-out temperature and transverse flow can result in a range of inverse spectral slopes. Here the density profiles were assumed to be uniform. It is clear that the same freeze-out parameters give rise to a variety of substantially different particle spectra. Conversely, fits to experimental data will only produce reliable information on the freeze-out conditions if and when we have a prior knowledge of the hadronization geometry and dynamics. Therefore, conclusions about statistical model fits, as well as arguments whether freeze-out occurs simultaneously for different particles or not, cannot be answered while the models used to fit the data are plagued by such uncertainties. We will now turn to study how the measurement of spectra of short-lived resonances might provide us with a way of making progress.

## III Momentum dependence of the resonance-particle ratios as a freeze-out probe

We have shown that the measurement of resonances can probe both the hadronization temperature, and the lifetime of the interacting hadron gas phase sqm2001 ; resonances . Ratios of a generic resonance (henceforward called $`Y^{}`$) to the light particle (which we will refer to as $`Y`$) with an identical number of valence quarks are particularly sensitive probes of freezeout temperature because chemical dependence cancels out within the ratio. If we examine this ratio within a given $`m_{}>m_Y^{}`$ range, we expect to disentangle flow and freeze-out conditions, since the ratio $`Y^{}/Y`$ should not depend on $`m_{}`$ for a purely static and thermal source.

We therefore take the most general Boost-invariant freeze-out hypersurface in the Boltzmann limit (see Table 1. Boost invariant implies this is a good approximation at midrapidity)

$$\frac{dN}{dm_T^2}S(m_{},p_{})=_\mathrm{\Sigma }r𝑑r𝒮(m_{},p_{},r),$$
(15)
where

$`𝒮(m_{},p_{},r)`$ $`=`$ $`m_{}K_1(\beta m_{})I_0(\alpha p_{})`$ (16)
$``$ $`{\displaystyle \frac{t_f}{r}}p_{}K_0(\beta m_{})I_1(\alpha p_{}),`$
with

$`\beta ={\displaystyle \frac{\mathrm{cosh}[y_{}(r)]}{T}},\alpha ={\displaystyle \frac{\mathrm{sinh}[y_{}(r)]}{T}}`$ (17)
and use it to calculate the ratio between two particles with the same chemical composition. The chemical factors cancel out, and we are left with

$$\frac{Y^{}}{Y}=\left(\frac{g^{}}{g}\right)\frac{S(m_{},p_{}^{})}{S(m_{},p_{})},$$
(18)
where $`g^{}`$ and $`g`$ refer to each particle’s degeneracy and the function $`S(m_{},p_{})`$ is given by Eq. (16). (Note that $`m_{}`$ is the same for $`Y^{}`$ and $`Y`$, but $`p_{}`$ varies).

Figure 3 shows the application of this procedure to the cases $`(K^{}+\overline{K^{}})/(K_S)`$ (top), $`\mathrm{\Sigma }^{}(1385)/\mathrm{\Lambda }`$ (middle), and $`\eta ^{}/\eta `$ (bottom) at two freeze-out temperatures and flows: $`T=140\text{MeV},v_{max}/c=0.55`$ on left and $`T=170\text{MeV},v_{max}/c=0.3`$ on the right. Significant deviations from simple constant values are observed, showing the sensitivity of the ratio to freeze-out geometry and dynamics. The analytically simple result in Eq. (18) is valid only if the light particle $`Y`$ has been corrected for feed down from resonances, including $`Y^{}`$. In other words, Eq. (18) as well as Fig. 3 require that decay products from reconstructed $`Y^{}`$ do not appear on the bottom of the ratio. Experiments usually do not do such feed down corrections fachini ; vanburen ; markert ; friese , since this would increase both statistical and systematic error on the ratio, and it is not always possible to do such corrections at all (undetected decays) or in the full range of experimental sensitivity.

Introducing the feed down corrections into Eq. (18), we obtain

$$\frac{Y_{observed}^{}}{Y_{observed}}=\frac{g^{}S(m_{},p_{}^{})}{gS(m_{},p_{})+_ig_i^{}b_{Y_i^{}Y}R(m_{},p_{Ti})}.$$
(19)
Here, $`S(m_{},p_{})`$ describes the directly produced particles and has the form given by Eq. (15) and each term $`R(m_{},p_{Ti}^{})`$ describes a feed down contribution.

In the case of an incoherent many-particle system, such as that we are dealing with, the dynamical (matrix element) part of the decay amplitude factors out respectra , and $`R(m_{},p_{Ti}^{})`$ is obtained by integrating the statistical hadronization distribution with a weight given by the phase space elements of the decay products. Thus, for a generic $`Y^{}Y`$ feed down is given by an $`N`$-body decay,

$`R(m_{},p_{})`$ $`=`$
$`{\displaystyle \underset{j=2}{\overset{N}{}}\frac{d^3p_j}{E_j}S(m_T^{},p_T^{})\delta \left(p_\mu ^{}p_\mu \underset{2}{\overset{N}{}}p_{j\mu }\right)},`$
where the integral is performed over the whole allowed region. If more than one feed down occurs, Eq. (III) can be used iteratively, with the left-hand side to be fed back to the right-hand side at each successive iteration.

In general, this expression can get very complicated, and the Monte Carlo integration becomes necessary. For most cases considered here, where there is one feed down and two or three body decays, Eq. III can be carried out semianalytically florkowski ; resonances ; fireballspectra .

Figure 4 shows the ratios, including feed down of resonances, for the same particles and statistical hadronization conditions as were studied in Fig. 3. In the $`\mathrm{\Sigma }^{}/(`$all $`\mathrm{\Lambda })`$ case we omitted the feed down from $`\mathrm{\Xi }`$ to $`\mathrm{\Lambda }`$ which is usually corrected for (if this is not done the ratio $`\mathrm{\Sigma }^{}/(`$all $`\mathrm{\Lambda })`$ would depend strongly on the chemical potentials). We did allow for the $`\varphi K_SK_L`$ feed down, since it is a strong decay that cannot so easily be corrected for. We note that the feed down from particles with a different chemical composition cannot always be corrected for, and thus some resonances ratios will also acquire a (mild) dependence on the chemical potentials. This is even true for ratios such as $`\eta ^{}/(`$all $`\eta )`$, given different $`s\overline{s}`$ content, in this paper, these type chemical corrections were set equal to unity.

To further study the sensitivity of resonance-particle $`m_{}`$-ratio to freeze-out dynamics, we also present the (feed down corrected) case as a function of $`p_{}`$ rather than $`m_{}`$ in Fig. 5. Unsurprisingly, we see grossly different behaviors, with many of the results coalescing. This of course is an expression of the fact that $`Y^{}`$ and $`Y`$ have dramatically different $`p_{}`$ at the same $`m_{}`$ and vice versa. We believe that the $`m_{}`$ ratio will in general be more sensitive to freeze-out dynamics, since its dependence on $`m_{}`$ is dominantly due to freeze-out geometry and dynamics. However, the $`p_{}`$ dependence seen in Fig. 5 provides an important self-consistency check for our previous results. We have found that the $`m_{}`$ ratios are often greater than unity even though there must be more ground state particles than resonances. Now it can be seen in the $`p_{}`$ ratio, that this requirement is satisfied.

## IV Discussion

In general the the $`m_T`$ and $`p_T`$ dependence of the ratios in Fig 3 and, respectively, Fig 4 depends on freeze-out geometry and dynamics. Changes in temperature and flow velocity alter the shape. The introduction of a steeper flow profile will further raise all of the considered ratios, since a considerable fraction of particles will be produced in regions that do not flow as much. The effect of freeze-out dynamics will generally go in the same direction as freeze-out approaches the explosive limit ($`dt_f/dr1`$). However, both the magnitude and the qualitative features of the two effects (flow and freeze-out velocity) will be considerably different. Especially, when more than one ratio is measured, it would appear that we will be able to determine the freeze-out condition. This is in contrast to the $`m_{}`$ distributions in Fig. 2, where the effects discussed in this paper result in linear corrections, which tend to compete, making the task of extracting the freeze-out dynamics much more ambiguous. Thus, there is considerable potential of resonance-particle $`m_{}`$-ratios as a freeze-out probe.

The presence of a long living hadronic gas rescattering phase can distort our freeze-out probe. In particular, the apparent $`Y^{}/Y`$ ratio will be altered due to the depletion of the detectable resonances through the rescattering of their decay products. It’s dependence on $`m_{}`$ will be affected in a non-trivial way, since faster (higher $`p_{}`$) resonances will have a greater chance to escape the fireball without decaying, thus avoiding the rescattering phase altogether. Regeneration of resonances in hadron scattering may add another $`m_{}`$ dependence which is different for the $`\mathrm{\Sigma }^{}/\mathrm{\Lambda }`$ and the $`K^{}/K`$ ratios bleicher . Other signals of the existence of such an interacting hadron gas phase have been considered sqm2001 ; resonances . Fortunately, there is no evidence that a rescattering phase plays a great role in particle distributions. Even so, it would seem that the “safest” probes for freeze-out are the particles and resonances most unlikely to rescatter.

For this reason we have included the $`\eta ^{}/\eta `$ ratio in our considerations. $`\eta \gamma \gamma `$ and $`\eta ^{}\gamma \gamma `$ have very different branching ratios, but have the same degeneracies and similar but rather small partial widths. The electromagnetic decay mode is practically insensitive to posthadronization dynamics. Regeneration effects are suppressed since the hadronic two body decay channel is suppressed. All these features make these particles interesting probes, allowing for the analysis considered here. $`\eta ,\eta ^{}`$ mesons have been measured at SPS energies in the $`\gamma \gamma `$ decay channel eta1 ; eta2 , and detectors such as PHENIX are capable of reconstructing the same decays at RHIC.

While a long rescattering phase would affect the $`\mathrm{\Sigma }^{}`$ distribution, the effect would be very easy to detect experimentally: $`95\%`$ of $`\mathrm{\Sigma }^{}`$ decay through the $`p`$-wave $`\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi `$ channel. However, regenerating $`\mathrm{\Sigma }^{}`$ s in a gas of $`\mathrm{\Lambda }`$ s and $`\pi `$ s is considerably more difficoult, since $`\mathrm{\Lambda }\pi `$ scattering will be dominated by the s-wave $`\mathrm{\Lambda }\pi \mathrm{\Sigma }^\pm `$. This situation will not occur for $`K^{}K\pi `$, since both decay and regeneration happen through the same process, leading to a very fast reequilibration time bleicher . Since both $`\mathrm{\Sigma }^{}/\mathrm{\Lambda }`$ and $`K^{}/K`$ ratios have been calculated within the thermal model sqm2001  (neglecting rescattering), a strongly depleted $`\mathrm{\Sigma }^{}/\mathrm{\Lambda }`$ ratio (compared with $`K^{}/K`$ ) suggests that a statistical freeze-out description, such as that given in this paper, is incomplete, and an interacting hadron gas phase is also necessary.

In summary, we have presented an overview of the different statistical freeze-out models used to fit heavy ion data. We have shown how the freeze-out geometry and freeze-out dynamics influences the hadron spectra. Our primary result is the finding that the $`m_{}`$ dependence of the resonance-particle ratios is a probe of freeze-out. We have presented these ratios for three particle species and two freeze-out conditions and have considered how our results could be altered by posthadronization phenomena.

###### Acknowledgements.

Supported by a grant from the U.S. Department of Energy, DE-FG03-95ER40937 . We thank Patricia Fachini (BNL) and Zhangbu Xu (Yale), from the STAR collaboration, as well as Marcus Bleicher (ITP-Frankfurt) for fruitful discussions.


text:
# Adaptive Methods for PDE’s Wavelets or Mesh Refinement?

## 1. Introduction

Among those relevant phenomenons which are modelled by partial differential or integral equations, countless are the instances where the mathematical solutions exhibit singularities. Perhaps the most classical examples are elliptic equations on domains with re-entrant corners, or nonlinear hyperbolic systems of conservation laws. While such singularities are sources of obvious theoretical difficulties—classical solutions should be abandonned to the profit of weak solutions—they are also an obstacle to the convergence of numerical approximation methods, in the sense that they deteriorate the rate of decay of the error with respect to the size of the discrete problem : achieving a prescribed accuracy will typically require finer resolution and therefore heavier computational cost and memory storage, in comparison to the approximation of smooth solutions. Let us remark that singularities often have a physical relevance : they represent the concentration of stress in elasticity, boundary layers in viscous fluid flows, shock waves in gas dynamics… It is therefore a legitimous requirement that they should be accurately resolved by the numerical method.

In this context, the use of adaptive methods, appears as a natural solution to improve the approximation at a reasonable computational cost. Here, the word adaptivity has a twofold meaning : (i) the discretization is allowed to be refined only locally, in particular near the singularities of the solution, and (ii) the resolution algorithm uses information gained during a given stage of the computation in order to derive a new refined discretization for the next stage. The most typical example is adaptive mesh refinement based on a-posteriori error estimates in the finite element context. While these methods have proved to be computationally successful, the theory describing their advantages over their non-adaptive counterpart is far from being complete. In particular, the rate of convergence of the adaptive algorithm, which describes the trade-off between the accuracy and complexity of the approximation, is not clearly understood.

In recent years, wavelet bases have been proposed as an alternative to adaptive mesh refinement, motivated by their good performances in data (more specifically image) compression. In wavelet-based adaptive schemes, the set of basis functions which describe the approximate solution is updated at each stage of the computation. Intuitively, the selection of the appropriate basis functions plays a similar role as the selection of the mesh points in adaptive finite element methods, and one could therefore expect similar performances from both approaches. On a more rigorous level, a specific feature of the wavelet approach is the emergence of a sound theoretical setting which allows to tackle fundational questions such as the rate of convergence of the adaptive method.

The goal of this paper is to give some elements of comparison between adaptive wavelet and mesh refinement methods from this perspective. We shall first describe in §2 a general setting which leads us in §3 to a first comparison between wavelets and adaptive finite elements from the point of view of approximation theory. We discuss in §4 the relation between these results and adaptive algorithms for PDE’s. After recalling in §5 the classical approach in the finite element context, we present in §6 an adaptive wavelet strategy which has been applied to various problems, and discuss its fundational specificities. Finally, we shall conclude in §7 by pointing out some intrinsic shortcoming of wavelet-based adaptive methods.

## 2. A general framework

Approximation theory is the branch of mathematics which studies the process of approximating general functions by simple functions such as polynomials, finite elements or Fourier series. It plays therefore a central role in the accuracy analysis of numerical methods. Numerous problems of approximation theory have in common the following general setting : we are given a family of subspaces $`(S_N)_{N0}`$ of a normed space $`X`$, and for $`fX`$, we consider the best approximation error

$$\sigma _N(f):=\underset{gS_N}{inf}fg_X.$$
(1)
Typically, $`N`$ represents the number of parameters which are needed to describe an element in $`S_N`$, and in most cases of interest, $`\sigma _N(f)`$ goes to zero as this number tends to infinity. If in addition $`\sigma _N(f)CN^s`$ for some $`s>0`$, we say that $`f`$ is approximated at rate $`s`$.

Given such a setting, the central problem of approximation theory is to characterize by some analytic (typically smoothness) condition those functions $`f`$ which are approximated at some prescribed rate $`s>0`$. Another important problem is how to design simple approximation procedures $`ff_N\mathrm{\Sigma }_N`$ which avoid solving the minimization problem (1), while remaining near optimal in the sense that

$$ff_N_XC\sigma _N(f),$$
(2)
for some constant $`C`$ independent of $`N`$ and $`f`$.

As an example, consider approximation by finite element spaces $`V_h`$ defined from regular conforming partitions $`𝒯_h`$ of a domain $`\mathrm{\Omega }\text{I}\text{R}^\mathrm{d}`$ into simplices with uniform mesh size $`h`$. The approximation theory for such spaces is quite classical, see e.g. , and can be summarized in the following way. If $`W^{t,p}`$ denotes the classical Sobolev space, consisting of those functions $`fL^p`$ such that $`D^\alpha fL^p`$ for $`|\alpha |t`$, we typically have

$$fW^{t+r,p}\underset{gV_h}{inf}fg_{W^{t,p}}Ch^r$$
(3)
provided that $`V_h`$ is contained in $`W^{t,p}`$ and that $`V_h`$ has approximation order larger than $`t+r`$, i.e. contains all polynomials of degree strictly less than $`t+r`$. Such classical results also hold for fractional smoothness. We can express them in terms of the number of parameters, remarking that $`N:=\mathrm{dim}(V_h)h^d`$, so that if we set $`X=W^{t,p}`$ and $`S_N:=V_h`$ with $`h:=N^{1/d}`$, we have obtained

$$fW^{t+r,p}\sigma _N(f)CN^{r/d}.$$
(4)
We have thus identified an analytic condition which ensures the rate $`s=r/d`$. Note that this is not a characterization (we only have an implication), yet a deeper analysis shows that an “if and only if” result holds if we slightly modify the notion of Sobolev smoothness (using Besov classes, see ). In summary, the rate of approximation in $`W^{s,p}`$ is governed by the approximation order of the $`V_h`$ spaces, the dimension $`d`$ and the level of smoothness of $`f`$ measured in $`L^p`$. Let us finally remark that near-optimal approximation procedures can be obtained if we can find a sequence of finite element projectors $`P_N:XS_N`$ such that $`P_N_{XX}K`$ with $`K`$ independent of $`N`$ : in this case, we simply take $`f_N=P_Nf`$ and remark that $`ff_N_X(1+K)\sigma _N(f)`$.

In the following we shall address the same questions in the cases of adaptive finite element and wavelet approximation. As we shall see, a specific feature to such cases is that the spaces $`\mathrm{\Sigma }_N`$ are not linear vector spaces.

## 3. Adaptive finite elements and wavelets

In the adaptive finite element setting, the number of parameters $`N`$ is proportional to the number of triangles, but for a given budget $`N`$ the partition $`𝒯`$ and the finite element space $`V_𝒯`$ are allowed to be locally refined in a way which depends on the function $`f`$ to be approximated. It is therefore natural to define the approximation spaces $`S_N`$ as

$$S_N:=_{\mathrm{\#}(𝒯)N}V_𝒯.$$
(5)
It should be well understood that the $`S_N`$ are not linear vector spaces (the sum of two elements does not in general fall in $`S_N`$ when their triangulation do not match) but any $`gS_N`$ is still described by $`𝒪(N)`$ parameters, which encode both its triangulation $`𝒯`$ and its coordinates in $`V_𝒯`$. The requirement of adaptivity has thus led us to the concept of nonlinear approximation.

Wavelet bases offer another track toward nonlinear adaptive approximation. The simplest prototype of a wavelet basis is the Haar system. Let us describe this system in the case of expanding a function $`f`$ defined on $`[0,1]`$ : the first component in this expansion is simply the average of $`f`$, i.e. the orthogonal projection $`f,e_0e_0`$ onto the function $`e_0=\chi _{[0,1]}`$. The approximation is then refined into the average of $`f`$ on the two half intervals of equal size. This refinement amounts in adding the orthogonal projection $`f,e_1e_1`$ onto the function $`e_1=\chi _{[0,1/2]}\chi _{[1/2,1]}`$. Iterating this refinement process, we see that the next components have the same form as $`e_1`$ up to a change of scale : at refinement level $`j`$, we are adding the projection onto the functions

$$\psi _{j,k}(x)=2^{j/2}\psi (2^jxk),k=0,\mathrm{},2^j1,$$
(6)
where $`\psi =e_1`$. Since all these functions are orthogonal to the previous ones, letting $`j`$ go to $`+\mathrm{}`$, we obtain the expansion of $`f`$ into an orthonormal basis of $`L^2([0,1])`$

$$f=\underset{\lambda }{}f_\lambda \psi _\lambda ,$$
(7)
with $`f_\lambda :=f,\psi _\lambda `$. In the above notation $`\lambda `$ concatenates the scale and space parameters $`j`$ and $`k`$, and $``$ is the set of all indices (including also the first function $`e_0`$). In order to keep track of the scale $`j`$ corresponding to an index $`\lambda =(j,k)`$ we shall use the notation $`|\lambda |=j`$. More general wavelet systems in one or several space dimension are built from similar nested approximation processes, involving e.g. spline functions or finite elements in place of piecewise constant functions (see  or  for a general presentation).

This brief description suggests that a natural construction of adaptive wavelet approximations is obtained by using only a limited set of indices $`\lambda `$ as the scale $`|\lambda |`$ grows, which depends on the function to be approximated and typically corresponds to those wavelets whose supports are close to its singularities. It is therefore natural to define the approximation spaces $`S_N`$ as the set of all $`N`$ terms combinations

$$S_N:=\{\underset{\lambda \mathrm{\Lambda }}{}d_\lambda \psi _\lambda ;\mathrm{\#}(\mathrm{\Lambda })N\}.$$
(8)
Again this is obviously not a linear space, since we allow to approximate a function by choosing the best $`N`$ terms which differ from one function to another. Note that we still do have $`S_N+S_N=S_{2N}`$.

Both adaptive finite element and wavelet framework have obvious similiarities. However, the answer to the two basic questions raised in the previous section—what are the properties of $`f`$ which govern the decay of $`\sigma _N(f)`$ and how to compute in a simple way a near optimal approximation of $`f`$ in $`\mathrm{\Sigma }_N`$ —is only fully understood in the wavelet framework. Concerning the first question, a striking result by DeVore and his collaborators  is the following : with $`X:=W^{t,p}`$, best $`N`$-term wavelet approximation satisfies

$$fW^{t+r,q}\sigma _N(f)CN^{r/d},$$
(9)
with $`q`$ and $`r`$ connected by the relation $`1/q=1/p+r/d`$, assuming that the multiresolution approximation spaces associated to the wavelet basis are in $`W^{t,p}`$ and contain the polynomials of degree strictly less than $`t+r`$.

Such an estimate should be compared with the linear estimate (4) : the same convergence rate is governed by a much weaker smoothness assumption on $`f`$ since $`q<p`$ (as in the linear case, an “iff and only if” result can be obtained up to slight technical modifications in the statement of (9)). This result gives a precise mathematical meaning to the spatial adaptation properties of best $`N`$-term wavelet approximation: a function $`f`$ having isolated discontinuity, has usually a smaller amount of smoothness $`s+t`$ when measured in $`L^p`$ than when measured in $`L^q`$ with $`1/q=1/p+t/d`$, and therefore $`\sigma _N(f)`$ might decrease significantly faster than $`\epsilon _N(f)`$.

The answer to the second question is given by a result due to Temlyakov : if $`f=_\lambda d_\lambda \psi _\lambda `$, and if we measure the approximation error in $`X=W^{t,p}`$, a near optimal strategy when $`1<p<\mathrm{}`$ consists in the thresholding procedure which retains the $`N`$ largest contributions $`d_\lambda \psi _\lambda _X`$ : if $`\mathrm{\Lambda }_N`$ is the corresponding set of indices, one can prove that there exists $`C>0`$ independent of $`N`$ and $`f`$ such that

$$f\underset{\lambda \mathrm{\Lambda }_N}{}d_\lambda \psi _\lambda _XC\sigma _N(f).$$
(10)
This fact is obvious when $`X=L^2`$ using the orthonormal basis property. It is a remarkable property of wavelet bases that it also holds for more general function spaces. In summary, thresholding plays for best $`N`$-term wavelet approximation an analogous role as projection for linear finite element approximation.

In the adaptive finite element framework, a similar theory is far from being complete. Partial answers to the basic questions are available if one chooses to consider adaptive partitions with shape constraints in terms of a uniform bound on the aspect ratio of the elements

$$\underset{K𝒯}{\mathrm{max}}\left([\mathrm{Diam}(K)]^d/\mathrm{vol}(K)\right)C.$$
(11)
Such a restriction means that the local refinement is isotropic, in a similar way to wavelets. In such a case, we therefore expect a rate of approximation similar to (9). Such a result is not available, yet the following can be proved  for Lagrange finite elements of degree $`m`$ : if for any given tolerance $`\epsilon >0`$, one is able to build a partition $`𝒯=𝒯(\epsilon )`$ of cardinality $`N=N(\epsilon )`$ such that on each $`K𝒯`$ the local error of approximation by polynomials satisfies

$$\frac{\epsilon }{2}\underset{p\mathrm{\Pi }_m}{inf}fp_{W^{t,p}(K)}\epsilon ,$$
(12)
then we can build global approximants $`f_NV_𝒯S_N`$ such that

$$fW^{t+r,q}ff_N_XCN^{r/d},$$
(13)
with $`q`$ and $`r`$ connected by the relation $`1/q=1/p+r/d`$ and assuming $`s+t<m`$. The effective construction of $`𝒯(\epsilon )`$ is not always feasible, in particular due to the conformity constraints on the partition which does not allow to connect very coarse and very fine elements without intermediate grading. However, this result shows that from an intuitive point of view, the adaptive finite element counterpart to wavelet thresholding amounts in equilibrating the local error over the partition. One can actually use these ideas in order to obtain the estimate (9) for adaptive finite elements under the more restrictive assumption that $`1/q<1/p+r/d`$. Let us finally mention that the approximation theory for adaptive finite elements without shape constraints is an open problem.

## 4. Nonlinear approximation and PDE’s

Nonlinear approximation theory has opened new lines of research on the theory of PDE’s and their numerical discretization. On the one hand, it is worth revisiting the regularity theory of certain PDE’s for which the solutions develop singularities but might possess significantly higher smoothness in the scale of function spaces which govern the rate of nonlinear approximation in a given norm than in the scale which govern the rate of linear approximation in the same norm. Results of this type have been proved in particular for elliptic problems on nonsmooth domains  and for scalar 1D conservation laws . These results show that if $`u`$ is the solution of such equations, the rate of decay of $`\sigma _N(u)`$ is significantly higher for best $`N`$-term approximation than for the projection on uniform finite element spaces, therefore advocating for the use of adaptive discretizations of such PDE’s.

On the other hand these results also provide with an ideal benchmark for adaptive discretizations of the equation, since $`\sigma _N(u)`$ represents the best accuracy which can be achieved by $`N`$ parameters. In the wavelet case these parameters are typically the $`N`$ largest coefficients of the exact solution $`u`$. However, in the practice of solving a PDE, these coefficients are not known, and neither is the set $`\mathrm{\Lambda }`$ corresponding to the indices of the $`N`$ largest contributions $`d_\lambda \psi _\lambda `$. It is therefore needed to develop appropriate adaptive resolution strategies as a substitute to the thresholding procedure. Such strategies aim at detecting the indices of the largest coefficients of the solutions and to compute them accurately, in a similar way that adaptive mesh refinement strategies aim at contructing the optimal mesh for finite element approximation. In both contexts, we could hope for an algorithm which builds approximations $`u_N\mathrm{\Sigma }_N`$ such that $`uu_N_X`$ is bounded up to a fixed multiplicative constant by $`\sigma _N(u)`$ for a given norm of interest, but this requirement is so far out of reach. A more reasonable goal is that the adaptive strategy exhibits the optimal rate of approximation : if $`\sigma _N(u)CN^s`$ for some $`s>0`$, then $`uu_N_XCN^s`$ up to a change in the constant. Another requirement is that the adaptive algorithm should be scalable, i.e. the number of elementary operations in order to compute $`u_N`$ remains proportional to $`N`$. Let us finally remark that the norm $`_X`$ for which error estimates can be obtained is often dictated by the nature of the equation (for example $`X=H^1`$ in the case of a second order elliptic problem) and that additional difficulties can be expected if one searches for estimates in a different norm.

## 5. The classical approach

The classical approach to numerically solving linear and nonlinear partial differential or integral equations $`(u)=0`$ by the finite element method is typically concerned with the following issues :

* Well-posedness of the equation, i.e. existence, uniqueness and stability of the solution.
* Discretization into a finite element problem $`_𝒯(u_𝒯)=0`$ by the Galerkin method with $`u_𝒯V_𝒯`$, analysis of well-posedness and of the approximation error $`uu_𝒯_X`$.
* Numerical resolution of the finite dimensional system.
* Mesh refinement based on a-posteriori error estimators in the case of adaptive finite element methods.

Several difficulties are associated to each of these steps. First of all, note that the well-posedness of the finite element problem is in general not a consequence of the well-posedness of the continuous problem. Typical examples even in the linear case are saddle point problems. For such problems, it is well known that, for Galerkin discretizations to be stable, the finite element spaces for the different solution components have to satisfy certain compatibility conditions (LBB or Inf-Sup condition), which are also crucial in the derivation of optimal error estimates. Thus the discrete problem does not necessarily inherit the “nice properties” of the original infinite dimensional problem. Concerning the numerical resolution of the discrete system, a typical source of trouble is its possible ill-conditioning, which interferes with the typical need to resort on iterative solvers in high dimension. An additional difficulty occuring in the case of integral equations is the manipulation of matrices which are densely populated.

Finally, let us elaborate more on the adaptivity step. Since more than two decades, the understanding and practical realization of adaptive refinement schemes in a finite element context has been documented in numerous publications . Key ingredients in most adaptive algorithms are a-posteriori error estimators which are typically derived from the current residual $`(u_𝒯)`$ : in the case where the Frechet derivative $`D(u)`$ is an isomorphism between Banach function spaces $`X`$ to $`Y`$, one can hope to estimate the error $`uu_𝒯_X`$ by the evaluation of $`(u_𝒯)_Y`$. The rule of thumb is then to decompose $`(u_𝒯)_Y`$ into computable local error indicators $`\eta _K`$ which aim to describe as accurately as possible the local error on each element $`K𝒯`$. In the case of elliptic problems, these indicators typically consist of local residuals and other quantities such as jumps of derivatives across the interface between adjacent elements. A typical refinement algorithm will subdivide those elements $`K`$ for which the error indicator $`\eta _K`$ is larger than a prescribed tolerance $`\epsilon `$ resulting in a new mesh $`\stackrel{~}{𝒯}`$. Note that this strategy is theoretically in accordance with our remarks in §3 on adaptive finite element approximation, since it tends to equilibrate the local error. Two other frequently used strategies consist in refining a fixed proportion of the elements corresponding to the largest $`\eta _K`$, or the smallest number of elements $`K`$ for which the $`\eta _K`$ contribute to the global error up to a fixed proportion. It is therefore hoped that the iteration of this process from an initial mesh $`𝒯_0`$ will produce optimal meshes $`(𝒯_n)_{n0}`$ in the sense that the associated solutions $`u_n:=u_{𝒯_n}V_{𝒯_n}`$ converge to $`u`$ at the optimal rate :

$$\sigma _N(u)CN^ruu_n_XC[\mathrm{\#}(𝒯_n)]^r,$$
(14)
up to a change in the constant $`C`$. Unfortunately, severe obstructions appear when trying to prove (14) even in the simplest model situations. One of them is that $`\eta _K`$ is in general not an estimate by above of the local error, reducing the chances to derive the optimal rate. For most adaptive refinement algorithms, the theoretical situation is actually even worse in the sense that it cannot even be proved that the refinement step actually results in a reduction of the error by a fixed amount and that $`u_n`$ converges to $`u`$ as $`n`$ grows. Only recently  have proof of convergence appeared for certain type of adaptive finite element methods, yet without convergence rate and therefore no guaranteed advantage over their non-adaptive counterparts.

## 6. A new paradigm

Wavelet methods vary from finite element method in that they can be viewed as solving systems that are finite sections of one fixed infinite dimensional system corresponding to the discretization of the equation in the full basis. This observation has led to a new paradigm which has been explored in  for linear variational problems. It aims at closely intertwining the analysis—discretization—solution process. The basic steps there read as follows :

* Well-posedness of the variational problem.
* Discretization into an equivalent infinite dimensional problem which is well posed in $`\mathrm{}^2`$.
* Devise an iterative scheme for the $`\mathrm{}_2`$-problem that exhibits a fixed error reduction per iteration step.
* Numerical realization of the iterative scheme by means of an adaptive application of the involved infinite dimensional operators within some dynamically updated accuracy tolerances.

Thus the starting point (n1) is the same. The main difference is that one aims at staying as long as possible with the infinite dimensional problem. Only at the very end, when it comes to applying the operators in the ideal iteration scheme (n4), one enters the finite dimensional realm. However, the finite number of degrees of freedom is determined at each stage by the adaptive application of the operator, so that at no stage any specific trial space is fixed.

The simplest example is provided by the Poisson equation $`\mathrm{\Delta }u=f`$ on a domain $`\mathrm{\Omega }`$ with homogeneous boundary conditions, for which the variational formulation in $`X=H_0^1`$ reads : find $`uX`$ such that

$$a(u,v)=L(v),\text{for all}vX,$$
(15)
with $`a(u,v):=_\mathrm{\Omega }uv`$ and $`L(v):=_\mathrm{\Omega }fv`$. The well-posedness for a data $`fX^{}=H^1`$ is ensured by the Lax-Milgram lemma. In the analysis of the wavelet discretization of this problem, we shall invoke the fact that wavelet bases provide norm equivalence for Sobolev spaces in terms of weighted $`\mathrm{}^2`$ norms of the coefficients : if $`u=_\lambda u_\lambda \psi _\lambda `$, one has

$$u_{H^s}^2\underset{\lambda }{}u_\lambda \psi _\lambda _{H^s}^2\underset{\lambda }{}2^{2s|\lambda |}|u_\lambda |^2.$$
(16)
We refer to  and  for the general mechanism allowing to derive these equivalences, in particular for Sobolev spaces on domains with boundary conditions such as $`H_0^1`$. Therefore, if we renormalize our system in such a way that $`\psi _\lambda _X=1`$, we obtain the norm equivalence

$$u_X^2U^2,$$
(17)
where $`U:=(u_\lambda )_\lambda `$ and $``$ denotes the $`\mathrm{}^2`$ norm. By duality, one also easily obtains

$$f_X^{}^2F^2,$$
(18)
with $`F:=(f,\psi _\lambda )_\lambda `$. The equivalent $`\mathrm{}^2`$ system is thus given by

$$AU=F,$$
(19)
where $`A(\lambda ,\mu )=a(\psi _\lambda ,\psi _\mu )`$ is a symmetric positive definite matrix which is continuous and coercive in $`\mathrm{}^2`$. In this case, a converging infinite dimensional algorithm can simply be obtained by the Richardson iteration

$$U^n:=U^{n1}+\tau (FAU^{n1})$$
(20)
with $`0<\tau <2[\lambda _{\mathrm{max}}(A)]^1`$ and $`U^0=0`$, which guarantees the reduction rate $`UU^n\rho UU^{n1}`$ with $`\rho =\mathrm{max}\{1\tau \lambda _{\mathrm{min}}(A),\tau \lambda _{\mathrm{max}}(A)1\}`$. Note that renormalizing the wavelet system plays the role of a multiscale preconditioning, similar to multigrid yet operated at the infinite dimensional level.

At this stage, one enters finite dimensional adaptive computation by modifying the Richardson iteration up to a prescribed tolerance according to

$$U^n:=U^{n1}+\tau (\mathrm{𝐂𝐎𝐀𝐑𝐒𝐄}(F,\epsilon )\mathrm{𝐀𝐏𝐏𝐑𝐎𝐗}(AU^{n1},\epsilon ))$$
(21)
where $`F\mathrm{𝐂𝐎𝐀𝐑𝐒𝐄}(F,\epsilon )\epsilon `$ and $`AU\mathrm{𝐀𝐏𝐏𝐑𝐎𝐗}(AU^{n1},\epsilon )\epsilon `$, and the $`U^n`$ are now finite dimensional vector supported by adaptive sets of indices $`\mathrm{\Lambda }_n`$. The procedure COARSE, which simply corresponds to thresholding the data vector $`F`$ at a level corresponding to accuracy $`\epsilon `$, can be practically achieved without the full knowledge of $`F`$ by using some a-priori bounds on the size of the coefficients $`f,\psi _\lambda `$, exploiting the local smoothness of $`f`$ and the oscillation properties of the wavelets. The procedure APPROX deserves more attention. In order to limitate the spreading effect of the matrix $`A`$, one invokes its compressibility properties, namely the possibily to truncate it into a matrix $`A_N`$ with $`N`$ non-zero entries per rows and columns in such a way that

$$AA_N_{\mathrm{}^2\mathrm{}^2}CN^s.$$
(22)
The rate of compressibility $`s`$ depends on the available a-priori estimates on the off-diagonal entries $`A(\lambda ,\mu ):=_\mathrm{\Omega }\psi _\lambda \psi _\mu `$ which are consequences of the smoothness and vanishing moment properties of the wavelet system, see . Once these properties are established, a first possibility is thus to choose

$$\mathrm{𝐀𝐏𝐏𝐑𝐎𝐗}(AU^{n1},\epsilon )=A_NU^{n1}$$
(23)
with $`N`$ large enough so that accuracy $`\epsilon `$ is ensured. Clearly the modified iteration (21) satisfies $`UU^n\rho UU^{n1}+2\tau \epsilon `$, and therefore ensures a fixed reduction rate until the error is of the order $`\frac{2\tau }{1\rho }\epsilon `$, or until the residual $`FAU^n`$ is of order $`\frac{2\tau A}{1\rho }\epsilon `$. A natural idea is therefore to update dynamically the tolerance $`\epsilon `$, which is first set to $`1`$ and divided by $`2`$ each time the approximate residual $`\mathrm{𝐂𝐎𝐀𝐑𝐒𝐄}(F,\epsilon )\mathrm{𝐀𝐏𝐏𝐑𝐎𝐗}(AU^{n1},\epsilon )`$ is below $`[\frac{2\tau A}{1\rho }+3]\epsilon `$ (which is ensured to happen after a fixed number of steps).

We therefore obtain a converging adaptive strategy, so far without information about the convergence rate. It turns out that the optimal convergence rate can also be proved, with a more careful tuning of the adaptive algorithm. Two additional ingredients are involved in this tuning.

Firstly, the adaptive matrix vector multiplication APPROX has to be designed in a more elaborate way than (23) which could have the effect of inflating too much the sets $`\mathrm{\Lambda }_n`$. Instead, one defines for a finite length vector $`V`$

$$\mathrm{𝐀𝐏𝐏𝐑𝐎𝐗}(AV,\epsilon )=\underset{l=0}{\overset{j}{}}A_{2^{jl}}[V_{2^l}V_{2^{l1}}]$$
(24)
where $`V_N`$ denotes the restriction of $`V`$ to its $`N`$ largest components (with the notation $`V_{1/2}=0`$), and $`j`$ is the smallest positive integer such that the residual $`_{l=0}^jAA_{2^{jl}}V_{2^l}V_{2^{l1}}+AVV_{2^j}`$ is less than $`\epsilon `$. In this procedure, the spreading of the operator is more important on the largest coefficients which are less in number, resulting in a significant gain in the complexity of the outcome.

Secondly, additional coarsening steps are needed in order to further limitate the spreading of the sets $`\mathrm{\Lambda }_n`$ and preserve the optimal rate of convergence. More precisely, the procedure COARSE is applied to $`U^n`$ with a tolerance proportional to $`\epsilon `$, for those $`n`$ such that $`\epsilon `$ will be updated at the next iteration.

With such additional ingredients, it was proved in  that the error has the optimal rate of decay in the sense that

$$\sigma _N(u)CN^suu_n_XUU^nC[\mathrm{\#}(\mathrm{\Lambda }_n)]^s,$$
(25)
and that moreover, the computational cost of producing $`u_n`$ remains proportional to $`\mathrm{\#}(\mathrm{\Lambda }_n)`$. It is interesting to note that this strategy extends to non-elliptic problems such as saddle-points problems, without the need for compatibility conditions, since one inherits the well-posedness of the continuous problem which allows to obtain a converging infinite dimensional iteration, such the Uzawa algorithm or a gradient descent applied to the least-square system (see also ). The extension to nonlinear variational problems, based on infinite dimensional relaxation or Newton iterations, has also been considered in . It requires a specific procedure for the application of the nonlinear operator in the wavelet coefficients domain which generalizes (23). It should also be mentioned that matrix compressibility also applies in the case of integral operators which have quasi-sparse wavelet discretizations. Therefore several of the obstructions from the classical approach—conditioning, compatibility, dense matrices—have disappeared in the wavelet approach.

Let us finally mention that the coarsening steps are not really needed in the practical implementations of the adaptive wavelet method (for those problems which have been considered so far) which still does exhibit optimal convergence rate. However, we do not know how to prove (25) without these coarsening steps. There seems to be a similar situation in the finite element context : it has recently been proved in  that (14) can be achieved by an adaptive mesh refinement algorithm which incorporates coarsening steps, while these steps are not needed in practice.

## 7. Conclusions and shortcomings

There exist other approaches for the development of efficient wavelet-based adaptive schemes. In particular, an substantial research activity has recently been devoted to multiresolution adaptive processing techniques, following the line of idea introduced in . In this approach, one starts from a classical and reliable scheme on a uniform grid (finite element, finite difference or finite volume) and applies a discrete multiresolution decomposition to the numerical data in order to compress the computational time and memory space while preserving the accuracy of the initial scheme. Here the adaptive sets $`\mathrm{\Lambda }_n`$ are therefore limited within the resolution level of the uniform grid where the classical scheme operates. This approach seems more appropriate for hyperbolic initial value problems , in which a straightforward wavelet discretization might fail to converge properly. It should again be compared to its adaptive mesh refinement counterpart such as in .

Let us conclude by saying that despite its theoretical success, in the sense of achieving for certain classes of problems the optimal convergence rate with respect to the number of degrees of freedom, the wavelet-based approach to adaptive numerical simulation suffers from three major curses.
The curse of geometry : while the construction of wavelet bases on a rectangular domains is fairly simple—one can use tensor product techniques and inherit the simplicity of the univariate construction—it is by far less trivial for domains with complicated geometries. Several approaches have been proposed to deal with this situation, in particular domain decomposition into rectangular patches or hierarchical finite element spaces, see , and concrete implementations are nowaday available, but they result in an unavoidable loss of structural simplicity in comparison to the basic Haar system of §3.
The curse of data structure : encoding and manipulating the adaptive wavelet approximations $`U^n`$ to the solution means that we both store the coefficients and the indices of the adaptive set $`\mathrm{\Lambda }_n`$ which should be dynamically updated. The same goes for the indices of the matrix $`A`$ which are used in the matrix-vector algorithm (24) at each step of the algorithm. This dynamical adaptation, which requires appropriate data structure, results in major overheads in the computational cost which are observed in practice : the numerical results in  reveal that while the wavelet adaptive algorithm indeed exhibits the optimal rate of convergence and slightly outperforms adaptive finite element algorithms from this perspective, the latter remains significantly more efficient from the point of view of computational time.
The curse of anisotropy : adaptive wavelet approximation has roughly speaking the same properties as isotropic refinement. However, many instances of singularities such as boundary layers and shock waves, have anisotropic features which suggests that the refinement should be more pronounced in one particular direction. From a theoretical point of view, the following example illustrate the weakness of wavelet bases in this situation : if $`f=\chi _\mathrm{\Omega }`$ with $`\mathrm{\Omega }\text{I}\text{R}^\mathrm{d}`$ a smooth domain, then the rate of best $`N`$-term approximation in $`X=L^2`$ is limited to $`r=1/(2d2)`$ and therefore deteriorates as the dimension grows. Wavelet bases should therefore be reconsidered if one wants to obtain better rates which take some advantage of the geometric smoothness of the curves of dicsontinuities. On the adaptive finite element side, anisotropic refinement has been considered and practically implemented, yet without a clean theory available for the design of an optimal mesh.

The significance of wavelets in numerical analysis remains therefore tied to these curses and future breakthrough are to be expected once simple and appropriate solutions are proposed in order to deal with them.