text:
# Some remarks on 𝐴₁⁽¹⁾ soliton cellular automata

## 1 Introduction

Several years ago, Torii, Takahashi and Satsuma \[TTS\] proved a conservation law for their box-ball system (soliton cellular automata) using the Robinson-Schensted-Knuth correspondence: we associate a permutation to each state $`p`$, which we call the stack permutation of the state, then the shape of the $`P`$-symbols of these stack permutations is conserved. We denote this partition by $`\lambda (p)`$.

Recently, it was observed that there exists a crystal structure behind this system, and the identification of this box-ball system with a box-ball system arising from $`A_1^{(1)}`$-crystal was made in \[HHIKTT\]. In this crystal picture, the time evolution is described by combinatorial row-to-row transfer matrices, and the energy functions $`E_l(p)`$ $`(l)`$ of this system naturally gives us another conservation law \[FOY\]. Further, it was conjectured how these laws were related. It is given by a simple formula:

$$E_l(p)E_{l1}(p)=\lambda _l(p)$$

where $`\lambda _l(p)`$ is the length of the $`l`$ th row of the partition $`\lambda (p)`$.

In terms of the lengths of solitons $`N_1,N_2,\mathrm{}`$, $`\lambda (p)`$ is the partition which has $`N_k`$ columns of length $`k`$ $`(k)`$, and $`E_l(p)=_k\mathrm{min}(l,k)N_k`$.

For example, if the state is an asymptotic soliton state,

$$\mathrm{}01^{k_1}0^{l_1}1^{k_2}0^{l_2}\mathrm{}(l_1,l_2,\mathrm{}>>0)$$

it is straightforward to verify it.

The purpose of this short note is to prove the formula. It is done by supplying conceptual explanation about the appearence of stack permutations and their $`P`$-symbols.

Main idea is to interprete the box-ball system as a discrete dynamical system of a path on $`\times `$, from which we naturally read off the evolution of the permutation poset of the stack permutation of the state of the original box-ball system. This gives us a natural explanation why the Torii-Takahashi-Satsuma law holds.

We then turn to the crystal picture, and describe the sites which contribute to the energy function by using stack permutations. This explains why these energy functions are related to stack permutations.

Our conclusion is that the depth of stacks explains both conservation laws, which proves the relation of these laws.

The author hopes that this explanation would be valid after modifications in the case of $`A_r^{(1)}`$ soliton cellular automata. In this case, Nagai’s conserved quantities remain mysterious from combinatorial point of view.

## 2 Fomin-Greene theory on posets

We start with the Fomin-Greene theory of posets. Good references are \[BF\] and \[F\]. Let $`(𝒫,)`$ be a poset. A chain is a totally ordered subset of $`𝒫`$. An antichain is a subset of $`𝒫`$ on which no two elements are comparable.

###### Definition 1

Let $`(𝒫,)`$ be a poset. We define $`I_k(𝒫),D_k(𝒫)`$ for $`k`$ as follows.

$$\begin{array}{c}I_k(𝒫):=\mathrm{max}\{|C_1\mathrm{}C_k||\text{ }C_i\text{:(possibly empty) chain}\}\\ D_k(𝒫):=\mathrm{max}\{|A_1\mathrm{}A_k||\text{ }A_i\text{:(possibly empty) antichain}\}\end{array}$$

We also define $`\lambda _k(𝒫),\lambda _k^{}(𝒫)`$ for $`k`$ by their differences:

$$\lambda _k(𝒫):=I_k(𝒫)I_{k1}(𝒫),\lambda _k^{}(𝒫):=D_k(𝒫)D_{k1}(𝒫)$$

We thus obtain two compositions

$$\begin{array}{c}\lambda (𝒫):=(\lambda _1(𝒫),\lambda _2(𝒫),\mathrm{})\\ \lambda ^{}(𝒫):=(\lambda _1^{}(𝒫),\lambda _2^{}(𝒫),\mathrm{})\end{array}$$

The following theorem is due to Greene and Fomin, which justifies the use of the notation $`\lambda ^{}(𝒫)`$.

###### Theorem 2

Let $`𝒫`$ be a poset. Then $`\lambda (𝒫)`$ and $`\lambda ^{}(𝒫)`$ are partitions. Further, $`\lambda ^{}(𝒫)`$ is the transpose of $`\lambda (𝒫)`$.

Let $`x=x_1\mathrm{}x_n`$ be a word in $`[1,r]^n`$, where $`[1,r]:=\{1,2,\mathrm{},r\}`$ is the set of alphabets. For a pair $`(T,k)`$ of a semistandard tableau and $`k[1,r]`$, we have the (row) insertion algorithm which produces another semistandard tableau. We denote this semistandard tableau by $`Tk`$.

###### Definition 3

Let $`x=x_1\mathrm{}x_n[1,r]^n`$ be a word. The semistandard tableau $`P(x)`$ is defined by

$$P(x)=\mathrm{}x_1x_2\mathrm{}x_n$$

and is called the $`P`$-symbol of $`x`$.

###### Definition 4

Let $`S_n`$ be the symmetric group of $`n`$ letters acting on $`[1,n]`$. For $`wS_n`$, its permutation poset $`(𝒫(w),)`$ is defined by

$$\begin{array}{c}𝒫(w)=\{(i,w(i))|i[1,n]\}\\ \\ (i,w(i))(j,w(j))ij,w(i)w(j)\end{array}$$

We identify $`wS_n`$ with the word $`w(1)w(2)\mathrm{}w(n)[1,n]^n`$.

Let $`x[1,r]^n`$ be a word, and assume that $`k`$ appears $`n_k`$ times in $`x`$. Then we see $`x`$ as a distinguished coset representative of $`S_n/S_{n_1}\times \mathrm{}\times S_{n_r}`$. Thus we can consider permutation posets for arbitrary $`x[1,r]^n`$, which we denote by $`𝒫(x)`$.

###### Example 5

Let $`x=312143[1,4]^6`$. Then to see it as a distinguished coset representative (an element of $`S_6`$) is the same as seeing it as $`3_11_121_243_2`$. Here, we use $`1_1<1_2<2<3_1<3_2<4`$ instead of $`1<2<3<4<5<6`$. The permutation poset $`𝒫(x)`$ is as follows.

We have $`I_1=3,I_2=5,I_3=6,I_4=6,\mathrm{},`$ and $`\lambda (𝒫(x))=(3,2,1)`$.

For permutation posets, the following is well known.

###### Theorem 6

Let $`x[1,r]^n`$ be a word, and $`𝒫(x)`$ be its permutation poset. Then $`\lambda (𝒫(x))`$ equals the shape of the $`P`$-symbol $`P(x)`$.

## 3 Box-Ball system

We now recall the box-ball system. Each state is given by an infinite sequence of $`\{0,1\}`$ which has finitely many $`1`$’s. We denote by $`1_1,\mathrm{},1_N`$ these $`1`$ read from left to right. The description of one step time evolution is very simple: for $`k=1,\mathrm{},N`$, we move $`1_k`$ to the leftmost $`0`$ among those which sit on the right hand side of $`1_k`$. We give an example.

###### Example 7

$$\begin{array}{cc}t:& \mathrm{}\mathrm{0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}0}\mathrm{}\hfill \\ & \\ t+1:& \mathrm{}\mathrm{0\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}0}\mathrm{}\hfill \end{array}$$

It is visualized as follows, and in fact this is the original description of the rule.

For a state, we shall define a finite sequence consisting of ”$`0`$”, ”$`(`$” and ”$`)`$”. We first choose subsequences of the form ”$`\mathrm{1\hspace{0.17em}0}`$”. These are called pairs of stack depth $`1`$. We change these pairs ”$`\mathrm{1\hspace{0.17em}0}`$” to ”$`()`$”. We then delete these pairs from the original state, and choose subsequences ”$`\mathrm{1\hspace{0.17em}0}`$” again. These pairs are called pairs of stack depth $`2`$. We change these pairs ”$`\mathrm{1\hspace{0.17em}0}`$” to ”$`()`$”, and delete these pairs again. We continue this procedure repeatedly until all $`1`$ are deleted. In the end, all $`1`$ are made into pairs with ”)” and their stack depths are defined. We say that an opening parenthesis and a closing parenthesis are matched if they belong to a same pair. We now define the stack permutation of the state.

###### Definition 8

For a state, we associate a finite sequence of ”$`0`$”, ”(” and ”)” as above. We read the opening parentheses from left to right, and number them $`1,2,`$ accordingly.

We number a closing parenthesis $`k`$ if it makes a pair with the $`k`$ th opening parenthesis. The permutation obtained by reading the numbering of closing parentheses from left to right is called the stack permutation of the state.

###### Example 9

For a state given by

$$\mathrm{}\mathrm{0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}0\hspace{0.17em}0}\mathrm{},$$

we make pairs as follows.

$$\begin{array}{c}\mathrm{}\mathrm{0\hspace{0.17em}0}()\mathrm{\hspace{0.33em}0\hspace{0.17em}1}()\mathrm{\hspace{0.33em}1}()\mathrm{\hspace{0.33em}0\hspace{0.17em}0\hspace{0.17em}0}\mathrm{}\\ \mathrm{}\mathrm{0\hspace{0.17em}0}()\mathrm{\hspace{0.33em}0\hspace{0.17em}1}()(())\mathrm{\hspace{0.33em}0\hspace{0.17em}0}\mathrm{}\\ \mathrm{}\mathrm{0\hspace{0.17em}0}()\mathrm{\hspace{0.33em}0}(()(()))\mathrm{\hspace{0.33em}0}\mathrm{}\end{array}$$

Thus the numbering of the closing parentheses is given by

$$\begin{array}{c}\mathrm{}\mathrm{0\hspace{0.17em}0}()\mathrm{\hspace{0.33em}0}(()(()))\mathrm{\hspace{0.17em}0}\mathrm{}\\ \mathrm{}\mathrm{\hspace{0.33em}1}\mathrm{\hspace{0.33em}3}\mathrm{\hspace{0.33em}5\hspace{0.17em}4\hspace{0.17em}2}\mathrm{}\end{array}$$

That is, the stack permutation is $`\mathrm{1\hspace{0.17em}3\hspace{0.17em}5\hspace{0.17em}4\hspace{0.17em}2}`$.

Note that an opening parenthesis does not move to the corresponding closing parenthesis, but the total set of the opening parentheses moves to the total set of the closing parentheses as a whole in each step of time evolution.

We now describe the box-ball system as a discrete dynamical system of a path on $`\times `$. Recall that we have defined a finite sequence of ”$`0`$”, ”(” and ”)” for each state. We read the sequence from the first ”(”, and associate a walk on $`\times `$ starting from $`(0,0)`$ by the rule that ”$`(`$” corresponds to ”$``$”, and ”$`0`$”, ”$`)`$” correspond to ”$``$”. We see it as a walk from $`(\mathrm{},0)`$ to $`(\mathrm{},N)`$, where $`N`$ is the number of balls in the box-ball system, by adding infinitely many ”$``$” to both sides. We give an example.

###### Example 10

For a state in the previous example, we have obtained the following sequence.

$$()\mathrm{\hspace{0.33em}0}(()(()))$$

Thus we have a walk as follows. The permutation poset of the stack permutation $`w=13542`$ is given by circles.

Let us think of a vertical edge corresponding to ”(” of a subsequence ”0 ( ”. Since the left neighbor of the vertical edge is ”$`0`$”, we have that a positive number of ”(” remain waiting to be matched as long as the number of ”(” is greater than the number of ”)”. It implies that if the vertical edge starts at $`(x_0,y_0)`$, and the walk is in the area $`yy_0>xx_0`$, then all edges are ”(” or ”)”, and no ”$`0`$” appears. Using this, we can give a simple evolution rule of the path: Let us start walking with $`(x_0,y_0)`$ and continue walking on the path until it hits the line $`yy_0=xx_0`$ again. We say that these edges constitute a group. We partiton the edges of the walk into such groups. For each group with $`(x_0,y_0)`$ as above, we reflect the edges in this group with respect to the line $`yy_0=xx_0`$. We then move the zero of the $`(x,y)`$-plane to the lower vertex of the first vertical edge. These complete one step time evolution. This evolution rule is best understood by an example.

Since the closing parentheses of the pairs of ”(” and ”)” define the stack permutation, we have natural correspondence between the vertices of the permutation poset of the stack permutation and the matching pairs: for each vertical edge, we choose the horizontal edge in the matching pair. If the vertical edge is on the $`j`$ th row and the horizontal edge is on the $`i`$ th column, then $`(i,j)`$ is an element of the permutation poset. In the above example, the leftmost three vertices correspond to the pairs of stack depth $`1`$, and the left vertex of the remaining two corresponds to the pair of stack depth $`2`$, and the rightmost vertex corresponds to the pair of stack depth $`3`$.

To describe the poset structure, it is convenient to describe the matching pairs by framed boxes. For the left example of the above, we have

( )( ( )( ( ) ) )

We say that two pairs are in outer relation if each of the corresponding framed box does not contain the other. If one framed box contains the other, we say that these are in inner relation. Note that two pairs are either in outer relation or inner relation by the definition of the pairs.

###### Lemma 11

Assume that two pairs of ”$`(`$” and ”$`)`$” are in outer (resp. inner) relation. Then the corresponding vertices in the permutation poset are comparable (resp.not comparable).

(Proof) If they are in outer relation, their positions in the sequence of ”(” and ”)” are given by

$$\mathrm{}(\mathrm{})\mathrm{}(\mathrm{})\mathrm{}.$$

If the left ”(” is the $`i`$ th opening parenthesis and the right ”(” is the $`j`$ th opening parenthsis, we have $`i<j`$ and $`w(i)<w(j)`$ by the definition of the stack permutation. Hence the corresponding vertices in the permutation poset are comparable. The argument for the inner case is similar. $`\mathrm{}`$

###### Proposition 12

To each state, we associate the permutation poset of the stack permutation of the state as above. Then its vertices corresponding to pairs of stack depth $`k`$ form a chain in the poset. We call it the depth $`k`$ chain and denote it by $`C_k`$. We then have that $`|C_1\mathrm{}C_k|`$ gives the maximal number of vertices covered by $`k`$ chains.

(Proof) By Lemma 11, $`C_k`$ is a chain. We show that this permutation poset admits decomposition into disjoint union of antichains $`A_k`$ such that each $`A_k`$ has the form $`\{v_1,\mathrm{},v_{l_k}\}`$ where $`v_i`$ is a vertex corresponding to a pair of stack depth $`i`$. Assume that we have already distributed vertices of stack depth smaller than $`k`$ into such antichains. In the definition of the pairs, it corresponds to the stage that we have deleted ”$`\mathrm{1\hspace{0.17em}0}`$”’s $`k1`$ times. By the definition of the pairs, each framed box of stack depth $`k`$ contains a framed box of stack depth $`k1`$, and these framed boxes of stack depth $`k`$ are in outer relation. The latter implies that we can choose distinct framed boxes of depth $`k1`$ for framed boxes of depth $`k`$. Since the vertices of stack depth $`k1`$ are distributed to distinct antichains, we can distribute the vertices of stack depth $`k`$ to antichains without violating the required property.

We now assume that $`C_1^{}\mathrm{}C_k^{}`$ gives the maximal number of vertices covered by $`k`$ chains. Since each antichain intersects $`C_1^{}\mathrm{}C_k^{}`$ at most $`k`$ times, we can move these vertices into $`C_1\mathrm{}C_k`$ keeping them mutually distict. This is possible by the existence of the antichain decomposition we have just proved. Hence, $`|C_1^{}\mathrm{}C_k^{}|`$ can not exceed $`|C_1\mathrm{}C_k|`$. $`\mathrm{}`$

If we denote by $`𝒫^t`$ the permutation poset at time $`t`$, and by $`C_k^t`$ the depth $`k`$ chain of $`𝒫^t`$, we have $`\lambda _k(𝒫^t)=|C_k^t|`$ by Proposition 12. By Theorem 6, we have that $`\lambda (𝒫^t)`$ is nothing but the shape of the $`P`$-symbol of the stack permutation of the state at time $`t`$. Hence, the following theorem is almost obvious. It simply says that the length of depth $`k`$ chain is conserved, which is easily seen from the evolution rule of the path as follows.

###### Theorem 13 (\[TTS\])

For each state, we compute its stack permutation. Then the shape of its $`P`$-symbol is conserved under time evolution.

(Proof) We show that $`|C_k^t|`$ is conserved. For $`k=1`$, the elements of the chain correspond to convex corners of the path. Hence it is obviously conserved by the evolution rule of the path. We then delete the depth $`1`$ chain from the posets. This is the same as deleting convex corners from the original path and the reflected path.

To know that the deletion of convex corners from the original path and the reflected path gives a same walk, it is enough to see that deleting convex corners gives the same walk as deleting concave corners. To compare the location of $`1`$’s in the walks, We divide the cases by looking at vertical lines (the middle lines of the figures below). For the location of $`0`$’s, we divide the cases by looking at horizontal lines and the argument is entirely similar, which we omit.

The leftmost figure represents the case that we have vertical lines on both sides. One may subdivide the case into four by separating the case that there is exactly one $`0`$ in the middle of $`1`$’s from the case that there are more than one $`0`$’s in the middle of $`1`$’s, if one wishs. The remaining two cases are the left end and the right end of the walk.

By comparing the results of the deletion of the concave corners and the convex corners, we know that the new walks are the same. In particular, deleting convex corners from the original path and the reflected path give a same walk. Since the depth $`2`$ chain becomes the depth $`1`$ chain of the new poset, we can apply the same argument to conclude that $`|C_2^t|`$ is conserved. By repeated use of the argument, we also have the conservation laws for all $`k`$. $`\mathrm{}`$

## 4 Energy functions and stack permutations

We now turn to the crystal description of the box-ball system. Let $`B:=\{\text{0},\text{1}\}`$ be the $`A_1`$ crystal associated with the vector representation whose highest weight vector is 0. Its affinization is denoted by $`Aff(B):=\times B`$. This is an $`A_1^{(1)}`$ crystal. Note that the numbering of 0 and 1 is different from the usual one. $`B`$ is identified with the subset $`\{0\}\times B`$. For each state, we cut sufficiently remote $`0`$’s and consider it as an element in $`B^n`$.

To describe the time evolution rule, we take the crystal of the $`l`$ th symmetric tensor $`B_l`$ and its affinization $`Aff(B_l)`$ with $`l`$ sufficiently large. The elements of $`B_l`$ are nondecreasing sequences of length $`l`$ whose entries are $`0`$ and $`1`$. We write $`0^{m_1}1^{m_2}`$ $`(m_1+m_2=l)`$ for these elements. We use combinatorial $`R`$ matrices to get isomorphism of affine crystals as follows.

After we embed $`B^n`$ to $`Aff(B)^n`$, we apply this combinatorial row-to-row transfer to the tensor product of $`0^l`$ with the upper $`Aff(B)^n`$ to get the lower $`Aff(B)^n`$ tensored by $`0^l`$. Then we forget the symmetric tensor part and the $``$ part of the affine crystal. The result is an element of $`B^n`$. This procedure gives one step time evolution of the box-ball system.

We consider the isomorphism for arbitrary $`l`$. Then for a state $`p`$, we have

$$\left(0\times \overline{)0^l}\right)pp^{}\left(E_l(p)\times \overline{)0^{m_1}1^{m_2}}\right)$$

for some $`p^{}Aff(B)^n`$ and $`m_1,m_2`$. These $`E_l(p)`$ are called energy functions. It is known \[NY\] that if we set $`E_l=0`$ and increase it by one at the sites of the following form, then the final value of $`E_l`$ coincides with $`E_l(p)`$.

By using the fact that time evolution is obtained from a crystal isomorphism of affine crystals, Fukuda, Okado and Yamada \[FOY, Theorem 3.2\] have proved that these $`E_l(p)`$ are conserved quantities of this box-ball system.

The purpose of this section is to relate these quantities to the stack permutation of the state $`p`$.

###### Theorem 14

For each state $`p`$, we define the sequence of ”$`0`$”, ”(” and ”)” as in the previous section. Then $`E_l`$ increases precisely at the sites corresponding to ”)”.

(Proof) Assume that ”)” corresponds to a pair of depth $`k`$. We shall show that if $`lk`$, $`E_l(p)`$ does increase at this site.

Let ”(” be the corresponding opening parenthesis, and $`1^{k_1}0^{l_1}\mathrm{}1^{k_N}0^{l_N}`$ be the walk starting from the vertical edge corresponding to the ”(” and ending at the horizontal edge corresponding to the ”)”. We write the evolution of the symmetric tensor along the path as follows.

Assume that $`E_l`$ does not increase at this site. Then the status of the symmetric tensor on the both ends of the last edge is $`0^l`$. We denote by $`e_1,\mathrm{},e_N`$ the last edges of $`1^{k_1},\mathrm{},1^{k_N}`$ respectively. Then we can prove the following by downward induction on $`i`$.

* The symmetric tensor on the upper end of $`e_i`$ has the form $`0^{ls}1^s`$ with $`s<l`$.
* Since the upper end of $`e_i`$ is not saturated, we have steady increase of the number of $`1`$ in the symmetric tensor during $`k_i`$ vertical edges, and no saturation occurs during these edges.
* Since $`l_{i1}k_i+\mathrm{}k_N+l_Nk`$, the symmetric tensor on the upper end of $`e_{i1}`$ has the form $`0^{ls}1^s`$ with $`s<l`$.

We then have the following.

* Since $`k_1l_1+\mathrm{}l_i0`$, we have steady decrease of the number of $`1`$ in the symmetric tensor during $`l_i`$ horizontal edges, and no saturation occurs.

Therefore, we conclude that the left end of the last edge of the walk has the symmetric tensor of the form $`0^{ls}1^s`$ with $`s>k_il_i0`$, which contradicts the assumption at the beginning. (In particular, we have that $`I_l(𝒫)E_l(p)`$ where $`𝒫`$ is the permutation poset of the stack permutation of the state $`p`$.)

We nextly show that if $`l<k`$, then $`E_l`$ does not increase at this site. To prove this, we show for arbitrary $`k,l`$ that the right end of the last edge of the walk has the symmetric tensor $`0^l`$ if the stack depth $`k`$ is equal or greater than $`l`$, and $`0^{m_1}1^{m_2}`$ $`(m_1k)`$ if $`k`$ is equal or smaller than $`l`$. We prove it by induction on $`k`$. If $`k=1`$, the proof is obvious. If $`kl`$, we choose the last closing parethesis of stack depth $`k1`$. Then by the induction hypothesis, the symmetric tensor has the form $`0^{m_1}1^{m_2}`$ $`(m_1k1)`$ at this site. Note that we have already proved that no saturation occurs during the walk if $`kl`$. Hence, if we start the walk with $`0^{m_1}1^{m_2}`$, we end the walk with $`0^{m_1}1^{m_2}`$. From this, we know that the symmetric tensor at the left end of the last edge of the walk is also $`0^{m_1}1^{m_2}`$ $`(m_1k1)`$. Hence, the right end of the last edge has the form $`0^{m_1^{}}1^{m_2^{}}`$ $`(m_1^{}k)`$ if $`m_2>0`$ and $`0^l`$ if $`m_2=0`$. But we also have $`lk`$ in the latter case.

We now assume that $`kl`$. We choose the last closing parenthesis of stack depth equal or greater than $`l`$. Since its stack depth is smaller than $`k`$, we can apply the induction hypothesis to know that the symmetric tensor has the form $`0^l`$ at this site. Further, since we have pairs of stack depth less than $`l`$ during this site and the last edge of the walk, we have that the symmetric tensor at the left end of the last edge has $`0^l`$. Thus the same is true for the right end of the last edge.

Therefore, we have proved that $`E_l`$ increases precisely at the sites corresponding to ”)”. (In particular, we have also proved that $`I_l(𝒫)=E_l(p)`$.) $`\mathrm{}`$

## 5 Conclusion

For a state $`p^t`$ at time $`t`$, we denote by $`𝒫^t`$ the permutation poset of the stack permutation of the state $`p^t`$. Then the energy function $`E_l(p^t)`$ counts the vertices of $`𝒫^t`$ whose stack depth are equal or less than $`l`$. On the other hand, the $`l`$ th row of the shape $`\lambda (p^t)`$ of the $`P`$-symbol of the stack permutation is equal to the number of vertices of $`𝒫^t`$ whose stack depth are $`l`$. Hence these quantities are naturally explained by the notion of stack depth, and we have $`I_l(𝒫^t)=E_l(p^t)`$.

Further, the evolution rule of a path naturally explains why these quantities are conserved.

Tokyo University of Mercantile Marine,

Tokyo 135-8533, Japan.


format: text

text:
# The CLEO III Ring Imaging Cherenkov Detector

## 1 Introduction

A crucial element of all the experiments exploring heavy flavour phenomenology is the particle identification system, with the primary goal of identifying $`\pi `$’s, $`K`$’s and $`p`$’s. CLEO III is instrumented with a Ring Imaging Cherenkov detector (RICH) implemented in a barrel geometry comprising a total of 20 cm radial space between the tracking system and the CsI electromagnetic calorimeter.

Our goal is to achieve a $`\pi /K`$ separation greater than $`3\sigma `$ up to the maximum momentum of particles produced in $`B`$ meson decays at the $`\mathrm{{\rm Y}}(4S)`$, $`2.65`$ GeV/c, the momentum characteristic of the very important decays $`B\pi \pi /K\pi `$. The $`\pi /K`$ separation at this momentum is 14.3 mr. Thus the required angular resolution per track is $`\sigma _t4.8`$ mr. Our design parameters are a resolution per photon, $`\sigma _\gamma `$, of 14 mr and approximately 12 detected photons. This corresponds to $`\sigma _t=4.2`$ mr. Note that we expect $`2\sigma \pi /K`$ separation from dE/dx information from the drift chamber enclosed in our system. Our Monte Carlo simulation, confirmed by an extensive set of test beam data , predicts that the CLEO RICH detector meets our goals. These expectations will be compared with preliminary data from an engineering run that took place during the first few months of detector operation.

## 2 Detector description

Figure 1 shows a $`r\varphi `$ section of the CLEO III RICH detector system. It involves an inner cylindrical shell comprising thin LiF crystal radiators , where charged particles moving with a speed greater than the speed of light in this medium produce the Cherenkov cone. Note that when the charged particle crosses the detector at normal incidence, the Cherenkov light is trapped inside the radiator because of total internal reflection. We have demonstrated  that this problem can be circumvented by shaping the outer surface of the radiator like the teeth of a saw, hence the name “sawtooth radiator”. The production of these radiators is more difficult because of the necessity of unconventional polishing techniques. Thus we are employing them only up to a dip angle of 22.

The outer cylindrical shell of the detector is composed by thin gap multiwire chambers with their cathode plane finely segmented into 7.5$`\times `$8 mm2 pads, to reconstruct the location of the Cherenkov photons at the detector surface with high precision. The photosensitive element in these chambers is triethylamine (TEA), dispersed in $`\mathrm{CH}_4`$. The sensitive bandwidth of $`\mathrm{CH}_4`$-TEA is centered around 150 nm. Therefore we need to detect vacuum ultraviolet photons. Thus the interior faces of these chambers are made of $`\mathrm{CaF}_2`$ windows, coated with metal strips that provide one of the two cathode planes. The use of $`\mathrm{CaF}_2`$ windows and the need to seal tightly this system posed a great challenge to the mechanical design.

The cathode pad signals are processed by low noise front end electronics devices located at the back of the cathode planes. The fine segmentation needed involves a relatively complex readout system, featuring digitization and sparsification of 230,400 charge signals. A low noise ASIC called VA\_RICH, custom made for this application by IDE AS, Norway, converts the charge induced in the cathode pad into a differential current signal. It features excellent noise performance, measured to be :

$$ENC(e^{})=130+9\times C_{in}(\mathrm{pF})$$
(1)
This corresponds to an equivalent noise charge expected to be of the order of 250 $`e^{}`$ for the input capacitance seen by the preamplifier, dominated by the relatively long traces connecting the chamber pad to the chip input. The goal of keeping the electronic noise to a minimum is dictated by the single photon response, exponential in shape for the moderate gains at which we are planning to operate our chambers. This chip features also a high dynamic range, of the order of a few 100,000 electrons, in order to measure also the pulse height produced by the minimum ionizing particle, 20-50 times higher than the single photon pulse. The differential current signal is digitized and sparsifed in the data boards residing in VME crates located a few meters away from the detector. The challenge of maintaining the low noise performance in such a complex system is mitigated by the implementation of an online common mode subtraction prior to the sparsification process.

## 3 Preliminary performance characterization

The RICH detector was installed in CLEO in early August 1999 and the first data were collected starting in November 1999. The detector performance is very satisfactory. The transparency of the expansion volume is above 99% at 150 nm. The high voltage system is very stable and has been successfully operated at several different voltage settings. The average pad gain in the data reported here is $`2.5\times 10^4e^{}`$. The electronics noise, prior to any working point optimization, has a mean value of 500 $`e^{}`$, upon common mode subtraction, and of 750 $`e^{}`$ including the coherent component. This performance is quite remarkable for such a large system, powered by boards with a mixed digital and analog environment.

During the engineering run, the detector has been characterized in a “stand alone” mode, without resorting to any kind of external tracking information, not yet available. The data set used for this preliminary study is composed primarily of Bhabha events. For these events, having the simple topology of back to back charged tracks of opposite charge, we can reconstruct the track helix parameters using the centroid of the clusters corresponding to the $`e^+`$ and $`e^{}`$ intersection with the RICH pad chambers, combined with constraints determined by the event kinematics and the beam properties. We selected the large pulse heights to be charged track clusters. In addition we set their momentum equal to the beam energy and we assumed that the two tracks originated from a point with coordinates $`x_0=y_0=0`$, exploiting the small transverse beam size, and $`z_0=(z_{e^+}+z_e^{})/2`$, where $`z`$ is the coordinate along the beam direction.

Fig. 2 shows the Cherenkov angle distribution for single photons. These data are fit with a signal shape plus a polynomial background function. We find $`\sigma _\gamma =14.7\pm 0.2`$ mrad for flat radiators and $`\sigma _\gamma =12.2\pm 0.4`$ mrad for sawtooth radiators. These results are preliminary because the real tracking information from the combined silicon strip detector and drift chamber was not available. Moreover the RICH alignment is still underway. The number of detected photoelectrons is consistent with our expectations of about 12 $`\gamma `$’s per track. Thus we can be confident that, once the tracking information is available and the alignment procedure is complete, we will achieve excellent particle separation at all the momenta of interest.

## 4 Conclusions

The first few months of operation of the CLEO III RICH have been quite successful. The detector has performed in a reliable manner and has proven to meet our specification. It will be a key element in the CLEO physics program in the next few years.

## 5 Acknowledgements

We would also like to thank T. Ypsilantis for the inspiring role that he had at the beginning of this work and for the great creative energy that he continues to bring to the field of particle identification.


text:
# Periodically kicked turbulence

## Abstract

Periodically kicked turbulence is theoretically analyzed within a mean field theory. For large enough kicking strength $`A`$ and kicking frequency $`f`$ the Reynolds number grows exponentially and then runs into some saturation. The saturation level $`Re^{sat}`$ can be calculated analytically; different regimes can be observed. For large enough $`Re`$ we find $`Re^{sat}Af`$, but intermittency can modify this scaling law. We suggest an experimental realization of periodically kicked turbulence to study the different regimes we theoretically predict and thus to better understand the effect of forcing on fully developed turbulence.

Periodically driven flow is ubiquitous. Faraday’s experiment  is an idealized version thereof, more relevant examples are the earth’s atmosphere, driven by periodical heating of the sun, or the blood flow in veins, driven by the beating heart. Another example is the gas flow inside a sonoluminescing bubble which is periodically kicked by the collapsing bubble wall . Another example of periodically kicked turbulence is the numerical realization of homogeneous shear flow  where periodical remeshing is necessary.

In this paper we set up a mean field theory for periodically kicked flow, based on the mean field theory for decaying turbulence  which was able to describe the experimentally measured energy decay in turbulent liquid helium flow with fixed external length scale . The goal of this paper is to theoretically understand the different flow regimes which are to be expected, to explore the effect of intermittency corrections on these regimes, and to ultimately initiate experiments.

Another motivation for the paper is to study the effect of a specified type of forcing on turbulence. In most theoretical studies on turbulence a Gaussian random noise, acting on the largest length scales, is assumed. Only recently experimentalists started to systematically vary the type of forcing . This work is a further step towards the analysis of a more specific type of forcing.

To define the model, we have to (i) calculate the energy input during the kick and (ii) know how the energy is dissipated in the time $`\mathrm{\Delta }t`$ between successive kicks.

(i) Kick: As an illustration, consider plane shear flow in the 1-direction; the flow is sheared in the 3-direction. The width of the channel is $`L`$, the velocity of the upper plate is $`U`$, the lower one is at rest. The average energy dissipation rate can be calculated to be (see e.g. )

$$ϵ=\frac{U}{L}\left(u_3u_1_A\nu _3u_1_A\right),$$
(1)
where $`_A`$ denotes the average over the x-y-plane and $`\nu `$ is the viscosity. For laminar flow the first term in the bracket is zero and the second one is $`\nu U/L`$. In turbulent flow in the middle of the channel the second term on the rhs will hardly contribute for large enough Reynolds numbers. The first term represents the total turbulent flow energy $`E`$, order of magnitude wise. Therefore, in general,

$$\dot{E}(t)=ϵ(t)\frac{U}{L}E(t)+\frac{U^2}{L^2}\nu .$$
(2)
Imagine now a short intense kick of time $`\mathrm{\Delta }t_{kick}\mathrm{\Delta }t`$ on the flow by rapidly moving the upper plate with $`U`$. After this kick the initial energy $`E_0`$ increases according to (2)

$$E_1=E_0+\left(E_0+\frac{U\nu }{L}\right)\frac{U}{L}\mathrm{\Delta }t_{kick},$$
(3)
where we have assumed $`\mathrm{\Delta }t_{kick}L/U`$. Assume isotropic turbulence in the flow center and define a Reynolds number ***Note that this is not the standard definition and gives lower values for the laminar-turbulent transition than what one is used to.

$$Re(t)=\frac{Lu_{1,rms}(t)}{\nu }=\sqrt{\frac{2}{3}}\frac{L\sqrt{E(t)}}{\nu }.$$
(4)
Then eq. (3) translates to

$$Re_1=Re_0\sqrt{1+2A+\frac{Re_{lam}^2}{Re_0^2}}$$
(5)
with the dimensionless kicking strength $`A=\frac{1}{2}\mathrm{\Delta }t_{kick}U/L1`$ and the “laminar” Reynolds number $`Re_{lam}=\frac{2}{3}\mathrm{\Delta }t_{kick}U^2/\nu `$. We choose this name because for very small $`Re_0`$ we have $`Re_1=Re_{lam}`$. For very large $`Re_0Re_{lam}`$ we have $`Re_1=(1+A)Re_0`$. Shear flow and the “derivation” of eq. (5) are only thought of as a motivation; there will be many other experimental situations where the energy input roughly corresponds to a law of type (5).

(ii) Decay: In ref.  we calculated how the turbulent activity decays within a time $`t`$ for flow with fixed external length scale $`L`$. The calculation was based on Effinger and Grossmann’s variable range mean field theory of turbulence  in which viscous subrange and inertial subrange can be treated equally well. The result of ref.  is that for given initial Reynolds number $`Re_i`$ the time dependence of $`Re(t)`$ (defined as in eq. (4)) is determined by the inverse function of

$$\frac{t(Re)}{\tau }=\frac{3}{c_{ϵ,\mathrm{}}}\left[F(Re)F(Re_i)\right]$$
(6)
where $`\tau =L^2/\nu `$ and $`F(Re)`$ is given by

$`F(Re)`$ $`=`$ $`{\displaystyle \frac{1}{2Re^2}}\left\{\gamma +\sqrt{\gamma ^2+Re^2}\right\}`$ (7)
$`+`$ $`{\displaystyle \frac{1}{2\gamma }}\mathrm{log}\left\{{\displaystyle \frac{\gamma +\sqrt{\gamma ^2+Re^2}}{Re}}\right\}.`$ (8)
$`c_{ϵ,\mathrm{}}=(6/b)^{3/2}`$ (in this theory) is the dimensionless energy dissipation in the large Reynolds number limit, $`\gamma =9/c_{ϵ,\mathrm{}}`$, and $`b`$ is the Kolmogorov constant  which is the only free parameter in the theory of ref. . From experiment  $`b=6.0`$, a value which we take in all calculations here. Consequently, $`c_{ϵ,\mathrm{}}=1.0`$ and $`\gamma =9.0`$. Rather than the Reynolds number eq. (4) one could also give the Taylor-Reynolds number $`Re_\lambda `$ ,

$$Re_\lambda =\sqrt{\frac{15Re^2}{c_{ϵ,\mathrm{}}\left(\gamma +\sqrt{\gamma ^2+Re^2}\right)}}.$$
(9)
Eqs. (5) - (8) with $`t=\mathrm{\Delta }t`$ (the time between successive kicks) fully define the present model. The two main physical parameters in the model are the kicking strength $`A`$ and the kicking frequency $`f=1/\mathrm{\Delta }t`$. The third physical parameter is $`Re_{lam}`$, the minimal Reynolds number after a kick. We pick $`Re_{lam}=1`$ throughout.

Figure 1 shows $`Re(t)`$ for fixed $`A`$ and three different kicking frequencies $`f`$ for the initial Reynolds number $`Re_0=Re_{lam}=1`$. During each cycle there is a kick $`Re_0Re_1`$ according to eq. (5) and a subsequent decay according to eq. (6). Overall, there is growth up to some saturation level $`Re^{sat}(A,f)`$, achieved after $`t_{sat}(A,f)`$. In this saturation state energy input and loss through decay in $`\mathrm{\Delta }t`$ balance and the degree of excitation fluctuates between a lower saturation level $`Re_l^{sat}`$ and an upper saturation level

$$Re_u^{sat}=Re_l^{sat}\sqrt{1+2A+(Re_{lam}/Re_l^{sat})^2}.$$
(10)
The (lower) level of saturation $`Re_l^{sat}(A,f)`$ is given through the implicit equation

$$\frac{1}{\tau f}=\frac{3}{c_{ϵ,\mathrm{}}}\left[F\left(Re_l^{sat}\right)F\left(Re_l^{sat}\sqrt{1+2A+\left(\frac{Re_{lam}}{Re_l^{sat}}\right)^2}\right)\right].$$
(11)
For large $`Re_u^{sat}`$, $`Re_l^{sat}\gamma `$ one has the explicit result

$$Re_u^{sat}=\frac{3\tau f}{c_{ϵ,\mathrm{}}}(\sqrt{1+2A}1)\frac{3\tau }{c_{ϵ,\mathrm{}}}Af.$$
(12)
In figure 2 we show a log-log plot of $`Re_{l,u}^{sat}`$ as a function of $`f`$ for different kicking strength $`A`$. Two regimes are seen: (i) For small $`Re_{l,u}^{sat}\text{ }\stackrel{<}{}\text{ }\gamma `$ a laminar regime in which after the decay the kick always brings back the level of excitation to the laminar value $`Re_{lam}=1`$. (ii) For large $`Re_{l,u}^{sat}\text{ }\stackrel{>}{}\text{ }\gamma `$ we have a turbulent scaling regime with $`Re_u^{sat}f`$. The transition from the laminar to the turbulence regime takes place around

$$f_{trans}=\frac{3}{\tau A}$$
(13)
Similarily, these two regimes are also seen in fig. 3 where we plotted $`Re_{l,u}^{sat}`$ vs $`A`$. In the turbulent regime for large $`Re_{l,u}^{sat}`$ it is $`Re_{l,u}^{sat}A`$.

We now come to the important question of how intermittency effects  change the exponents calculated within this mean field theory. In refs.  intermittency effects have been included into the mean field theory of ref.  on a phenomenological basis. One possibility for their effect is that the dimensionless energy dissipation rate $`c_ϵ`$ becomes slightly Reynolds number dependent even in the large Reynolds number limit , $`c_ϵRe^\kappa `$, with $`\kappa =(9/8)\delta \zeta _2/(1+3\delta \zeta _2/8)`$. Here, $`\delta \zeta _20.03`$ is the experimentally found deviation from the classical scaling exponent $`\zeta _2=2/3`$ of the second order velocity structure function. The consequences of this small ($`\kappa 0.03`$) scaling correction can straightforwardly be embodied in the present mean field approach to periodically kicked turbulence. The result is that in the turbulent regime the saturation level now obeys

$$Re_{l,u}^{sat}(Af)^{1/(1\kappa )}$$
(14)
rather than eq. (12). Equation (14) may offer a new and independent way to experimentally determine intermittency exponents.

Another effect related to intermittency is the following: We expect the total energy to build up again and again over time scales larger than $`\mathrm{\Delta }t`$ and then to suddenly drop because a energy pulse is traveling downscale. Such behavior has been observed in numerical simulations of periodically remeshed homogeneous shear flow  and in simulations of periodically kicked shell modells of turbulence \[results to be publihsed\]. As based on a mean field theory the model of this paper is only applicable to the mean energy and not to these fluctuations.

We now come back to an experimental realization. It will be easier to perform experiments in a closed system rather than in a channel flow. A particular suited experimental setup for periodically kicked turbulence would be the flow in a cylinder between two counterrotating disks . Also, Rayleigh-Benard convection may be well suited. Here, as an example we take Taylor-Couette flow : If the radii of the inner and outer cylinder are similar, the energy input will still roughly follow eq. (5). Take water as fluid which has $`\nu =10^6m^2/s`$ and take $`L=1cm`$. Then $`\tau =100s`$. Realistically achievable kick-strengths would be $`A=\mathrm{\Delta }t_{kick}U/(2L)0.1s0.1m/s/(21cm)=0.5`$, e.g., the whole range $`A<1`$ in which the theory is applicable. For this value the scaling regime sets in at $`f_{trans}=0.06Hz`$. Roughly two decades of scaling are necessary to explore intermittency corrections and to distinguish between eqs. (12) and (14). I.e., one has to go up to frequencies of around $`6Hz`$ which should be achievable. At these relatively large kicking frequencies measurements can only reveal instantaneous values. To get statements on the averaged quantities dealt with in this paper, ensemble averaging is necessary. This is best done by repeatedly probing the flow at some fixed phase after the respective kick. Averaging over the results at phase $`0^+`$ will give $`Re_u^{sat}`$, averaging over the results at phase $`\mathrm{\Delta }t^{}`$ will give $`Re_l^{sat}`$, etc.

We suggest to perform a periodically kicked turbulence experiment and to measure $`Re^{sat}`$ as a function of both $`A`$ and $`f`$. To our knowledge it would be one of the first ones with some active control on the type of forcing. Immediate questions to ask are: Does the level of saturation for large $`Re`$ indeed only depend on the product $`Af`$ as suggested by eqs. (12) and (14) or do boundary effects carry on into the strongly turbulent central regime and cause a more subtle relation? If so, an application of the so commonly used volume forcing for turbulence becomes more questionable. Do the (scaling) relations between the quantities introduced in this paper, e.g., $`Re_{l,u}^{sat}(f,A)`$ offer a new way to measure intermittency effects? What modifications arise if forcing and decay do not decouple as assumed in this simple model?

Acknowledgements: The author thanks L. Biferale, M. Brenner, B. Eckhardt, and F. Toschi for discussions. The work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). This research was also supported in part by the European Union under contract HPRN-CT-2000-00162 and by the National Science Foundation under Grant No. PHY94-07194 and we thank the Institute of Theoretical Physics in Santa Barabara for its hospitality.

e-mail address: lohse@tn.utwente.nl


text:
# Particle-Mesh Simulations of the Lyman-Alpha Forest

## 1 Introduction

Numerical simulations of structure formation in the universe incorporating hydrodynamics in Cold Dark Matter (CDM) dominated cosmologies have proven very successful in reproducing the statistical properties of the Ly$`\alpha `$forest as measured in high redshift Quasi-Stellar Object (QSO) spectra (Cen et al. 1994; Zhang, Anninos & Norman 1995; Hernquist et al. 1996; Zhang et al. 1997; Bond & Wadsley 1997; Theuns, Leonard & Efstathiou 1998a). The high level of statistical agreement suggests that the models are capturing the essential physical nature of the absorbing structures (Meiksin et al. 2000). Much of the structure of the Intergalactic Medium (IGM) may be understood as a consequence of the spatial coherence in the statistical properties of the density fluctuations of the dark matter alone in a “cosmic web”, as emphasized by Bond & Wadsley and supported by the slow comoving evolution of the integalactic filaments found by Zhang et al. (1998). Gnedin & Hui (1998) claim that in fact it is unnecessary to incorporate full hydrodynamics in the simulations to reproduce the properties of the Ly$`\alpha `$forest. Gnedin & Hui introduce instead a pseudo-hydrodynamical scheme into an N-body calculation, reporting comparable results in the distributions of flux and $`\mathrm{I}`$column density to full hydrodynamical $`+`$ N-body simulations for a standard CDM model. Pseudo-hydrodynamical methods based on PM alone have been used by Petitjean, Mücket & Kates (1995) and Croft et al. (1998).

Because of the large computational overhead in incorporating proper hydrodynamics into an N-body code, considerable benefit would be obtained if the hydrodynamics computations could be avoided. It would then be possible to perform a wide range of cosmological models at only moderate computational expense, as well as test systematics on the properties of the forest: the effect of box size, resolution, and even physical effects like variations in the character of the background UV radiation field that ionizes the gas. It would be of additional theoretical interest to infer the degree to which the properties of the forest are hydrodynamical in origin, as opposed to being solely a consequence of the clustering of dark matter.

To these ends, we perform simulations of the Ly$`\alpha `$forest for a variety of cosmological models using two approaches: a pure gravitational computation using a Particle-Mesh N-body code (PM simulations), and a pseudo-hydrodynamical technique, based on the method of Gnedin & Hui, incorporating an effective gas pressure force computed from the local dark matter density (HPM simulations).

## 2 The Models and Simulations

We consider two sets of models. The first is based on simulations in a flat CDM model with a cosmological constant ($`\mathrm{\Lambda }`$CDM ) using both the PM and HPM codes. These are performed in order to test convergence of the results for differing box sizes, particle numbers, and grid resolutions. A second set is based on 4 CDM models in Machacek et al. (2000), who performed computations using full hydrodynamics. We duplicate their models in order to compare our resulting statistics for the absorption features with theirs. The models investigated, in additional to the $`\mathrm{\Lambda }`$CDM model, are: a standard critical-density flat CDM model without a cosmological constant (SCDM), an open CDM model (OCDM), and the standard CDM model but with the power spectrum of the density perturbations tilted (tCDM) to match the normalization on large scales as determined from the COBE measurements of the Cosmic Microwave Background (Bunn & White 1997). The model parameters are shown in Table 1. As in Machacek et al. , the initial data were generated using the BBKS transfer function (Bardeen et al. 1986) to compute the starting redshifts and the unconstrained initial particle positions and velocity perturbations appropriate for each model. The non-dimensional linear power spectra $`\mathrm{\Delta }^2(k)`$ are shown at $`z=3`$ in Figure 1.

To account for the degree of ionization of the hydrogen, it is also necessary to include a background UV radiation field. We essentially circumvent this by fixing the mean optical depth, as defined by $`\overline{\tau }_\alpha =\mathrm{log}<\mathrm{flux}>`$, at $`\overline{\tau }_\alpha =0.30`$ at $`z=3`$, as in Machacek et al.

The evolution is computed using the PM code described in detail in Meiksin, White & Peacock  and White . The computational volume is chosen to be periodic with side length 1, and we take as our time coordinate the log of the scale factor, $`\mathrm{log}a`$, where $`a(1+z)^1`$. Velocities are measured in units of the expansion velocity across the box $`aHL_{\mathrm{box}}`$, where $`H`$ is the Hubble parameter fixed by the Friedmann equation

$$H^2H_0^2\left(\frac{\mathrm{\Omega }_{\mathrm{mat}}}{a^3}+\mathrm{\Omega }_\mathrm{\Lambda }+\frac{\mathrm{\Omega }_K}{a^2}\right)$$
(1)
in terms of the Hubble constant $`H_0`$ and the densities in matter, cosmological constant and curvature, in units of the critical density. The density is defined from the particle positions by assigning particles to a regular cartesian grid using the Cloud-In-Cell (CIC) scheme (Hockney & Eastwood ). The Fourier Transform is taken and the force computed using the kernel $`\stackrel{}{k}/k^2`$. We use a second order leap frog method to integrate the equations. The relevant positions are predicted at a half time step and used to calculate the accelerations which modify the velocities. The time step is dynamically chosen as a small fraction of the free-fall time with a maximum size of $`\mathrm{\Delta }\mathrm{log}a=3\%`$.

The initial conditions were created by displacing the particles from a uniform grid using the Zel’dovich approximation.

Our implementation of the HPM code differs slightly from that in Gnedin & Hui. Firstly we use 2 sets of particles, one of which feels only the gravitational force (the ‘dark matter’) and one of which additionally feels a pseudo-pressure force (the ‘gas’). The gravity is computed as before. To compute the ‘hydro’ forces we bin the ‘gas’ particles onto a grid using CIC as for the gravity. We compute the temperature using the equation of state

$$T=T_0\left(\frac{\rho }{\overline{\rho }}\right)^{\gamma 1}$$
(2)
where $`T_0`$ and $`\gamma `$ are obtained from Hui & Gnedin (1997). The enthalpy per particle

$$\phi =\frac{\gamma }{\gamma 1}T$$
(3)
is then computed, Fourier transformed and smoothed by a gaussian of half a grid cell (to supress high frequency noise). The ‘hydro’ force is then calculated from the inverse Fourier transform of $`i\stackrel{}{k}\stackrel{~}{\phi }`$ and added to the gravitational force for the ‘gas’ particles.

In practice, $`T_0`$ will depend on the reionization history of the IGM. In particular, if He $`\mathrm{II}`$is reionized (to He $`\mathrm{III}`$) late, lower values of $`T_0`$ will result prior to He $`\mathrm{II}`$reionization and higher values in rarified gas after (Meiksin 1994). The effect of reionization on the thermal properties of the gas will only be known once the nature of the ionizing sources has been established, and may eventually require the incorporation of radiative transfer into the simulations.

We consider a variety of particle number and grid zone combinations, as well as several box sizes. The PM simulations are summarised in Table 2. For the HPM simulations, two combinations are used: $`(N_p,N_g)=(128,256)`$ and $`(256,512)`$. For HPM the dark matter and the gas are represented by an identical number of particles, each equal to $`N_p^3`$.

Given a set of final particle positions and velocities, we compute the spectra as follows. First the density and density-weighted line-of-sight velocity are computed on a grid (using CIC interpolation as above) and Gaussian smoothed using FFT techniques. This forms the fundamental data set. A grid of sightlines is drawn through the box, parallel to the box sides and along each sightline a new 1D grid of density and velocity are obtained from the fundamental grid using CIC interpolation. The position in velocity units and the temperature are computed on this 1D grid using the cosmological model and equation (2) for the equation of state.

Using this 1D grid we then integrate in real space to find $`\tau (u)`$ at a given velocity $`u`$. Specifically we define

$$\tau (u)=A𝑑x\left[\frac{\rho (x)}{\overline{\rho }}\right]^2T(x)^{0.7}b^1e^{(uu_0)^2/b^2}$$
(4)
where $`u_0=xaHL_{\mathrm{box}}+v_{\mathrm{los}}`$ and $`b=\sqrt{2k_BT/m_\mathrm{H}}`$ is the Doppler parameter, where $`m_\mathrm{H}`$ is the mass of a hydrogen atom. The flux at velocity $`u`$ is $`\mathrm{exp}[\tau ]`$. The integration variable $`x`$ indicates the distance along the box in terms of the expansion velocity across the box. The constant $`A`$ is then iteratively adjusted to obtain a predefined $`\overline{\tau }_\alpha =\mathrm{log}\mathrm{exp}(\tau )`$.

We generate sample spectra from the simulated data along random lines–of–sight through the computational boxes over the length of each box. In order to make a fair comparison with the data in M00, both sets of spectra are rebinned to a constant pixel resolution of $`\lambda /\mathrm{\Delta }\lambda =74000`$ and gaussian smoothed to mimic a spectral resolution of 8 km s-1 . These parameters were chosen to approximate what may be achieved using the Keck HIRES.

There are several measurements of the Ly$`\alpha `$forest which may be used as a basis for comparison. The most fundamental is simply the distribution of flux per pixel in the synthesized spectra. The fluctuations in the flux on different velocity scales provide additional constraints, which may be related both to the density fluctuation power spectrum and to the thermal widths of the lines. We quantify the fluctuations in three ways. We compute the power spectrum $`P_F(k)`$ of the flux, an analog of the density power spectrum. (For the $`P_F(k)`$ estimates we use the raw spectra from the PM and HPM runs.) The widths of the spectral features may be directly characterised using wavelets, which quantify the changes in a given spectrum on being smoothed from one velocity resolution to another (Meiksin 2000). A multiscale analysis based on the Daubechies wavelets effectively performs a weighted smoothing of the spectrum over successive doublings of the pixel width. We use a Daubechies wavelet of order 20. Lastly, the most traditional method to quantify the flux fluctuations is by decomposing each spectrum into a set of Voigt line profiles. We utilise AutoVP (Davé et al. 1997) to perform such an analysis. Each sample spectrum produces an average of on the order of 900 lines per unit redshift. Typically $`(36)\times 10^4`$ lines are analysed per model. The absorption lines are characterised by their $`\mathrm{I}`$column densities and Doppler parameters. The $`\mathrm{I}`$column densities reflect both the density and physical thickness of the structures which give rise to the absorption features. The Doppler parameters reflect the degree of line broadening due both to the thermal motion of atoms and any internal peculiar velocities within the structures, and may provide probes of the temperature and velocity history of the universe at high redshifts.

To ensure our results are not affected by sample variance, we performed two independent simulations of some of the models. We found convergence on the statistical properties could be achieved with a single realization provided the statistics were averaged over a sufficient number of lines-of-sight. We found 3072 line-of-sight (each one box side in length) was adequate for the PM simulations, and 1024 for HPM.

## 3 Convergence Tests

Simulations of the Ly$`\alpha `$forest are restricted by two size constraints. High spatial resolution is required to resolve the structure of the absorbers, while a large box size is required both to capture the large-scale power which will affect the velocity widths of the absorption lines as well as to obtain a fair sample of the universe. Although a small box may be adequate to resolve the line structure, if the box is too small the fluctuations corresponding to the scale of the box will become nonlinear, in which case the results are no longer representative of the cosmological model simulated.

Resolution effects were previously studied using full hydrodynamics by Theuns et al. (1998b) and Bryan et al. (1999). Using P3M-SPH, Theuns et al. found for SCDM that a (comoving) box size of at least 5.5 Mpc using $`64^3`$ SPH particles and $`64^3`$ dark matter particles is required to converge on $`\overline{\tau }_\alpha `$ and the Doppler parameter distributions at $`z=3`$. This corresponds to a mean proper interparticle separation of 21 kpc, or about half this in structures with an overdensity of $`10`$, typical of the filaments with $`\mathrm{I}`$column densities of $`10^{14}10^{15}`$ cm-2 (Zhang et al. 1998). Bryan et al. found that a box size of at least 9.6 Mpc was needed to reach convergence on the Doppler parameter distribution at $`z=3`$ using a PM-PPM code (Kronus) with $`128^3`$ grid cells, corresponding to a proper resolution of 19 kpc.

Neither group was able to maintain a resolution of 10 kpc while increasing box sizes beyond 5 Mpc due to memory limitations, so that the contribution of larger scale power at high resolution was left unexplored. We also note that SCDM has less large-scale power than more viable cosmological models like $`\mathrm{\Lambda }`$CDM , so that the resolution limitations are least restrictive for SCDM. It is possible to achieve much higher resolution using PM or HPM. In this section we present the results of a series of tests of the convergence of the PM and HPM simulations for $`\mathrm{\Lambda }`$CDM as the box size, particle number, and number of grid zones are varied. The simulation results are all shown at $`z=3`$.

### 3.1 Flux distribution

A comparison of the cumulative distributions of flux per pixel for several resolutions and box sizes is shown in Figure 2. The PM-generated distributions agree at low flux values (high optical depths), but diverge for fluxes above 0.8, with differences in the cumulative distributions as large as 0.05. As the resolution of the $`9.6h^1\mathrm{Mpc}`$ box is increased to $`(N_p,N_g)=(256,512)`$, the amount of fine-scale structure continues to increase. The HPM-generated flux distributions differ markedly from the PM results, showing a smaller amount of fine-scale structure. The corresponding flux power spectra $`P_F(k)`$ are show in Figure 3. Again, as the grid resolution is increased in the $`9.6h^1\mathrm{Mpc}`$ box, the amount of large-scale power continues to increase, while there is a decrease in the amount of power at shorter scales ($`k>0.04h\mathrm{Mpc}^1`$). The $`9.6h^1\mathrm{Mpc}`$ HPM simulations show a deficit in power compared with the $`9.6h^1\mathrm{Mpc}`$ PM simulations for $`k>0.08h\mathrm{Mpc}^1`$, suggesting that the gas pressure forces have suppressed power at short scales.

### 3.2 Wavelet decomposition

The cumulative distributions of wavelet coefficients are shown in Figure 4. The groups of curve from left to right correspond to the respective velocity scales $`48`$ km s-1 , $`816`$ km s-1 , $`1632`$ km s-1 , $`3264`$ km s-1 , and $`64128`$ km s-1 . The small coefficients at the lowest scales show that no features are present on these scales. The physically most interesting scales are represented by the $`1632`$ km s-1 curves, corresponding to typical velocity widths of the absorption lines. The largest differences between the simulations are found on these scales. The convergence for the $`9.6h^1\mathrm{Mpc}`$ box simulations is non-monotonic with $`N_p`$ and $`N_g`$, with the $`(N_p,N_g)=(128,256)`$ and $`(256,512)`$ simulations agreeing well, but not with the $`(256,256)`$ simulation. The latter however agrees closely with the $`(256,512)`$ $`12.48h^1\mathrm{Mpc}`$ simulation, suggesting that convergence has been reached with box size. The largest coefficients are found for the smallest box ($`6.24h^1\mathrm{Mpc}`$) simulation, indicating a greater amount of structure in the spectra on these velocity scales than found for the larger boxes. By contrast, the HPM simulation shows the least amount of structure on these scales.

### 3.3 Line parameter decomposition

A comparison of the cumulative distributions of $`N_{\mathrm{HI}}`$ and $`b`$ determined by AutoVP is shown in Figure 5 for PM simulations for 3 different pairs of $`N_p`$ and $`N_g`$ and 3 box sizes. As was found for the wavelet coefficient distributions, the behaviours for both the $`N_{\mathrm{HI}}`$ and b distributions are non-monotonic. The $`(N_p,N_g)=(128,256)`$ and $`(256,512)`$ $`9.6h^1`$ Mpc simulations agree well in the b-distribution, but disagree with the $`(N_p,N_g)=(256,256)`$ simulations. It would appear that AutoVP is picking up unreal features resulting from noise in the $`N_p=128`$ simulations that decrease in going to $`N_p=256`$. Increasing the grid resolution, holding $`N_p=256`$ fixed, however, appears to recover real features that correspond statistically to the previously unresolved fluctuations in the $`(N_p,N_g)=(128,256)`$ simulations. Since the same trend is found for the wavelet coefficients, it seems unlikely that the non-monotonic behaviour is an artefact of the line finding and fitting procedures used by AutoVP. There are no obvious systematic trends with box size, the results otherwise appearing to have converged to better than 5%.

Also shown are the results for HPM for the $`9.6h^1\mathrm{Mpc}`$ box and $`(N_p,N_g)=(256,512)`$. Both the column density distribution and the Doppler parameter distribution agree well with the $`N_p=N_g=256`$ $`9.6h^1\mathrm{Mpc}`$ PM simulation, demonstrating that incorporating the pseudo-pressure forces in the HPM simulation has very little effect on the resulting absorption line parameter distributions for this model. We note that the HPM simulation tends to produce distributions shifted somewhat to higher column density and Doppler parameter compared with the analogous $`(N_p,N_g)=(256,512)`$ PM simulation.

A comparison of the results for HPM for different box sizes and particle/grid numbers is shown in Figure 6. The $`N_{\mathrm{HI}}`$ distributions are well-converged. The $`b`$-distributions would have appeared to converge, except the distributions shift toward larger $`b`$-values in the highest resolution, small box simulation. Curiously, this is opposite the trends reported in the resolution studies of Theuns et al. and Bryan et al. , who found that the $`b`$-values tended to decrease both with decreasing box size and increasing resolution, although at lower resolutions compared with our simulations. The trends found in the full hydrodynamics simulations are physically sensible: a larger box adds longer wavelength power to the simulations, which may increase the peculiar velocities of the dark matter and the gas, while higher particle and grid numbers permit higher resolution of density fluctuations and shocks, which may result in a narrowing of absorption features. Possibly the opposing trend found here is a consequence of not allowing for shocks. For a typical sound speed of $`10`$ km s-1 and peculiar velocities in the gas of $`100`$ km s-1 , shocks may be expected to play an important role in the velocity structure of the lines.

## 4 Comparison with Hydrodynamical Simulations

In this section we compare the results of the PM and HPM simulations against the Kronus simulation results of Machacek et al. (2000) (M00) for a variety of cosmological models. The simulations were run for the same cosmological models as in M00, matching the box sizes and the numbers of hydrodynamical cells and grid points for the gravitational force calculations.

### 4.1 Flux distribution

The distributions of flux for the 4 cosmological models considered are shown in Figure 7. Excellent agreement is found between HPM and the full hydrodynamical simulations (M00), with the agreement least good for SCDM. The agreement with PM is not as good, but the cumulative distributions still agree to within 10%. The agreement between both sets of $`N_p`$ and $`N_g`$ pairs for the HPM simulations shows that the results have converged at the relevant box sizes.

### 4.2 Wavelet decomposition

The cumulative distributions for the wavelet coefficients are shown in Figure 8. Much better agreement is generally found between the PM simulations and the full hydrodynamical simulations than is found for HPM. The agreement is even poorer for the HPM simulations with $`(N_p,N_g)=(256,512)`$ in all cases. The smaller values for the coefficients in the HPM spectrum indicates that the artificial hydrodynamics in HPM suppresses the full amount of velocity structure in the spectra compared with that found in the full hydrodynamics calculations.

### 4.3 Line parameter decomposition

The results of the Voigt line parameter analyses using AutoVP are shown in Figures 9 to 12. Both the PM and HPM results agree reasonably well with the results of M00 for the $`\mathrm{I}`$column density distributions. The best agreement for the Doppler parameter distributions is found for the PM simulations: the HPM simulations produce features that are too broad. While the PM and HPM simulations agree with the results of M00 to $`510`$% on the $`N_{\mathrm{HI}}`$ distributions, the agreement between the Doppler parameter distributions is worse: $`10`$% for PM and $`20`$% for HPM. The disagreement between the Doppler parameter distributions tends to increase with the magnitude of the power spectra at $`z=3`$.

## 5 Summary

We have performed several simulations of the Ly$`\alpha `$forest using different background cosmological models, numerical codes and grid resolutions. Convergence to the distributions of flux, wavelet coefficients, and absorption line parameters tends to be non-monotonic for the PM simulations as the number of particles and grid zones is increased, although agreements in the cumulative distributions to within 10% are achieved.

A comparison with full hydrodynamical calculations shows that there is little advantage in HPM over PM, except possibly for the flux distributions. The HPM simulations generally result in broader absorption lines than PM. Both the PM and HPM simulations tend to yield broader lines than found in the full hydrodynamics calculations. This may be a consequence of the neglect of shocks in the PM and HPM simulations, which may sharpen the velocity structure of the absorption lines.

The simulations show that the principal spectral properties of the Ly$`\alpha `$forest may be qualitatively reproduced from the dark matter structures alone, to a remarkable accuracy of 10% in the cumulative distributions of flux, wavelet coefficients, $`\mathrm{I}`$column density, and Doppler parameter. Tests of the models, however, require agreement at the 5% level per spectrum or better at the resolution of the Keck HIRES (Meiksin et al. 2000), so that ultimately full hydrodynamical calculations will be required to confront models for the structure and physical properties of the Ly$`\alpha `$forest with direct observations. In order to match the observed spectra, it may ultimately be necessary to incorporate radiative transfer into the simulations both to allow for spatial fluctuations in the UV ionizing background and to reproduce the correct gas temperatures (and line widths), particularly if helium is completely photoionized at only moderate redshifts. Several other physical effects may also be important, such as heating sources in addition to QSOs or galactic feedback like winds. These effects may have sufficient impact on the statistics of the Ly$`\alpha `$forest that 10% accuracy is adequate to constrain them. Our results suggest that PM simulations alone may provide a reliable tool for investigating some of these issues without the considerable overhead of solving the full set of hydrodynamics equations.

## acknowledgments

The authors thank G. Bryan for permission to use the results of his Kronos simulations, and R. Davé for permission to use AutoVP. M. White was supported by the US National Science Foundation and a Sloan Fellowship. Parts of this work were done on the Origin2000 system at the National Center for Supercomputing Applications, University of Illinois, Urbana-Champaign.


text:
# Resolving Gamma-Ray Burst 000301C with a Gravitational Microlens

## 1 INTRODUCTION

The rapid localization of gamma-ray bursts (GRBs) has brought a new dimension to GRB research by allowing many events to be followed up at longer wavelengths. Afterglows of bursts have been detected at X-ray (Costa et al. 1997), optical (van Paradijs et al. 1997) and radio (Frail et al. 1997) wavelengths. Precise positions have allowed redshifts to be measured for a number of GRBs (Metzger et al. 1997), providing a definitive proof of their cosmological origin.

The afterglow of a GRB is thought to be synchrotron radiation from a relativistic shock driven into the circumburst environment (Meszaros & Rees 1993, 1997; Paczynski & Rhoads 1994; Katz 1994; Waxman 1997a). The light curves and spectral energy distributions are well fit by power-laws as expected from the shock model. There is recent evidence that at least some GRB are not spherical explosions. A broad-band break in the light curve power-law index was predicted for shocks produced by collimated jets (Rhoads 1997) and such breaks have been seen in GRB 990510 (Stanek et al. 1999; Harrison et al. 1999), GRB 991216 (Halpern et al. 2000), and GRB 000301C (Sagar et al. 2000; Masetti et al. 2000; Jensen et al. 2000; Berger et al. 2000). The ratio of the spectral index to the light curve index also suggests non-spherical energy ejection for some events (e.g. GRB 991216: Garnavich et al. 2000a). In general, the synchrotron afterglow model has been very successful in matching most of the observations. However, it is heavily strained explaining the well-studied afterglow of GRB 000301C which shows a peculiar achromatic fluctuation which deviates significantly from the broken power-law fit to the lightcurve (Sagar et al. 2000; Berger et al. 2000).

Here we propose an elegant solution to the GRB 000301C problem: the GRB was microlensed. Previously, Loeb & Perna (1998) predicted that microlensing by a solar mass lens at a cosmological distance would produce a nearly achromatic fluctuation of similar amplitude and duration to that observed. In §2 we describe the optical and near-IR observations available for our analysis. In §3 we describe the simplest parameterization of a GRB microlensing event and discuss our technique in fitting the data. In §4 we analyze the results and estimate the probability that a microlensing event can occur given the observational constraints. Finally, we summarize our conclusions in §5.

## 2 THE BURST

GRB 000301C was detected by ASM, Ulysses and NEAR on 2000 March 1.41085 (UT) and localized to a region of about 50 arcmin2 on the sky (Smith, Hurley & Cline 2000). An optical transient was detected about $`1.5`$days after the burst (Fynbo et al. 2000) and its presence confirmed in the near-IR (Stecklum et al. 2000) and at radio wavelengths (Bertoldi 2000). The initial optical decline was steep with a power-law index of $`\alpha =1.6\pm 0.3`$ (Halpern, Mirabal & Lawrence 2000), but by $`3`$days after the burst the decay rate slowed down (Garnavich et al. 2000b), a very unusual occurrence in an optical afterglow. This “standstill” did not last long and two weeks after the burst the power-law index was $`\alpha =2.7`$ (Veillet 2000b). Rhoads & Fruchter (2000) found that their near-IR data behaved differently than the optical although the time sampling was more sparse in the IR.

Spectra of the afterglow showed a strong Lyman cutoff in the UV (Smette et al. 2000) and absorption features of Ly$`\alpha `$ and metal lines were consistent with a large absorbing column at $`z=2.04`$ (Jensen et al. 2000; Feng, Wang & Wheeler 2000). This could be gas in the host or an intervening galaxy, but the lack of a Lyman break longward of 318 nm (Feng et al. 2000) places an upper limit on the redshift of the GRB of $`z_\mathrm{s}<2.5`$.

The most complete photometric compilation is by Sagar et al. (2000) (shown in their Fig. 3) and is roughly described by two power-laws with a break around $`7`$days after the burst. But there are clearly significant deviations from the broken power-law model, especially between 3 and $`6`$days after the burst. A similar break was seen in GRB 990510 (e.g. Stanek et al. 1999) and interpreted as the lateral spreading of the jet as its Lorentz factor dropped below the inverse of its opening angle. But GRB 000301C is unique in showing variability on the relatively short time-scales of hours to a few days superimposed on the more typical power-law trends. As noted by Berger et al. (2000) and Sagar et al. (2000), the rapid fluctuation is achromatic over a range of more than a factor of 5 in wavelength. Achromatic brightening of a source can be a sign of gravitational lensing, but the time-scale for microlensing a point source at $`z>2`$ is substantially longer than a few days due to the low velocities ($`10^3c`$) of conventional sources and lenses. However, Loeb & Perna (1998) showed that the superluminal motion of a GRB source on the sky generically results in microlensing events with durations of hours to days.

For our analysis, we use the compilation of Sagar et al. (2000) with additional $`UBVI`$ photometry (one point in each band) from Stanek et al. (2000). One $`R`$-band point has been removed because of its large deviation from other data taken near the same time. This leaves a total of 104 photometric data points in seven photometric bands, with the following distribution of points: $`N(U,B,V,R,I,J,K)=(6,18,8,46,16,3,7)`$.

## 3 ANALYSIS

### 3.1 Microlensing of a GRB

A spherical GRB fireball appears on the sky as a thin ring that exhibits superluminal expansion at a speed $`\gamma c`$, where $`\gamma `$ is the Lorentz factor of the GRB shock (Waxman 1997b). Figure 3 of Sagar et al. (2000) shows that the achromatic fluctuation in the lightcurve of GRB 0000301C occurred well before the power-law break, and so we assume that the during the microlensing event the jet still behaves as if it is part of a spherical fireball.

We use $`t`$ to denote the observed time since the GRB trigger in units of days. There are two sets of parameters that will define the light curves. The first set describes the intrinsic GRB light decays: power-law slopes, $`\alpha _1`$ and $`\alpha _2`$, parameter $`\beta `$ that describes the smoothness of the transition between them, and the transition time $`t_b`$. Each wavelength band also requires a parameter that sets the zero point of the power-law.

Loeb & Perna (1998) showed that the simplest microlensing model can be described by three parameters (which are all constants), namely: $`R_0`$, $`b`$ and $`W`$. They are defined as follows:

(i) $`R_0`$ is the ring radius at $`t=1\mathrm{day}`$ in units of the Einstein radius of the lens. At other times the ring radius evolves as (Waxman 1997b)

$$R_s(t)=R_0t^{5/8}.$$
(1)
Here $`R_0=\rho _0/r_E`$, and

$$\rho _0=4\times 10^{16}\left(\frac{E_{53}}{n_1}\right)^{\frac{1}{8}}(1+z_\mathrm{s})^{5/8}\mathrm{cm},$$
(2)
where the factor involving the source redshift $`z_\mathrm{s}`$ is due to the cosmic time dilation, $`E_{53}`$ is the “isotropic-equivalent” of the energy release in units of $`10^{53}\mathrm{erg}\mathrm{s}^1`$, and $`n_1`$ is the ambient gas density in units of $`1\mathrm{cm}^3`$. The actual energy release could, of course, be much smaller due to the small solid angle occupied by the jet, but this does not affect the source size until the break in the lightcurve. The Einstein radius of a lens of mass $`M_{\mathrm{lens}}`$ is,

$$r_E=\left[\left(\frac{4GM_{\mathrm{lens}}}{c^2}\right)D\right]^{1/2}=7.7\times 10^{16}\left(\frac{M_{\mathrm{lens}}}{1M_{}}\right)^{1/2}\left(\frac{D}{10^{28}\mathrm{cm}}\right)^{1/2}\mathrm{cm}$$
(3)
where $`D(D_\mathrm{s}D_{\mathrm{ls}}/D_\mathrm{l})`$ is the ratio of the angular-diameter distances between the observer and the source, the lens and the source, and the lens and the observer. The value of $`D`$ depends on the lens and source redshifts and the cosmological parameters.

(ii) $`b`$ is the lens-source separation on the sky (or equivalently, the “impact parameter”) in units of the Einstein angle, $`\theta _E(r_E/D_\mathrm{s})`$. Since the apparent source radius is changing faster than the speed of light, all other non-relativistic astrophysical motions are irrelevant. We may therefore consider the source-center and lens positions as fixed on the sky. Initially, as long as $`R_\mathrm{s}b`$ the source is pointlike and is magnified by a time-independent factor. The maximum in the magnification curve is reached when $`R_\mathrm{s}b`$; the timing of this maximum can be used to fix the ratio $`R_0/b`$.

(iii) $`W`$ is the brightness-weighted width of the ring divided by its radius. The width determines the height and the duration of the magnification event. The larger the width, the smaller the height is and the broader the duration is. This width was estimated theoretically to be $`10\%`$ in the optical-IR (see, e.g. Waxman 1997b; Sari 1998; Panaitescu & Meszaros 1998; Granot et al. 1999) but is still subject to uncertainties concerning the extent of the emission region behind the GRB shock front. Note that the assumption that the ring has sharp boundaries and that $`R_0`$ and $`W`$ are the same for all wavelengths, is an over-simplification (see ring profiles in Granot et al. 1999). Ring properties that vary significantly with wavelength would produce different light curves in the various observed bands (Wambsganss & Paczyński 1991). But, more parameters are required to fit a more elaborate model and this is not justified by the quality of the available data on GRB 000301C.

If $`W`$ were close to unity, then the lensing signal would be very weak. Thus, GRB 000301C provides the first evidence that a GRB produces a ring on the sky. Previous data, such as the scintillations of GRB 970508, constrained the source size but not its shape (Waxman, Kulkarni, & Frail 1998).

### 3.2 Magnification Factor

The magnification factor $`\mu `$ is a function of $`R_s`$, $`W`$, and $`b`$. As described by Loeb & Perna (1998),

$$\mu (R_s,W,b)=\frac{\mathrm{\Psi }[R_s,b](1W)^2\mathrm{\Psi }[(1W)R_s,b]}{1(1W)^2},$$
(4)
where $`\mathrm{\Psi }(R_s,b)`$ is the magnification for a uniform disk of radius $`R_s`$ (Schneider, Falco, & Ehlers 1992),

$$\mathrm{\Psi }[R_s,b]=\frac{2}{\pi R_s^2}\left[_{|bR_s|}^{b+R_s}𝑑r\frac{r^2+2}{\sqrt{r^2+4}}\mathrm{arccos}\frac{b^2+r^2R_s^2}{2rb}+H(R_sb)\frac{\pi }{2}(R_sb)\sqrt{(R_sb)^2+4}\right].$$
(5)
Here $`H(x)`$ is the Heaviside step function. The integral in equation (5) can be expressed more explicitly as a sum of elliptic integrals (Witt & Mao 1994). In general, an arbitrary ring profile can be incorporated as a sum over a set of infinitesimal rings.

### 3.3 Fitting the Data

Using all the data for GRB 000301C (104 points), we performed $`\chi ^2`$ minimization fits using broken power-law model (11 free parameters: 4 shape parameters $`+`$ 7 photometric zero points), in the form described by Sagar et al. (2000), and using broken power-law plus microlensing model (14 free parameters: 11 parameters $`+`$ 3 microlensing parameters). We repeated the calculation using only the $`R`$-band data (5 and 8 free parameters, respectively), since it had the most complete temporal coverage with the total of 46 points.

For the broken power-law model, the fits were clearly not good, with $`\chi ^2/DOF=2.76`$ for all seven bands and $`\chi ^2/DOF=2.99`$ for the $`R`$-band data only, values similar to those found in the other papers describing GRB 000301C (e.g. Jensen et al. 2000). Adding the microlensing magnification has substantially improved the fits, with $`\chi ^2/DOF=1.77`$ for all seven bands and $`\chi ^2/DOF=1.48`$ for the $`R`$-band data only.

The best-fit microlensing models are shown in Figure 1 for all seven bands. Also shown, in an insert, is the best-fit model for the $`R`$-band data only. The best-fit to all data provides an impact parameter $`b=1.04\pm 0.02`$ and a fractional ring width of $`W=0.16\pm 0.02`$. The best-fit to the $`R`$-band data only yields an impact parameter $`b=1.09\pm 0.02`$ and a substantially smaller width of the ring, $`W=0.07\pm 0.02`$. In both cases $`R_0=0.49\pm 0.02`$, and microlensing provides an angular resolution better than a micro-arcsecond in probing the GRB fireball. The values for the broken power-law parameters obtained by us are very similar to those obtained by Sagar et al. (2000): $`\alpha _1=1.1,\alpha _2=2.9`$ and $`t_b=7.6`$days, regardless if all the data or $`R`$-band only were fitted. The value of $`\beta `$ is not well constrained, but it is always large, $`\beta >5`$, indicating a sharp break. The effect of lensing is shown in Figure 2 after our best fit broken power-law is removed.

The errors in the microlensing parameters estimations are based on conditional probability distributions, obtained by fixing the rest of the parameters at their most probable values, and should be treated only as rough estimates. It is clear from Figure 1 that the actual errors of the photometry are not Gaussian, but are dominated by systematic errors, probably resulting from different reduction procedures applied to the various data. Given the special nature of GRB 000301C, it would be well-worth the effort for one of the groups, possibly one which obtained much of the photometric data for this afterglow, to reduce the CCD data obtained by the other groups in order to ensure a uniform reduction procedure. Such uniform reduction would allow for better statistical testing of the GRB 000301C microlensing hypothesis.

## 4 DISCUSSION

The global fit to the data matches the light curves of the individual bands well. The fit to the $`R`$-band alone is also of high quality and gives similar values for the microlensing parameters. The main difference between the global and $`R`$-band is in the sharpness of the microlensing peak which is narrower and therefore favors smaller values of $`W`$ in the $`R`$-band.

We can estimate the mass of the lens from our measurement of the ring radius parameter, $`R_0`$. The lens mass can be written as

$$M_{\mathrm{lens}}=0.13M_{}\times \left(\frac{E_{53}}{n_1}\right)^{1/4}\left(\frac{D_\mathrm{l}}{D_{\mathrm{ls}}}\right)R_0^2,$$
(6)
assuming a source redshift of $`z_\mathrm{s}=2.04`$ and a flat universe with $`h=0.7`$ and $`\mathrm{\Omega }_m=0.3`$. Hence, for $`R_0=0.5`$, the lens mass is of order $`0.5M_{}`$ i.e. stellar size, for a wide range of reasonable GRB energies and circumburst densities. This is reassuring for the microlensing interpretation since it does not rely on some peculiar and rare object to do the lensing. In a case analogous to GRB sources, Koopmans et al. (2000) have recently interpreted variability of the macro-lensed superluminal radio source B1600+434 as being due to microlensing.

### 4.1 Probability for Microlensing

Deep Hubble Space Telescope (HST) imaging of the GRB and field were obtained by Fruchter, Metzger & Petro (2000) but no host or intervening galaxy was detected to a limit of 28.5 mag. A galaxy with $`R=24.3\pm 0.3`$ mag (Veillet 2000a) is 2′′ from the position of the GRB. Our analysis of the deep STIS exposures obtained 2000, April 19, is shown in Figure 3 and suggests that the surface brightness profile along the major axis of the galaxy is best fit by an exponential with a scale length of $`0.22^{\prime \prime }`$. This implies that the light is from a rather compact disk. The redshift of this nearby galaxy is unknown. The galaxy is not detectable out to the position of the GRB and extrapolation of the surface brightness profile implies a very low surface brightness of $`<10^{23}`$ $`\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1\mathrm{\AA }^1\mathrm{arcsec}^2`$ at an effective wavelength of 585 nm. Of course, a dark halo may extend well out beyond the visible disk. From the background variations in the $`HST`$ images, we place a limit on the surface brightness at the GRB of $`>28.5`$ $`\mathrm{mag}\mathrm{arcsec}^2`$ in approximately the $`V`$-band.

What is the probability for microlensing by a solar mass star at a cosmological distance? This probability can be written as $`P_{\mathrm{ml}}=\mathrm{\Sigma }_{}/\mathrm{\Sigma }_{\mathrm{crit}}`$, where $`\mathrm{\Sigma }_{}`$ is the surface mass density of stars in an intervening galaxy close to the line-of-sight and $`\mathrm{\Sigma }_{\mathrm{crit}}=(c^2/4\pi G)(D_\mathrm{s}/D_\mathrm{l}D_{\mathrm{ls}})b^2`$. If the mass-to-light ratio of these stars is $`(M/L)`$, then their apparent surface brightness is $`\mu _{\mathrm{SB}}=(M/L)^1(\mathrm{\Sigma }_{}/4\pi (1+z_\mathrm{l})^4)`$. Hence the expected surface brightness of stars around the location of the afterglow is given by,

$$\mu _{\mathrm{SB}}=\left(\frac{c^2}{16\pi ^2Gb^2}\right)\left(\frac{M}{L}\right)^1\left(\frac{D_\mathrm{s}}{D_\mathrm{l}D_{\mathrm{ls}}}\right)\frac{P_{\mathrm{ml}}}{(1+z_\mathrm{l})^4}.$$
(7)
The simple point-lens model is adequate for $`P_{\mathrm{ml}}1`$, so that the caustics induced by the external shear occupy a region much smaller than $`r_E`$ (Chang & Refsdal 1984). For $`P_{\mathrm{ml}}0.1`$, $`z_\mathrm{s}=2`$, $`b=1`$, $`(M/L)=5`$ in solar units, and a flat universe with $`h=0.7`$ and $`\mathrm{\Omega }_m=0.3`$, we find $`\mu _{SB}29`$ magnitudes per arcsec2 in the V-band, for a lens galaxy at $`1z_\mathrm{l}1.7`$, assuming a $`K`$-correction of $`4`$ V-magnitudes as required for an elliptical galaxy in this redshift interval (see Fig. 10 of Fukugita et al. 1995). This value of $`\mu _{SB}`$ is below the inferred upper limit on the surface brightness from the HST image. Coincidentally, it is comparable to the $`R^{1/4}`$–law extrapolation of the surface brightness in Figure 3, but well above the exponential extrapolation. The above constraint is much weaker if the lensing star belongs to a halo population of compact objects (MACHOs) which have a high $`M/L`$ ratio. The microlensing probability would inevitably be large if the HI column density $`10^{21}\mathrm{cm}^2`$, detected at $`z=2`$ (Jensen et al. 2000; Smette et al. 2000), is due to the galactic host of the lens (Perna & Loeb 1997); however the observed damped Ly$`\alpha `$ absorber is more likely to be the host galaxy of the GRB.

### 4.2 Alternative Interpretations

Fluctuations in the lightcurve may result from shell collisions within the fireball (Dai & Lu 2000), but such collisions would re-energize the fireball and change its spectrum by increasing the peak frequency of its emission. There is no evidence for a chromatic change of this type during the event (Berger et al 2000). The activity of the central engine of the GRB has to be fine-tuned in an ad-hoc manner so as to produce a collision only after $`4`$days when the Lorentz factor of the expanding shell ($`\gamma 5`$) already declined by more than an order of magnitude relative to its initial value.

Berger et al. (2000) suggested an alternative interpretation of the achromatic fluctuation in terms of an inhomogeneity of the ambient density into which the fireball is propagating, $`n`$. Naively, one may argue that the flux from a fireball propagating into a uniform ambient medium scales as $`n^{1/2}`$ at all wavelengths, and so a brightening of the flux by a factor of 1.7 may be achieved if the fireball encounters a clump of gas which is a factor of $`3`$ denser than the mean. The required clump location ($`5\times 10^{17}\mathrm{cm}`$), transverse size ($`10^{17}\mathrm{cm}`$), depth ($`5\times 10^{17}\mathrm{cm}`$) and overdensity ($`3`$), can all be chosen ad-hoc so as to fit the data. However, the standard model for the fireball dynamics and emission implies that the flux at a particular observed time is emitted with different weights on a particular spatial region for frequencies above and below the peak spectral frequency \[see Eqs. (15) and (16) in Wang & Loeb 1999\]. This results in two effects: (i) the ring is narrower and has a higher contrast (center to limb variation in brightness) at higher frequencies \[see Figs. (11) and (12) in Granot et al. 1999\]; (ii) the limb of the ring is influenced by the clump region at a later observed time than its center because of the geometric time delay. Since the peak frequency is $`300`$GHz at $`4`$days (Berger et al. 2000; Sagar et al. 2000), the fluctuation is therefore expected to brighten to a maximum earlier in the radio than in the optical-IR. The chromaticity should be even more pronounced below the synchrotron self-absorption frequency (Berger et al. 2000), as this frequency depends on $`n`$. Although it may not be possible to rule out the above chromatic behavior with the sparse radio data for GRB 000301C, achromaticity does not trivially follow from this model at wavelengths ranging all the way from the optical to the radio.

The ring parameters $`W`$ and $`R_0`$ which define the microlensing lightcurves should also depend on wavelength (Panaitescu & Meszaros 1998; Granot et al. 1999) but in a generic predictable way. Hence, with better analysis of the existing data for GRB 000301C or with detailed monitoring of future GRBs it should be possible to test it against the above interpretations.

## 5 CONCLUSIONS

We successfully model the light curve of GRB 000301C as a standard broken power-law plus a gravitational microlensing event. We find the lens mass required for the event is $`0.5`$M, thus requiring only a normal star to explain the lensing. No galaxy or diffuse light from the stellar population is detected along the line of sight, but reasonable lensing probability is expected even for surface brightness below the deep $`HST`$ images we analyzed. Ground based images with a large telescope may detect the lensing population.

New GRB satellites such as $`HETE2`$ and $`SWIFT`$ will provide many opportunities to study GRB afterglows. This is the first of what will be a number of microlensed GRBs and we show that analysis of high-quality data covering a large wavelength range can be a useful tool for probing the physical structure of GRB afterglows with sub micro-arcsecond resolution. Our estimate of $`W`$ confirms the prediction that afterglow shocks appear as thin rings on the sky (Waxman 1997b; Sari 1998; Panaitescu & Meszaros 1998; Granot et al. 1999). Better broad-band data could check for variation with wavelength of the ring width and radius, and polarimetric observations could confirm the predictions of the polarization variations made by Loeb & Perna (1998). Future observations could also constrain the properties of the dark matter by taking a census of the number and masses of microlenses.

We thank Bohdan Paczyński and Eli Waxman for useful comments on the manuscript and Dave Bennett for helpful discussions. Support for KZS was provided by NASA through Hubble Fellowship grant HF-01124.01-A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This work was supported in part by NSF grant AST-9900877, the Israel-US BSF, and the NASA grant NAG5-7039 for AL. Support for PMG was provided by NASA LTSA grant NAG5-9364.


text:
# Legendrian Mirrors and Legendrian Isotopy

## Acknowledgments

The author would like to thank Tom Mrowka and Isadore Singer for their encouragement and support.