text:
# An Improved Exact Riemann Solver for Relativistic Hydrodynamics

## 1 Introduction

The general Riemann problem is based on the calculation of the one-dimensional temporal evolution of a fluid which, at some given initial time, has two adjacent states characterized by different values of uniform velocity, pressure and density. From a mathematical point of view, the Riemann problem is an initial value problem for a hyperbolic system of partial differential equations, with initial conditions characterized by a discontinuity between the two initial states. These initial conditions establish the way in which the discontinuity will decay after removal of the barrier separating the initial “left” and “right” states. The schematic evolution of a general Riemann problem can be represented as (Martí & Müller 1994)

$$L𝒲_{}L_{}𝒞R_{}𝒲_{}R,$$
(1.1)
where $`𝒲`$ denotes a shock or a rarefaction wave that propagates towards the left $`()`$ or the right $`()`$ with respect to the initial discontinuity, $`L`$ and $`R`$ are the known initial left and right states, while $`L_{}`$ and $`R_{}`$ are the new hydrodynamic states that form behind the two waves propagating in opposite directions. These waves are separated by a contact discontinuity $`𝒞`$ and therefore have the same values of the pressure and velocity, but different values of the density∗\*∗\*$``$Note that the contact discontinuity might be also trivial so that the density may not be discontinuous.. The Riemann problem is said to be solved when the velocity, pressure and density in the new states $`L_{}`$ and $`R_{}`$ have been computed, as well as the positions of the waves separating the four states. The solution of the one dimensional Riemann problem in relativistic hydrodynamics was discussed in the general case by Martí & Müller (1994) and the reader is referred to their work for further details (see also Pons et al. 2000, for the extension to multidimensions).

The numerical solution of a Riemann problem is the fundamental building block of hydrodynamical codes based on Godunov-type finite difference methods. In such methods, the computational domain is discretized and each interface between two adjacent grid-zones is used to construct the initial left and right states of a “local Riemann problem”. The evolution of the hydrodynamical equations is then obtained through the solution across the computational grid of the sequence of local Riemann problems set up at the interfaces between successive grid-zones (see Martí & Müller 1996, but also Godunov 1959 and Colella & Woodward 1984). The “core” of each of these local Riemann problems consists of determining the fluid pressure across the contact discontinuity, which can be calculated by imposing the continuity of the normal component of the fluid velocity across $`𝒞`$

$$v_L_{}(p_{})=v_R_{}(p_{}).$$
(1.2)
In general, (1.2) is a nonlinear algebraic equation in the unknown pressure $`p_{}`$ and requires a numerical solution. Depending on the different wave patterns forming after the decay of the discontinuity, a different nonlinear equation will need to be solved††††††$`\mathrm{}`$This approach is usually referred to as an “exact” Riemann solver to distinguish it from the family of so called “approximate” Riemann solvers, where the system of equations to be solved is reduced to quasi-linear form, thus avoiding any iterative procedure. See, for example, the approximate Riemann solver of Roe (1981).. This initial “ambiguity” in the wave pattern produced corresponds to the fact that the interval in pressure bracketing the solution $`p_{}`$ is not known a priori. In practice this lack of information is compensated by the use of efficient numerical algorithms which, via a process of trial and error, determine the correct wave pattern and then proceed to the solution of the corresponding nonlinear equation (Martí, & Müller 1999).

In this paper, we show that the relativistic expression for the relative velocity between the two initial states is a function of the unknown pressure $`p_{}`$ and so a new procedure for numerically solving the exact Riemann problem can be proposed in which the pressure $`p_{}`$ is no longer obtained by the solution of equation (1.2). Rather, $`p_{}`$ is calculated by equating the relativistic invariant expression for the relative velocity between the two initial states with the value given by the initial conditions.

When compared to equivalent approaches, our exact Riemann solver has some advantages. Firstly, we can remove the ambiguity mentioned above and determine the generated flow pattern by simply comparing the relative velocity between the two initial states with reference values built from the initial conditions of the Riemann problem. Doing so provides immediate information about which of the nonlinear equations (one for every wave pattern) needs to be solved. Secondly, by knowing the wave pattern we can produce an immediate bracketing of the solution. Doing so gives improved efficiency in the numerical root finding procedure. Finally, for one of the wave patterns (i.e. for two rarefaction waves moving in opposite directions) our method provides the solution of the relativistic Riemann problem in a closed analytic form.

Because of its simplicity, the numerical implementation of our method is straightforward and can be accomplished with a much smaller number of lines of code. When compared with other exact Riemann solvers (e. g. Martí & Müller 1999) it has also proved to be computationally more efficient. In particular, when solving a generic hydrodynamical problem (in which one solves for very simple Riemann problems) the approach proposed here brackets the solution very closely and this produces substantial computational improvements of up to $`30\%`$. For the cases discussed here however, where very strong shocks are considered, the speed-up is below $`10\%`$.

This paper briefly introduces our idea and is organized as follows: in Section 2 we define the mathematical set-up for the formulation of the Riemann problem in relativistic hydrodynamics. In Section 3 we determine the basic relativistic expressions linking the velocities ahead of and behind a shock or a rarefaction wave. These expressions will be used repeatedly in the following Sections 46, where will write explicit criteria for the occurrence of the three possible wave patterns. Section 7 discusses how the criteria can be used for making an efficient numerical implementation of a Riemann solver and briefly presents a comparison of performances with more traditional algorithms. Finally, Section 8 presents our conclusions. Throughout, we use a system of units in which $`c=1`$.

## 2 Setting up the Relativistic Riemann Problem

Consider a perfect fluid with four-velocity $`u^\mu =W(1,v,0,0)`$ with $`W(1v^2)^{1/2}`$ being the Lorentz factor and with a stress-energy tensor

$$T^{\mu \nu }=(e+p)u^\mu u^\nu +p\eta ^{\mu \nu }=\rho hu^\mu u^\nu +p\eta ^{\mu \nu },$$
(2.3)
where $`\mu =0,\mathrm{},3`$, $`\eta ^{\mu \nu }=\mathrm{diag}(1,1,1,1)`$ and $`e`$, $`p`$ are the proper energy density and pressure, respectively. Assume moreover that the fluid obeys a polytropic equation of state

$$p=k(s)\rho ^\gamma =(\gamma 1)\rho ϵ,$$
(2.4)
where $`\rho `$ is the proper rest mass density, $`\gamma `$ is the adiabatic index, and $`k(s)`$ is the polytropic constant, dependent only on the specific entropy $`s`$ (the latter is generally assumed to be different in the two initial states). Straightforward expressions can be written relating $`e`$ and $`p`$ to the specific enthalpy $`h`$ and to the specific internal energy $`ϵ`$ of the fluid

$`e`$ $`=`$ $`\rho (1+ϵ)=\rho +{\displaystyle \frac{p}{\gamma 1}},`$ (2.5)
$`h`$ $`=`$ $`1+ϵ+{\displaystyle \frac{p}{\rho }}=1+{\displaystyle \frac{p}{\rho }}\left({\displaystyle \frac{\gamma }{\gamma 1}}\right).`$ (2.6)
Consider now the fluid to consist of an initial left state (indicated with an index $`1`$) and an initial right state (indicated with an index $`2`$), each having prescribed and different values of uniform pressure, density and velocity‡‡‡‡‡‡$`\mathrm{}`$Note that hereafter we will consider the fluid to be of the same type in the two initial states. In principle, in fact, the fluid in the two initial states might be governed by two different equations of state or by polytropic equations of state with different polytropic indices. The formulation of the problem in that case is equivalent to the one presented here for a single type of fluid but particular attention must be paid to distinguishing the different components in the relevant expressions.. The discontinuity between the two states which has been constructed in this way will then decay producing one of the three wave patterns listed below

(i)  two shock waves, one moving towards the initial left state, and the other towards the initial right state: $`L𝒮_{}L_{}𝒞R_{}𝒮_{}R`$.

(ii)  one shock wave and one rarefaction wave, the shock moving towards the initial right state, and the rarefaction towards the initial left state (or viceversa if $`p_2>p_1`$): $`L_{}L_{}𝒞R_{}𝒮_{}R`$.

(iii)  two rarefaction waves, one moving towards the initial left state, and the other towards the initial right state: $`L_{}L_{}𝒞R_{}_{}R`$. A special case of this wave pattern is produced when the two rarefaction waves leave a vacuum region behind them.

The basis of our approach lies in the possibility of determining “a priori” which of these three wave patterns will actually result by simply comparing the value of the special relativistic relative velocity between the initial left and right states

$$v_{12}\frac{v_1v_2}{1v_1v_2},$$
(2.7)
with the values of the limiting relative velocities for the occurrence of the wave patterns (i)–(iii). In this respect, our approach represents the relativistic generalization of a similar analysis proposed in Newtonian hydrodynamics by Landau and Lifshitz (1987) (see also Gheller 1997 for a numerical implementation). Mathematically, the occurrence of the different cases can be distinguished because the relative velocity§§§§§§$`\mathrm{\S }`$For compactness we will hereafter refer to as “relative velocity” the relativistic invariant expression. (2.7) is a continuous monotonic function of the pressure $`p_{}`$ across the contact discontinuity. This is shown in Fig. 1 where we have plotted $`v_{12}`$ as a function of $`p_{}`$ and it is proved mathematically in Appendix A. Note that the curve shown in Fig. 1 is effectively produced by the continuous joining of three different curves describing the relativistic relative velocity for the three different wave patterns corresponding respectively to two shocks (2S, indicated as a dashed line), one shock and one rarefaction wave (SR, indicated as a dotted line), and two rarefaction waves (2R, indicated as a continuous line). The joining of the different curves occurs at precise values of the relative velocity which we have indicated as $`(\stackrel{~}{v}_{12})_{_{2S}}`$, $`(\stackrel{~}{v}_{12})_{_{SR}}`$, $`(\stackrel{~}{v}_{12})_{_{2R}}`$ and which depend uniquely on the initial conditions of the two unperturbed states. These three values of the relative velocity mark the extremes of three different intervals on the vertical axis and correspond to the three different cases (i)–(iii).

Within this framework then, it is sufficient to compare the relative velocity between the initial states (2.7) with the calculated values of the limiting relative velocities $`(\stackrel{~}{v}_{12})_{_{2S}}`$, $`(\stackrel{~}{v}_{12})_{_{SR}}`$, $`(\stackrel{~}{v}_{12})_{_{2R}}`$, to determine, in advance, which wave pattern will be produced. This, in turn, determines the correct equation to solve and the correct bracketing of the solution. In the following Sections we will discuss in detail how to derive analytic expressions for the limiting relative velocities and use them efficiently within a numerical code.

## 3 Limiting Relative Velocities: Fundamental Expressions

The expression for the relative velocity (2.7) is a relativistic invariant, but its calculation can be considerably simplified when performed in an appropriate reference frame. In practice, we will consider each of the three wave patterns (i)–(iii) as being composed of two “discontinuity fronts” (two shocks, two rarefaction waves, or one of each) moving in opposite directions and separated by a region where a contact discontinuity is present. In the case of a shock, in particular, it is useful to use a reference frame comoving with the shock front, in which the relativistic expression for the relative velocities ahead of ($`a`$) and behind ($`b`$) the shock takes the form (Taub 1978)

$$v_{ab}\frac{v_av_b}{1v_av_b}=\sqrt{\frac{(p_bp_a)(e_be_a)}{(e_a+p_b)(e_b+p_a)}}.$$
(3.8)
In the case of a rarefaction wave, on the other hand, it is more convenient to use the Eulerian frame in which the initial states are measured and which has one of the axes aligned with the direction of propagation of the wave front. In such a frame, the flow velocity at the back of a rarefaction wave (i.e. behind the tail of the rarefaction wave) can be expressed as a function of pressure at the back of the wave as

$$v_b=\frac{(1+v_a)A_\pm (p_b)(1v_a)}{(1+v_a)A_\pm (p_b)+(1v_a)}.$$
(3.9)
The quantity $`A_\pm (p)`$ in (3.9) is defined as (Martí & Müller 1994)

$$A_\pm (p)\left\{\left[\frac{(\gamma 1)^{1/2}c_s(p)}{(\gamma 1)^{1/2}+c_s(p)}\right]\left[\frac{(\gamma 1)^{1/2}+c_s(p_a)}{(\gamma 1)^{1/2}c_s(p_a)}\right]\right\}^{\pm 2/(\gamma 1)^{1/2}},$$
(3.10)
with the $`\pm `$ signs corresponding to rarefaction waves propagating to the left ($`_{}`$) and to the right ($`_{}`$) of the contact discontinuity, respectively. The quantity $`c_s(p)`$ in (3.10) is the local sound speed which, for a polytropic equation of state, can be written as

$$c_s=\sqrt{\frac{\gamma (\gamma 1)p}{(\gamma 1)\rho +\gamma p}}.$$
(3.11)
We can now use expression (3.9) to write an invariant expression for the relative velocity across a rarefaction wave (i.e. the relative velocity of the fluid ahead of the rarefaction wave and behind the tail of the rarefaction wave) as

$$v_{ab}=\frac{1A_\pm (p_b)}{1+A_\pm (p_b)}.$$
(3.12)
Expressions (3.8) and (3.12) are not yet in a useful form since they cannot be combined to give (2.7). However, we can exploit the fact they are relativistic invariants to evaluate them in the rest frame of the contact discontinuity. The latter is, by definition, comoving with the fluid behind the discontinuity fronts and within such a frame we can set $`v_b=0`$ and use equations (3.8), (3.12) to obtain an explicit expression for the velocity ahead of the wave front. The expressions for the velocities ahead of both discontinuity fronts obtained in this way can then be combined so as to evaluate expression (2.7).

In the following Sections we will apply the procedure outlined above to derive the relative velocity between the left and right states for the three different wave patterns (i)–(iii). We will denote with indices $`3`$ and $`3^{}`$ the quantities evaluated in the left ($`L_{}`$) and right ($`R_{}`$) states behind the discontinuity fronts, and where $`p_3=p_3^{}`$, $`v_3=v_3^{}`$, and $`\rho _3\rho _3^{}`$. Note also that in all of the different cases we will assume that $`p_1>p_2`$, and take the positive $`x`$-direction as being the one from the region $`1`$ to the region $`2`$ (An alternative opposite choice is also possible.).

## 4 $`L𝒮_{}L_{}𝒞R_{}𝒮_{}R`$: Two Shock Fronts

We start by considering the wave pattern produced by two shocks propagating in opposite directions (see Fig. 2). This situation is characterized by a value of the pressure downstream of the shocks which is larger than the pressures in the unperturbed states, i.e. $`p_3>p_1>p_2`$. Applying equation (3.8) to the shock front moving towards the left and evaluating it in the reference frame of the contact discontinuity, we can write the velocity ahead of the left propagating shock as

$$v_1=\sqrt{\frac{(p_3p_1)(e_3e_1)}{(e_1+p_3)(e_3+p_1)}}.$$
(4.13)
Similarly, we can apply equation (3.8) to the shock front moving towards the right and evaluate it in the frame comoving with the contact discontinuity to obtain that the velocity ahead of the right propagating shock is

$$v_2=\sqrt{\frac{(p_3p_2)(e_3^{}e_2)}{(e_2+p_3)(e_3^{}+p_2)}}$$
(4.14)
Equations (4.13) and (4.14) can now be used to derive the relativistic expression for the relative velocity of the flow ahead of the two shocks $`(v_{12})_{_{2S}}`$. As proved in Appendix A, the expression for the relative velocity between the unperturbed states is a monotonic function of $`p_3`$ for all possible wave patterns. In particular, for the present choice of initial data, this expression is a monotonically increasing function of $`p_3`$. As a result, the value of $`(v_{12})_{_{2S}}`$ can be used to build a criterion for the occurrence of two shocks propagating in opposite directions. In fact, since $`p_1`$ is the smallest value that $`p_3`$ can take, two shocks will form if

$$v_{12}>(\stackrel{~}{v}_{12})_{_{2S}}\sqrt{\frac{(p_1p_2)(\widehat{e}e_2)}{(\widehat{e}+p_2)(e_2+p_1)}},$$
(4.15)
where

$$\widehat{e}=\widehat{h}\widehat{\rho }p_1=\widehat{h}\frac{\gamma p_1}{(\gamma 1)(\widehat{h}1)}p_1,$$
(4.16)
and $`\widehat{h}`$ is the only positive root of the Taub adiabat (Taub 1978, Martí & Müller 1994)

$$\left[1+\frac{(\gamma 1)(p_2p_3)}{\gamma p_3}\right]\widehat{h}^2\frac{(\gamma 1)(p_2p_3)}{\gamma p_3}\widehat{h}+\frac{h_2(p_2p_3)}{\rho _2}h_2^2=0.$$
(4.17)
when $`p_3p_1`$. Simple calculations reported in Appendix B show that the Newtonian limit of $`(\stackrel{~}{v}_{12})_{_{2S}}`$ corresponds to the expression derived by Landau and Lifshitz (1987).

## 5 $`L𝒮_{}L_{}𝒞R_{}_{}R`$: One Shock and one Rarefaction Wave

We next consider the wave pattern produced by one shock front propagating towards the right and one rarefaction wave propagating in the opposite direction (see Fig. 3). This situation is therefore characterized by $`p_1>p_3>p_2`$.

Evaluating expression (3.12) in the reference frame comoving with the contact discontinuity, we can evaluate the flow velocity ahead of the rarefaction wave to be

$$v_1=\frac{1A_+(p_3)}{1+A_+(p_3)},$$
(5.18)
where

$$A_+(p_3)\left\{\left[\frac{(\gamma 1)^{1/2}c_s(p_3)}{(\gamma 1)^{1/2}+c_s(p_3)}\right]\left[\frac{(\gamma 1)^{1/2}+c_s(p_1)}{(\gamma 1)^{1/2}c_s(p_1)}\right]\right\}^{2/(\gamma 1)^{1/2}}.$$
(5.19)
The flow velocity ahead of the shock front can be derived as in Section 4 by evaluating equation (3.8) in the reference frame of the contact discontinuity to get

$$v_2=\sqrt{\frac{(p_3p_2)(e_3^{}e_2)}{(e_3^{}+p_2)(e_2+p_3)}},$$
(5.20)
which, combined with expression (5.18), can be used to derive the relativistic expression for the relative velocity of the fluids ahead of the shock and ahead of the rarefaction wave $`(v_{12})_{_{SR}}`$. As for $`(v_{12})_{_{2S}}`$, it can be shown that $`(v_{12})_{_{SR}}`$ is a monotonically increasing function of $`p_3`$ (see Fig. 1 and Appendix A for an analytic proof). Exploiting now the knowledge that for this wave pattern the pressure in the region between the two waves must satisfy $`p_2<p_3<p_1`$, we can establish that the criterion on the relative velocity for having one shock and one rarefaction wave is

$$(\stackrel{~}{v}_{12})_{_{SR}}=\frac{1A_+(p_3)}{1+A_+(p_3)}|_{p_3=p_2}<v_{12}\sqrt{\frac{(p_1p_2)(\widehat{e}e_2)}{(e_2+p_1)(\widehat{e}+p_2)}}=(\stackrel{~}{v}_{12})_{_{2S}}.$$
(5.21)
Note that the upper limit of (5.21) coincides with $`(\stackrel{~}{v}_{12})_{_{2S}}`$, which is the lower limit for the occurrence of 2 shock waves (4.15) and whose Newtonian limit coincides with the equivalent one found by Landau and Lifshitz (1987) (see Appendix B). Note also that in the limit, $`p_3p_2`$, regions 1 and 2 are connected by a single rarefaction wave. In this case the sound speed can be computed using $`p_3=p_2`$ but with $`\rho _3=\rho _1(p_2/p_1)^{1/\gamma }`$. Finally, note that this is the only wave pattern in which $`v_1`$ and $`v_2`$ have the same sign and it therefore includes the classical shock-tube problem, where $`v_1=v_2=0`$.

## 6 $`L_{}L_{}𝒞R_{}_{}R`$: Two Rarefaction Waves

We now consider the wave pattern produced by two rarefaction waves propagating in opposite directions (see Fig. 4). This situation is characterized by $`p_1>p_2>p_3`$ and when the waves are sufficiently strong it might lead to a vacuum region ($`\rho _3=0`$) behind the rarefaction waves.

Following again the same procedure discussed in the previous Section, we can determine the values of the fluid velocities ahead of the two rarefaction waves as, respectively,

$`v_1`$ $`=`$ $`{\displaystyle \frac{A_+(p_3)1}{A_+(p_3)+1}},`$ (6.22)
$`v_2`$ $`=`$ $`{\displaystyle \frac{1A_{}(p_3^{})}{1+A_{}(p_3^{})}},`$ (6.23)
where

$$A_{}(p_3^{})\left\{\left[\frac{(\gamma 1)^{1/2}c_s(p_3^{})}{(\gamma 1)^{1/2}+c_s(p_3^{})}\right]\left[\frac{(\gamma 1)^{1/2}+c_s(p_2)}{(\gamma 1)^{1/2}c_s(p_2)}\right]\right\}^{2/(\gamma 1)^{1/2}}.$$
(6.24)
We have indicated with $`c_s(p_3^{})`$ the sound speed in the region $`3^{}`$, which differs from the one in region $`3`$ because of the jump in the densities $`\rho _3`$ and $`\rho _3^{}`$. The relative velocity built using (6.22) and (6.23) is then

$$(v_{12})_{_{2R}}=\frac{A_+(p_3)A_{}(p_3^{})}{A_+(p_3)+A_{}(p_3^{})}.$$
(6.25)
As for the relative velocities of the previous wave patterns, it can be shown that $`(v_{12})_{_{2R}}`$, is a monotonically increasing function of $`p_3`$ (see Fig. 1 and Appendix A for an analytic proof) so that the criterion for the occurrence of two rarefaction waves can be expressed as

$$(v_{12})_{_{2R}}|_{p_3=0}<v_{12}(v_{12})_{_{2R}}|_{p_3=p_2}=(\stackrel{~}{v}_{12})_{_{SR}},$$
(6.26)
where $`A_{}(p_3=p_2)=1`$ and therefore the upper limit for (6.26) coincides with the lower limit for (5.21).

The condition (6.26) can also be expressed in a more useful form as

$$\frac{S_1S_2}{S_1+S_2}<v_{12}\frac{1A_+(p_2)}{1+A_+(p_2)},$$
(6.27)
where, the constants $`S_1`$ and $`S_2`$ are shorthand for

$`S_1`$ $``$ $`\left[{\displaystyle \frac{(\gamma 1)^{1/2}+c_s(p_1)}{(\gamma 1)^{1/2}c_s(p_1)}}\right]^{2/(\gamma 1)^{1/2}}`$ (6.28)
$`S_2`$ $``$ $`\left[{\displaystyle \frac{(\gamma 1)^{1/2}+c_s(p_2)}{(\gamma 1)^{1/2}c_s(p_2)}}\right]^{2/(\gamma 1)^{1/2}}.`$ (6.29)
An important property of equation (6.25) is that it can be inverted analytically. This involves rewriting it in terms of a quartic equation in the unknown sound speed in the region $`L_{}`$ (see Appendix C for the explicit form of the equation). Once the relevant real root of this equation has been calculated analytically, the value of the pressure $`p_{}`$ can be found through a simple algebraic expression. In this case, therefore, the solution of the exact relativistic Riemann problem can be found in an analytic closed form.

We conclude the analysis of the different wave patterns with a comment on the case of two rarefaction waves propagating in opposite directions and leaving behind them a region with zero density and pressure. This situation occurs when the fluids in the regions 1 and 2 are moving sufficiently fast in opposite directions. This is the case whenever the relative velocity between the two states ahead of the rarefaction waves is less than or equal to the lower limit for $`(v_{12})_{_{2R}}`$, i.e.

$$v_{12}(\stackrel{~}{v}_{12})_{_{2R}}\frac{S_1S_2}{S_1+S_2}.$$
(6.30)
Also in this case, taking the Newtonian limit of $`(\stackrel{~}{v}_{12})_{_{2R}}`$ we obtain the corresponding expression derived by Landau and Lifshitz (1987) (see Appendix B). Finally, note that when a vacuum is produced, $`v_{12}`$ is no longer dependent on $`p_3`$ and this branch of the curve cannot be plotted in Fig.1.

## 7 Numerical Implementation

The core of most exact Riemann solvers, in both Newtonian and relativistic hydrodynamics, is based on the numerical computation of the pressure in the regions $`L_{}`$ and $`R_{}`$ that form behind the waves. The key property exploited when performing the numerical calculation is that the velocity in such regions can be expressed as a monotonic function of the pressure i.e. $`v_L_{}=v_L_{}(p_L_{})`$, and $`v_R_{}=v_R_{}(p_R_{})`$. Since there is no jump across a contact discontinuity in either the velocity or in the pressure, the numerical solution of the Riemann problem consists then of finding the root $`p_{}=p_L_{}=p_R_{}`$ of the nonlinear equation

$$v_L_{}(p_{})v_R_{}(p_{})=0,$$
(7.31)
where $`v_L_{},v_R_{}`$ have different functional forms according to the different wave patterns produced. This method has two obvious disadvantages: (1) it cannot determine, using the initial conditions, the wave pattern produced and thus which of the functional forms to use for $`v_L_{},v_R_{}`$; (2) it cannot provide a straightforward bracketing interval for the root. In practice, however, these difficulties are effectively balanced by efficient algorithms based on a sequence of trial and error attempts that rapidly bracket the root and determine the correct equation to solve (see, for instance, Martí, & Müller 1999 and the algorithm presented therein).

The exact Riemann solver which we propose here differs from the one discussed above mostly because it avoids the disadvantages (1) and (2). In fact, as discussed in Section 1, by comparing the relative velocity between the initial left and right states $`(v_{12})_0`$ with the relevant limiting values constructed from the initial conditions $`(\stackrel{~}{v}_{12})_{_{2S}}`$, $`(\stackrel{~}{v}_{12})_{_{SR}}`$, $`(\stackrel{~}{v}_{12})_{_{2R}}`$, we can determine both the wave pattern which will be produced and the correct bracketing range in the pressure. Once this information has been obtained, the Riemann problem can be solved either through the solution of equation (7.31) or, equivalently, by looking for the value of the pressure $`p_{}`$ which would produce a relative velocity $`(v_{12})_0`$. This latter approach, which we will be discussing in the following, involves then the solution of the nonlinear equation

$$v_{12}(p_{})(v_{12})_0=0,$$
(7.32)
where $`v_{12}(p_{})`$ is given by the expressions for $`(v_{12})_{_{2S}}`$, or $`(v_{12})_{_{SR}}`$, or $`(v_{12})_{_{2R}}`$ derived in Sections 46. Furthermore, in the case of two rarefaction waves, equation (7.32) can also be solved analytically.

Besides providing direct information about the wave pattern produced, about the correct equation to solve and the relevant bracketing interval, our approach is also very simple to implement numerically. In practice, the basic steps for the solution of the Riemann problem can be summarised as follows:

1. Evaluate from the initial conditions the three limiting relative velocities $`(\stackrel{~}{v}_{12})_{_{2S}}`$, $`(\stackrel{~}{v}_{12})_{_{SR}}`$, $`(\stackrel{~}{v}_{12})_{_{2R}}`$.
2. Determine the wave pattern and the functional form of $`v_{12}(p_{})`$ by comparing $`(v_{12})_0`$ with the limiting values calculated in (a) and according to the scheme below.

| $`(v_{12})_0`$ | $`>`$ | $`(\stackrel{~}{v}_{12})_{_{2S}}`$: | $`𝒮_{}𝒞𝒮_{},`$ | $`v_{12}(p_{})=(v_{12})_{_{2S}}`$ |
| --- | --- | --- | --- | --- |
| $`(\stackrel{~}{v}_{12})_{_{SR}}<(v_{12})_0`$ | $``$ | $`(\stackrel{~}{v}_{12})_{_{2S}}`$: | $`L_{}L_{}𝒞R_{}𝒮_{}R,`$ | $`v_{12}(p_{})=(v_{12})_{_{SR}}`$ |
| $`(\stackrel{~}{v}_{12})_{_{2R}}<(v_{12})_0`$ | $``$ | $`(\stackrel{~}{v}_{12})_{_{SR}}`$: | $`L_{}L_{}𝒞R_{}_{}R,`$ | $`v_{12}(p_{})=(v_{12})_{_{2R}}`$ |
| $`(v_{12})_0`$ | $``$ | $`(\stackrel{~}{v}_{12})_{_{2R}}`$: | $`L_{}L_{}𝒞R_{}_{}R`$ with vacuum, | $``$ |
3. According to the wave pattern found, determine the extremes $`p_{max}`$ and $`p_{min}`$ of the pressure interval bracketing $`p_{}`$. Within our conventions this is equivalent to setting¶¶¶¶¶¶$`\mathrm{}`$Note that in practice, the upper limit for the pressure in the case of two shocks is found by starting from a reasonable value above $`p_{min}`$, which is incremented until the solution is effectively bracketed.

|  | $`𝒮_{}𝒞𝒮_{}`$ | $`L_{}L_{}𝒞R_{}𝒮_{}R`$ | $`L_{}L_{}𝒞R_{}_{}R`$ |
| --- | --- | --- | --- |
| $`p_{min}`$ | max$`(p_1,p_2)`$ | min$`(p_1,p_2)`$ | 0 |
| $`p_{max}`$ | $`\mathrm{}`$ | max$`(p_1,p_2)`$ | min$`(p_1,p_2)`$ |
4. Solve equation (7.32) and determine $`p_{}`$.
5. Complete the solution of the Riemann problem by computing the remaining variables of the intermediate states $`L_{}`$ and $`R_{}`$.

We have implemented our algorithm for an exact Riemann solver and have tested it for a range of Riemann problems. We have also compared the performance of our algorithm with the “standard” approach presented by Martí & Müller 1999 and have found a systematic reduction in the computational costs for the same level of accuracy in the solution. The quantitative efficiency improvement depends on the type of problem under consideration. In the case of a generic hydrodynamical problem (in which very simple Riemann problems are solved), our approach brackets the solution very closely and this produces substantial computational improvements of up to $`30\%`$. However, for the specific cases discussed in this paper, where very strong shocks have been considered, the speed-up is smaller and of the order of $`10\%`$. It is worth noting that such an improvement could reduce appreciably the computational costs in three-dimensional relativistic hydrodynamics codes, where this operation is repeated a very large number of times.

## 8 Conclusion

We have presented a new procedure for the numerical solution of the exact Riemann problem in relativistic hydrodynamics. In this approach special attention is paid to the relativistic invariant expression for the relative velocity $`v_{12}`$ between the unperturbed left and right states. This has been shown to be a monotonic function of the pressure $`p_{}`$ in the region formed between the wave fronts. The determination of this pressure is the basic step in the solution of the exact Riemann problem.

The use of the relative velocity has a number of advantages over alternative exact Riemann solvers discussed in the literature. In particular, it extracts the information implicitly contained in the data for two initial states to deduce the wave pattern that will be produced by the decay of the discontinuity between these two states. This, in turn, allows an “a priori” determination to be made of the interval in pressure bracketing $`p_{}`$ and the correct functional form for the nonlinear equation whose root will solve the exact Riemann problem. All of these advantages translate, in practice, into a simpler algorithm to implement and an improved efficiency in the numerical solution of the Riemann problem. Furthermore, in the case of two rarefaction waves propagating in opposite directions (and in strict analogy with what happens in Newtonian hydrodynamics), the use of the relative velocity allows for the solution of the Riemann problem in an analytic closed form.

Because of all of the advantages discussed above, its intrinsic simplicity of implementation and the numerical efficiency gain it produces, this new exact Riemann solver should be considered as an interesting alternative to the traditional exact Riemann solver presently discussed in the literature. Investigations about the extension of this approach to multidimensions are currently in progress.

###### Acknowledgements.

It is a pleasure to thank J. A. Font, C. Gheller, and J. C. Miller for useful discussions. Financial support for this research has been provided by the Italian Ministero dell’Università e della Ricerca Scientifica and by the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).”

## Appendix A: Monotonicity of the relative velocity as function of $`p_{}`$

We here prove that $`v_{12}`$ is a monotonic function of $`p_{}`$ for all of the possible wave patterns. In particular, because of our choice of referring the initial left state as to “1” and the initial right state as to “2”, we will here show that $`v_{12}`$ is a monotonically increasing function of $`p_{}`$.

Denoting by $`a^{}`$ the first derivative of the quantity $`a`$ with respect to $`p_{}=p_3=p_3^{}`$, it is straightforward to obtain that the first derivative of (2.7) is

$$v_{12}^{}=\frac{v_1^{}(1v_2^2)v_2^{}(1v_1^2)}{(1v_1v_2)^2}.$$
(8.33)
Since $`v_11,v_21`$, the two terms in the parentheses of (8.33) are always positive, and the proof that $`v_{12}`$ is monotonically increasing will follow if it can be shown that $`v_1^{}`$ and $`v_2^{}`$ are both positive. We will do so for each of the three possible wave patterns.

### Two Shock Fronts

Taking the derivative of the fluid velocity ahead of the left propagating shock front (measured from the contact discontinuity) \[cf. equation (4.13)\] we obtain

$$v_1^{}=\frac{(e_1+p_1)\left[(e_3e_1)(e_3+p_1)+(p_3p_1)(e_1+p_3)e_3^{}\right]}{2v_1(e_1+p_3)^2(e_3+p_1)^2}.$$
(8.34)
Since the energy density is an increasing function of pressure, $`e_3^{}>0`$; furthermore, $`p_3p_1`$ and $`e_3e_1`$ for the wave pattern considered and $`v_1^{}>0`$ as a result. The equivalent expression for the derivative of the fluid velocity ahead of the right propagating shock front (measured from the contact discontinuity) \[cf. eq (4.14)\] can be obtained by replacing in equation (8.34) the indices $`1`$ and $`3`$ by $`2`$ and $`3^{}`$ respectively, i.e.

$$v_2^{}=\frac{(e_2+p_2)\left[(e_3^{}e_2)(e_3^{}+p_2)+(p_3^{}p_2)(e_2+p_3^{})e_3^{}^{}\right]}{2v_2(e_2+p_3^{})^2(e_3^{}+p_2)^2}.$$
(8.35)
Since now for the wave pattern considered: $`v_2<0`$, $`p_3^{}p_2`$ and $`e_3^{}e_2`$, we are led to conclude that $`v_2^{}>0`$, thus making the overall $`v_{12}^{}`$ positive for any value of $`p_3`$.

### One Shock and one Rarefaction Wave

In this case we only need to show that $`v_1^{}>0`$ since for the velocity ahead of the right propagating shock front we can use the results derived in (8.35). Taking the derivative of expression (5.18) then yields

$$v_1^{}=\frac{2A_+^{}(p_3)}{[1+A_+(p_3)]^2},$$
(8.36)
where

$`A_+^{}(p_3)`$ $`=`$ $`4\left|{\displaystyle \frac{(\gamma 1)^{1/2}+c_s(p_1)}{(\gamma 1)^{1/2}c_s(p_1)}}\right|{\displaystyle \frac{A_+(p_3)^{(2\sqrt{\gamma 1})/2}}{[(\gamma 1)^{1/2}+c_s(p_3)]^2}}c_s^{}(p_3)`$ (8.37)
$``$ $`C_1c_s^{}(p_3),`$
where $`C_1>0`$. When the sound speed $`c_s(p_3)`$ is an increasing function of pressure, as is the case for a polytropic equation of state \[cf. eq (2.4)\], $`v_1^{}>0`$ and therefore $`v_{12}^{}>0`$.

### Two Rarefaction Waves

What we need to show in this case is that $`v_2^{}<0`$ since we can exploit the previous result that $`v_1^{}>0`$ where $`v_1`$ is the fluid velocity ahead of the left propagating rarefaction wave as measured from the contact discontinuity. In this case, taking the derivative of expression (6.23) one obtains

$$v_2^{}=\frac{2A_{}^{}(p_3^{})}{[1+A_{}(p_3^{})]^2},$$
(8.38)
where now

$`A_{}^{}(p_3^{})`$ $`=`$ $`4\left|{\displaystyle \frac{(\gamma 1)^{1/2}+c_s(p_2)}{(\gamma 1)^{1/2}c_s(p_2)}}\right|{\displaystyle \frac{A_{}(p_3^{})^{(2+\sqrt{\gamma 1})/2}}{[(\gamma 1)^{1/2}+c_s(p_3^{})]^2}}c_s^{}(p_3^{})`$ (8.39)
$``$ $`C_2c_s^{}(p_3^{})`$
with $`C_2>0`$ and therefore $`v_2^{}<0`$.

We have therefore shown that for all of the wave patterns considered $`v_1^{}>0`$ and $`v_2^{}>0`$, thus proving that $`v_{12}`$ is always a monotonically increasing function of $`p_{}`$.

## Appendix B: Newtonian limits of $`(\stackrel{~}{v}_{12})_{_{2S}}`$, $`(\stackrel{~}{v}_{12})_{_{SR}}`$, $`(\stackrel{~}{v}_{12})_{_{2R}}`$

We here show that the three limiting values of $`(\stackrel{~}{v}_{12})_{_{2S}}`$, $`(\stackrel{~}{v}_{12})_{_{SR}}`$ and $`(\stackrel{~}{v}_{12})_{_{2R}}`$ reduce to their Newtonian counterparts in the limit of $`v,c_s0`$, and $`h1`$. In particular, we will restrict ourselves to considering the case of a polytropic equation of state (2.4).

We start by considering the Newtonian limit of $`(\stackrel{~}{v}_{12})_{_{2S}}`$ which is obtained when $`p/e1`$ and $`e1/V`$, with $`V=1/\rho `$ being the specific volume. In this case, then

$`(\stackrel{~}{v}_{12})_{_{2S}}|_{\mathrm{Newt}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{(p_1p_2)(\widehat{e}e_2)}{\widehat{e}e_2}}}`$ (8.40)
$`=`$ $`\sqrt{(p_1p_2)(1/e_21/\widehat{e})}.`$
which coincides with the corresponding expression derived by Landau and Lifshitz (1987) but with inverted indices.

We next consider the Newtonian limit of $`(\stackrel{~}{v}_{12})_{_{SR}}`$ which is obtained when both $`c_s(p_3)1`$ and $`c_s(p_1)1`$. In this case,

$`A_+(p_3)`$ $``$ $`\left(1{\displaystyle \frac{c_s(p_3)}{\sqrt{\gamma 1}}}\right)^{4/\sqrt{\gamma 1}}\left(1+{\displaystyle \frac{c_s(p_1)}{\sqrt{\gamma 1}}}\right)^{4/\sqrt{\gamma 1}}`$ (8.41)
$``$ $`\left(1{\displaystyle \frac{4c_s(p_3)}{\gamma 1}}\right)\left(1+{\displaystyle \frac{4c_s(p_1)}{\gamma 1}}\right)`$
$``$ $`1{\displaystyle \frac{4}{\gamma 1}}(c_s(p_3)c_s(p_1)).`$
so that the Newtonian limit of $`(\stackrel{~}{v}_{12})_{_{SR}}`$ is given by

$`(\stackrel{~}{v}_{12})_{_{SR}}|_{\mathrm{Newt}}`$ $`=`$ $`{\displaystyle \frac{1A_+(p_3)}{2}}|_{p_3=p_2}`$ (8.42)
$``$ $`{\displaystyle \frac{2}{\gamma 1}}c_s(p_1)\left[1{\displaystyle \frac{c_s(p_3)}{c_s(p_1)}}\right]|_{p_3=p_2}.`$
Bearing in mind that

$$\frac{c_s(p_3)}{c_s(p_1)}|_{p_3=p_2}=\left(\frac{p_2}{p_1}\right)^{(\gamma 1)/2\gamma },$$
(8.43)
we obtain

$$(\stackrel{~}{v}_{12})_{_{SR}}|_{\mathrm{Newt}}=\frac{2}{\gamma 1}c_s(p_1)\left[1\left(\frac{p_2}{p_1}\right)^{(\gamma 1)/2\gamma }\right],$$
(8.44)
which again coincides with the corresponding expression derived by Landau and Lifshitz (1987) but with inverted indices.

Finally, we consider the Newtonian limit of $`(\stackrel{~}{v}_{12})_{_{2R}}`$ for $`c_s(p_1),c_s(p_2)1`$. In this case

$`S_1`$ $``$ $`1+{\displaystyle \frac{4c_s(p_1)}{\gamma 1}},`$
$`S_2`$ $``$ $`1{\displaystyle \frac{4c_s(p_2)}{\gamma 1}}.`$ (8.45)
so that the Newtonian limit is given by

$`(\stackrel{~}{v}_{12})_{_{2R}}|_{\mathrm{Newt}}={\displaystyle \frac{S_2S_1}{S_1+S_2}}`$ $`=`$ $`\left[{\displaystyle \frac{4c_s(p_1)}{\gamma 1}}{\displaystyle \frac{4c_s(p_2)}{\gamma 1}}\right]\left[2+{\displaystyle \frac{4}{\gamma 1}}[c_s(p_1)c_s(p_2)]\right]^1`$ (8.46)
$`=`$ $`\left[{\displaystyle \frac{2c_s(p_1)}{\gamma 1}}{\displaystyle \frac{2c_s(p_2)}{\gamma 1}}\right]\left[1{\displaystyle \frac{2}{\gamma 1}}[c_s(p_1)c_s(p_2)]\right]`$
$`=`$ $`{\displaystyle \frac{2c_s(p_1)}{\gamma 1}}{\displaystyle \frac{2c_s(p_2)}{\gamma 1}}.`$
Once more, expression (8.46) coincides with the corresponding expression derived by Landau and Lifshitz (1987) but with inverted indices.

## Appendix C: A closed form solution in the case of Two Rarefaction Waves

As discussed in Sections 6 and 7, when $`(\stackrel{~}{v}_{12})_{_{2R}}<(v_{12})_0<(\stackrel{~}{v}_{12})_{_{SR}}`$, the initial conditions give rise to two rarefaction waves and it is possible to derive a closed form solution for the unknown pressure $`p_{}`$. In this way we can, at least in principle, avoid any numerical root finding procedure and determine the solution exactly. In this Appendix we first derive this analytic solution in the context of relativistic hydrodynamics and then calculate its Newtonian limit. We will restrict ourselves to considering the particular case of a polytropic equation of state (2.4).

Using expression (3.11) we can write the pressures $`p_3`$ and $`p_3^{}`$ as functions of the sound speeds $`c_s(p_3)`$ and $`c_s(p_3^{})`$ which, for convenience, we will hereafter refer to as $`x`$ and $`x^{}`$ respectively

$`p_3`$ $`=`$ $`k_1^{1/(\gamma 1)}\left[{\displaystyle \frac{x^2(\gamma 1)}{\gamma (\gamma 1)\gamma x^2}}\right]^{\gamma /(\gamma 1)},`$ (8.47)
$`p_3^{}`$ $`=`$ $`k_2^{1/(\gamma 1)}\left[{\displaystyle \frac{(x^{})^2(\gamma 1)}{\gamma (\gamma 1)\gamma (x^{})^2}}\right]^{\gamma /(\gamma 1)},`$ (8.48)
where $`k_1=p_1/\rho _1^\gamma `$ and $`k_2=p_2/\rho _2^\gamma `$ are the two constants entering the polytropic equation. Since $`p_3=p_3^{}=p_{}`$, we can obtain the following relation between $`x^{}`$ and $`x`$:

$$(x^{})^2\frac{\gamma (\gamma 1)x^2}{\gamma x^2(1\alpha )+\alpha \gamma (\gamma 1)},$$
(8.49)
where $`\alpha (k_1/k_2)^{1/\gamma }`$. The expression for the relative velocity (6.25) can also be written as

$$\frac{A_+(p_3)}{A_{}(p_3^{})}=\frac{1(v_{12})_0}{1+(v_{12})_0}.$$
(8.50)
and we then use expressions (5.19) and (6.24) to expand the left hand side of (8.50). After some algebra we are then left with

$$\left(\frac{\mathrm{\Gamma }x}{\mathrm{\Gamma }+x}\right)\left(\frac{\mathrm{\Gamma }x^{}}{\mathrm{\Gamma }+x^{}}\right)=\left[\frac{\mathrm{\Gamma }c_s(p_2)}{\mathrm{\Gamma }+c_s(p_2)}\right]\left[\frac{\mathrm{\Gamma }c_s(p_1)}{\mathrm{\Gamma }+c_s(p_1)}\right]\left[\frac{1(v_{12})_0}{1+(v_{12})_0}\right]^{\mathrm{\Gamma }/2},$$
(8.51)
where $`\mathrm{\Gamma }^2\gamma 1`$ and the right hand side of (8.51) is a constant which we rename as

$$\mathrm{\Pi }\left[\frac{\mathrm{\Gamma }c_s(p_2)}{\mathrm{\Gamma }+c_s(p_2)}\right]\left[\frac{\mathrm{\Gamma }c_s(p_1)}{\mathrm{\Gamma }+c_s(p_1)}\right]\left[\frac{1(v_{12})_0}{1+(v_{12})_0}\right]^{\mathrm{\Gamma }/2}.$$
(8.52)
Introducing now the auxiliary quantity

$$\beta \frac{1+\mathrm{\Pi }}{1\mathrm{\Pi }},$$
(8.53)
expression (8.51) can be written as

$$(x^{})^2=\left[\frac{\mathrm{\Gamma }(x\beta \mathrm{\Gamma })}{x\mathrm{\Gamma }\beta }\right]^2.$$
(8.54)
Comparing now (8.49) and (8.54) gives a $`4`$-th order equation in the unknown sound velocity $`x`$

$$a_0x^4+a_1x^3+a_2x^2+a_3x+a_4=0,$$
(8.55)
where

$`a_0`$ $``$ $`1\beta ^2(1\alpha ),`$ (8.56)
$`a_1`$ $``$ $`2\mathrm{\Gamma }\alpha \beta ,`$ (8.57)
$`a_2`$ $``$ $`\mathrm{\Gamma }^2(1\alpha )(\beta ^21),`$ (8.58)
$`a_3`$ $``$ $`2\mathrm{\Gamma }^3\alpha \beta ,`$ (8.59)
$`a_4`$ $``$ $`\alpha \mathrm{\Gamma }^4.`$ (8.60)
The analytic solution of equation (8.55) will yield at least two real roots, one of which will be the physically acceptable one: i.e. positive, less than one, and such that the pressure $`p_{}`$ falls in the relevant bracketing interval.

In its Newtonian limit, equation (8.55) reduces to a second order equation in the unknown sound velocity

$$\left(\frac{1}{\alpha }1\right)x^2+2\mathrm{\Sigma }x\mathrm{\Sigma }^2=0,$$
(8.61)
where

$$\mathrm{\Sigma }c_s(p_1)+c_s(p_2)+\frac{(\gamma 1)}{2}v_{12},$$
(8.62)
and $`v_{12}=v_1v_2`$. The fact that the Newtonian Riemann problem in the case of two rarefaction waves can be solved analytically is well known and is at the basis of the so called “Two Rarefaction Approximate Riemann Solver” (Toro 1997).


format: text

text:
# Time Scales in Long GRBs

## 1 Introduction

The light curve of a long duration Gamma-Ray Burst (GRB) (Kouveliotou 1993) is usually complex, it is composed from several dozens of short (about 1sec) distinguishable pulses. The variability of these bursts has been already used (Sari & Piran, 1997; Fenimore, Ramirez & Sumner 1997) to suggest that internal shocks , rather than external shocks produce most GRBs (simple models of the latter cannot produce efficiently variable bursts). The recent discovery of correlation relating properties of the temporal structure with the burst luminosity offer possibility of deriving independent estimates of the redshift of GRBs (Stern, Poutanen & Svensson 1997, 1999; Norris , Marani & Bonnell 2000; Fenimore & Ramirez-Ruiz 2000; Reichart et al 2000).

The temporal properties of long duration GRBs , and specifically the pulses and intervals properties, were investigated in previous works (e.g. McBreen et.al. 1994, Li & Fenimore 1996, Norris et. al. 1996, Quilligan et. al. 1999). We present here a further study of the temporal properties of GRB light curves. Using a new algorithm (Nakar & Piran, 2001b), which is based on Li & Fenimore (1996), we analyze the distribution of the time intervals between pulses and the pulses width in long bright bursts. We find that the pulses width, $`\delta t`$, distribution is consistent with a lognormal distribution. The distribution of intervals between pulses, $`\mathrm{\Delta }t`$, becomes lognormal, with a comparable width to the $`\delta t`$ distribution, provided that we eliminate intervals that contain quiescent times. Namely, long periods with no evidence of photon counts above the background, between periods of strong $`\gamma `$-ray emission (Ramirez & Merloni 2001, Nakar & Piran 2001a). In nature such a deviation from a lognormal distribution occurs when a different mechanism govern the high end tail of the distribution. The similarity between the pulses width and the intervals distributions suggest that both distributions are influenced by the same (source) physical parameters. In this case, we expect also a correlation between the pulse width and its adjacent intervals. We search for this correlation using 12 bursts which contain enough pulses for such analysis to be applicable. In all these bursts there is a significant correlation between the pulse width and the duration of the interval preceding it. A significant correlation between the pulse width and the following interval was found in 60% of the bursts. In all the bursts the significance of the correlation between the pulse width and the preceding interval was higher or equal to this of the following interval.

## 2 The Algorithm and the Data

We have applied a modified Li & Fenimore (1996) peak finding algorithm (Nakar & Piran 2001b) to a sample of 68 long bursts ($`T_{90}>2sec`$, where $`T_{90}`$ is the time required to accumulate from 5 to 95 per cent of the total counts). This algorithm searches for individual peaks in BATSE 64ms concatenate data. In the analysis we consider the sum of all four energy channels, i.e. E $`>`$ 25Kev. This data type includes the photon counts, in a 64ms time bins, from a few seconds before the bursts trigger until a few hundred seconds after the trigger. The minimal $`\mathrm{\Delta }t`$ and $`\delta t`$ are 0.128sec in this resolution. Each peak the algorithm find corresponds to one pulse. We determine the pulse width and the interval between two maxima, which we define as the interval between pulses (Norris et. al. 1996, Li & Fenimore 1996). The algorithm also searches for quiescent times. The definition of minimal duration of a quiescent time is somewhat arbitrary (whether a single time bin with a count level of the background is a quiescent time or not). We demand that the average of ten bins would be at the background level to be considered as a quiescent time. Hence the minimal quiescent time in our analysis is about 1sec.

The sample we consider is the 68 brightest long bursts in the BATSE 4B catalog (peak flux in 64ms$`>`$10.19$`ph/(seccm^2)`$). This resulted in 1330 pulses and 1262 intervals. The analysis below is based on the width of these pulses and the corresponding intervals between them.

One may wonder if the brightest bursts are representative. Figure 1 depicts the histogram of $`T_{90}`$ of our sample compared to the histogram of all the long bursts in 4B catalog. The histograms are similar. Therefore our sample, at least as far as duration is concerned, is representative for the entire 4B catalog.

We also analyze 24 dimmer bursts, which are still bright enough for the analysis to be carried out. These bursts are chosen randomly from the bursts with peak flux at 64ms greater than 4$`ph/(seccm^2)`$ but smaller than 10$`ph/(seccm^2)`$. This sample included 438 intervals. The results for this sample are consistent with those obtained for the brightest sample. Again indicating that the brightest sample is representative.

## 3 Temporal properties

### 3.1 Pulse Width Distribution

The distribution of the 1330 pulses’ widths is in excellent agreement with a lognormal distribution. The $`\chi ^2`$ test gives a probability of $`0.52`$ that the data was taken from a lognormal distribution, with $`\mu =0.065\pm 0.04`$ ($`\delta t_{avg}1sec`$) and $`\sigma =0.77\pm 0.03`$ (corresponding at $`1\sigma `$ to 0.5- 2.3sec). Figure 2 depicts the histogram and the cumulative distribution of the pulses’ widths in the bright sample. The analysis of the dimmer sample (see sec 2) show also an excellent agreement with a lognormal distribution. The probability that the dimmer sample’s pulses width are taken from a lognormal distribution is 0.9.

### 3.2 Intervals Distribution and Quiescent Times

McBreen (1994) and Li & Fenimore (1996) suggested that the distribution of the intervals between pulses, $`\mathrm{\Delta }t`$, is lognormal. We have found that the $`\delta t`$ distribution is lognormal. These results suggest that we consider first the null hypothesis that $`\mathrm{\Delta }t`$ distribution is also lognormal.

Figure 3a depicts the histogram of the time intervals between neighboring peaks, $`\mathrm{\Delta }t`$. Figure 3b shows the cumulative distribution of $`log_{10}(\mathrm{\Delta }t)`$, compared to a best-fit Gaussian. Both figures show a clear excess of many long intervals. Using the $`\chi ^2`$ test, we find a probability of $`1.210^{10}`$ that the data was taken from a lognormal distribution. The null hypothesis clearly fails.

Li & Fenimore (1996) already noticed such a deviation. A hints to this deviation could also be seen in the results presented by Norris et. al. (1996). Li & Fenimore (1996) and McBreen (personal communication 2000) have suggested that this deviation arises due to the limited resolution of the observations, and therefore the intervals distribution is consistent with a lognormal one. To test this hypothesis we examine in Figure 4 the cumulative probability of a mirror image of the right half of the $`\mathrm{\Delta }t`$ histogram. This half is insensitive to the limited resolution. Again the figure show a deviation from a lognormal distribution and an excess of long intervals (and short ones which are of course the mirror of the long intervals). This indicates that $`\mathrm{\Delta }t`$ distribution is not a lognormal.

The analysis of the dimmer bursts sample (see sec. 2) gives similar results. Figure 5 shows the cumulative distribution of $`log_{10}(\mathrm{\Delta }t)`$ in the dimmer sample compared to the best-fit Gaussian. The excess of long intervals is clear. The $`\chi ^2`$ test gives a probability of $`910^6`$ that this data was taken from a lognormal distribution.

We found that the long intervals between neighboring peaks are often dominated by quiescent times. Such periods are seen in 35 of the bursts in our sample (all together there are 50 quiescent times). Most of the bursts contain one or two quiescent times, some contain three. The quiescent times last typically several tens of seconds and they range between a second (our arbitrary lower limit) to hundreds of seconds. In some bursts the quiescent times are a significant fraction of the total duration.

Since the quiescent times contribute to the duration of the longest intervals we examined the possibility that they cause the excess of long intervals (over a lognormal distribution). We examine the interval distribution excluding now all intervals that contained a quiescent time. Figure 6 shows the histogram and the cumulative distribution of $`log_{10}(\mathrm{\Delta }t)`$ excluding the intervals that contained a quiescent time, compared to the best-fit Gaussian with $`\mu =0.257\pm 0.051`$ ($`\mathrm{\Delta }t_{avg}1.3sec`$) and $`\sigma =0.90\pm 0.04`$ (corresponding at $`1\sigma `$ to 0.53-3.1sec). The fit is good. The $`\chi ^2`$ test gives a probability of $`0.27`$ that this data was taken from a lognormal distribution. Moreover, the modified $`\mathrm{\Delta }t`$ distribution has a comparable width to the $`\delta t`$ distribution (which is unaffected by quiescent times).

### 3.3 The Correlation Between Pulses and Neighboring Intervals

The similarity between the two distributions motivated us to explore the correlation between pulses and their neighboring intervals. We calculated Pearson’s linear correlation coefficient, $`r`$, and Spearman’s ranked order coefficient, $`r_s`$, between a pulse duration and the intervals just preceding it, and just following it (excluding intervals that contain a quiescent time). In order to avoid any influence of the algorithm on the correlation we considered here only well separated pulses, namely, pulses whose half maximum (used to determine the width) is above the minimum between the pulse and its neighbors. In this case there is no worry about overlap between the pulses that may introduce a spurious correlation.

We considered only bursts with more than 12 separated pulses, all together 12 bursts of our sample with an average of 23 pulses per burst111The comparison was done burst by burst in order to eliminate redshift or intrinsic effects that would have produced spurious correlations if we had considered the whole data as one set . Our null hypothesis is of no correlation between the pulses width and the intervals. The probabilities of rejection of the null hypothesis found to be similar for both methods, $`r`$ and $`r_s`$. All the considered bursts showed a positive correlation between a pulse width and the preceding interval. The rejection of the null hypothesis is above 90% in half of the bursts and above 70% in the rest. The correlation between a pulse width and the following interval is less significant. In 5 of the bursts the null hypothesis is rejected with a confidence of more then 90%, and In 7 of the bursts the rejection confidence level is above 70%. In 4 bursts there is no evidence for correlation. In all the considered bursts the significance of the correlation between the pulse width and the preceding interval is higher or equal to the correlation with the following interval.

Figure 7 shows the pulses width and the preceding interval duration of the most correlated burst and the least correlated one. In both bursts a correlation is clearly seen. The best linear fit parameters of each burst are different, but the variance of the slope parameter is not large, compared with the typical slope. The average linear fit is $`log(\mathrm{\Delta }t)=(0.49\pm 0.19)log(\delta t)+(0.25\pm 0.11)`$. Figure 8 shows the pulses width and the following interval duration of the most correlated burst and the least correlated burst. While the correlation is clearly seen in trigger 2156, there is no evidence for correlation in trigger 3408. The average Linear fit for the correlated bursts (significance level above 70%) is $`log(\mathrm{\Delta }t)=(0.55\pm 0.11)log(\delta t)+(0.26\pm 0.11)`$.

The correlation disappears if we consider non neighboring intervals and pulses. For example, only half of the considered bursts show a positive correlation between a pulse duration and the interval after the following pulse. In no burst the rejection of the null hypothesis of no correlation is above 90%.

## 4 Discussion

The pulse duration distribution is consistent with a lognormal one. However, the distribution of the intervals between pulses is inconsistent with a lognormal distribution. Removal of the intervals that include quiescent times results in a distribution consistent with lognormal. This suggests that the distribution of interval between pulses is made of the sum of two different contributions. A lognormal distribution that is similar to the pulse width distribution, and the quiescent times distribution. The last one is dominant at the high end tail of the intervals distribution and produce the deviation from a lognormal distribution. We stress (see also Ramirez & Merloni 2001) that the quiescent times that we find here are not the same as the gaps between the precursors and the main pulses (Koshut et al, 1995) found in 3% of the bursts. In that case the gaps are between a softer and weaker component (the precursor) and a harder and much stronger one (the main burst). Here we observe gaps between pulses with the same characteristics. Moreover, in many cases we observe several quiescent times within the same burst.

The existence of quiescent times, and of two different mechanisms that control the gaps between pulses pose a new challenge to any model attempting to explain the temporal structure of GRBs. This suggests that there are three time scales of different nature within GRBs : (i) The shortest time scale, the variability scale which determines the pulses duration and the interval between pulses (both have a similar time-scale). (ii) An intermediate time scale producing long periods within the bursts with no activity (quiescent times). (iii) The longest time scale which corresponds to the duration of the burst.

Within the popular internal shock model the observed temporal structure reflects the activity of the “inner engine” (e.g. Kobayashi, Piran & Sari 1997; Ramirez & Fenimore, 2000). Within this model quiescent times reflect, most likely, periods in which the “inner engine” is not active at all. Alternatively they could reflect periods in which the “inner engine” emits a sequence of shells that do not collide (e.g. shells with a decreasing Lorentz factors) (Ramirez, Merloni & Rees 2000). In both cases the quiescent times correspond to periods in which the activity of the “inner engine” differs from its usual activity. Within the internal shocks model the shortest time scale corresponds to the variability of the “inner engine”, the quiescent periods time scale corresponds to periods of different activity of the “inner engine” and the duration of the burst corresponds to the total duration of activity of the “inner engine”.

We also find that the pulse duration is correlated to the preceding interval duration in all the analyzed bursts. A weaker but still significant correlation is found between the pulse duration and the following interval duration. This correlation was found only in 60% of the analyzed bursts.

While, different mechanisms can produce the observed lognormal distributions, the similarity between the pulse and interval distributions and the correlation between the pulses and neighboring intervals require that both processes would be determined by the same mechanism (or at least by two mechanisms governed by the same (source) physical parameters). This provides a new constraint on GRB models that should be considered in any model that attempts to reconstruct GRBs light curves.

For example, our results should be taken into account when choosing parameters in simulations of internal shocks light curves. Specifically, the parameters that determine the temporal characteristics of the relativistic wind should be chosen such that the resulting pulses and intervals would have lognormal distributions, with the appropriate parameters as well as quiescent times. In another paper (in preparation) we discuss the correlation between pulses and neighboring intervals within the internal shock model.

## Acknowledgment

This research was supported by a US-Israel BSF grant.


text:
# Implementing Unitarity in Perturbation Theory

## Abstract

Unitarity cannot be perserved order by order in ordinary perturbation theory because the constraint $`UU^{}=\mathrm{𝟏}`$ is nonlinear. However, the corresponding constraint for $`K=\mathrm{ln}U`$, being $`K=K^{}`$, is linear so it can be maintained in every order in a perturbative expansion of $`K`$. The perturbative expansion of $`K`$ may be considered as a non-abelian generalization of the linked-cluster expansion in probability theory and in statistical mechanics, and possesses similar advantages resulting from separating the short-range correlations from long-range effects. This point is illustrated in two QCD examples, in which delicate cancellations encountered in summing Feynman diagrams of are avoided when they are calculated via the perturbative expansion of $`K`$. Applications to other problems are briefly discussed.

Probability conservation is not maintained order by order in ordinary perturbation theory. This happens because the unitarity relation $`UU^{}=\mathrm{𝟏}`$ (for the time-evolution operator $`U`$) is nonlinear, whereas the constraints for the other exact conservation laws, such as energy, momentum, and charge, are linear. Order-by-order probability conservation can be restored if we expand instead $`K=\mathrm{ln}U_{n1}K_n/n`$, for then the unitarity constraint becomes the linear constraint $`K=K^{}`$. As long as every $`K_n`$ is kept to be anti-hermitean, $`U=\mathrm{exp}(K)`$ will be unitary no matter where the $`K`$-expansion is truncated. We shall henceforth refer to the perturbation theory in which $`K=\mathrm{ln}U`$ is expanded as the unitary perturbation theory. Mathematically, such an expansion is known as the Magnus expansion .

Besides unitarity, there is another reason to deal with the perturbative expansion $`K_n`$ of $`K=\mathrm{ln}U`$, rather than $`U=\mathrm{𝟏}+_{n1}U_n`$, because the former is the generalization of the linked-cluster expansion in probability theory and in statistical mechanics , with the associated advantages in separating short-range correlations from the long-range effects. It is also analogous to the expansion of the effective potential in quantum field theory. The added complication here is that we are dealing with non-commutative operators rather than c-number functions. In QCD, short-range correlations refer to correlations in time and in color. This point will be explained later with two QCD examples.

Unitary perturbation theory can also be regarded as an extension of the familiar eikonal approximation. For elastic scattering of a high-energy particle from a static potential, the scattering amplitude is diagonal in energy and the impact parameter. The scattering amplitude $`U`$ to all orders can be summed up to the form $`\mathrm{exp}(K_1)=\mathrm{exp}(K_1)`$, where $`K_1`$ is equal to the first Born approximation $`U_1`$, and $`\mathrm{}`$ is the matrix element in impact-parameter space . This eikonal approximation is also valid for electron-electron scattering when fermion loops are neglected . For high-energy scattering of a quark from a color potential, amplitudes of all orders can again be summed up into the form $`\mathrm{exp}(K)=\mathrm{exp}(K)`$, but now $`K`$ is a color matrix, with $`K_n`$ proportional to the $`n`$th nested commutator of the color matrices $`t_a`$ . In particular, $`K_n0`$ even for $`n>1`$. Therefore unitary perturbation theory can be thought of as a generalization of the non-abelian eikonal approximation.

The relation between $`U_n`$ and $`K_n`$ can be obtained by expanding the exponential $`\mathrm{exp}(K_n/n)`$ to the $`n`$th order. One obtains

$`U_n={\displaystyle \underset{(m)}{}}{\displaystyle \frac{1}{k!}}{\displaystyle \frac{1}{_{i=1}^km_i}}K_{m_1}K_{m_2}\mathrm{}K_{m_k},\text{ }\text{unk}`$ (1)
where the sum is taken over all partition $`(m)=(m_1m_2\mathrm{}m_k)`$ of the integer $`n`$. If $`U`$ is the time-evolution operator between times $`T^{}`$ and $`T`$, ordinary perturbation theory tells us

$`U_n`$ $`=`$ $`(i)^n{\displaystyle __n}d^ntH_I(t_1)\mathrm{}H_I(t_n),\text{ }\text{uh}`$ (2)
where $`_n`$ is the hyper-triangular integration region $`\{Tt_1\mathrm{}t_nT^{}\}`$. Together with (1), this may be used to derive a dynamical expression for $`K_n`$. In particular,

$`K_1`$ $`=`$ $`U_1=i{\displaystyle _T^{}^T}𝑑tH_I(t),`$ (3)
$`K_2`$ $`=`$ $`2U_2U_1^2=(i)^2{\displaystyle _T^{}^T}𝑑t_1{\displaystyle _T^{}^{t_1}}𝑑t_2[H_I(t_1),H_I(t_2)].\text{ }\text{k12}`$ (4)
The difference between ordinary and unitary perturbation theories can be illustrated by a simple spin-$`\frac{1}{2}`$ example , with interaction $`\frac{1}{2}\stackrel{}{\sigma }\stackrel{}{B}(t)`$, and an external magnetic field $`\stackrel{}{B}(t)=(B_{}\mathrm{cos}\omega t,B_{}\mathrm{sin}\omega t,B_0)`$. The longitudinal field $`B_0`$ causes a Zeeman splitting of magnitude $`B_0=\mathrm{}\omega _0`$; the transverse field $`B_{}`$ induces a transition between the lower and the upper states. Assuming the system to occupy the lower state initially, the probability to be in the upper state is plotted in Fig. 1, as a function of the scaled frequency variable $`\mathrm{\Omega }=\mathrm{}(\omega \omega _0)/B_{}`$ at the scaled time $`\tau =B_{}t/\mathrm{}=\pi `$. The solid, deahed, dash-dot, and dotted curves are respectively the exact solution, the first Born approximation $`P_1=U_1`$, the first unitary approximation $`UP_1=\mathrm{exp}(K_1)`$, and the second unitary approximation $`UP_2=\mathrm{exp}(K_1+K_2/2)`$. The notation $`\mathrm{}`$ represents the matrix element between the upper and the lower states. We see from Fig. 1 that the probability given by the Born approximation is larger than 1 for small $`|\mathrm{\Omega }|`$, violating unitarity, but both unitary approximations stay within the unitarity bound. As a function of $`\tau `$ at the resonance frequency $`\mathrm{\Omega }=0`$, the exact solution as well as the two unitary approximations are periodic, with a period $`2\pi `$ oscillating between 0 and 1, whereas the Born approximation grows monotonically like $`\tau ^2`$, making it a worse and worse approximation at large $`\tau `$.

For QCD, the unitary perturbation theory can be used to classify unitary parton-parton elastic amplitudes and to produce unitary model for Pomeron amplitudes .

Besides unitarity, it is often more profitable to compute $`U_n`$ via $`K_n`$, because $`K_n`$ contains short-time and color correlations via the appearance of the commutator of the $`H_I`$’s, instead of their product in $`U_n`$. I will illustrate this remark with two second-order QCD examples. Higher-order situtations will be discussed later.

The first example concerns high-energy electron-electron and quark-quark elastic scattering near the forward direction, with Mandelstam variables $`s`$ and $`t`$. The only difference between the two cases is the presence of a color matrix $`t_a`$ at the vertices of QCD diagrams, Fig. 2(a),(b),(c). The formulas for quark-quark amplitudes  are equally applicable to electron-electron amplitudes if we replace $`t_a`$ by 1. All amplitudes have a common factor $`is/2m^2`$ which will be factored out, leaving the rest to be $`U`$. Fig. 2(a) gives the first Born amplitude $`U_1=(ig^2/\mathrm{\Delta }^2)G_1`$, where $`\mathrm{\Delta }`$ is the (transverse) momentum transfer with $`t=\mathrm{\Delta }^2`$, and $`G_1=t_a\times t_a`$ is the color factor of Fig. 2a. The amplitude for Fig. 2(b) is $`(ig^4/2\pi )I_2\mathrm{ln}(se^{\pi i})G_2`$, and that for Fig. 2(c) is $`(ig^4/2\pi )I_2\mathrm{ln}s(G_2+G_1N_c/2)`$, where $`I_2=d^2k_{}/[(2\pi )^2k_{}^2(\mathrm{\Delta }k_{})^2]`$. The color factor $`G_2=G_1\times G_1`$ is that of Fig. 2(b), with $`N_c`$ being the number of colors in the theory. This number is to be set equal to 0 for QED. The second Born amplitude, being the sum of Figs. 2(b) and 2(c), is $`U_2=(ig^4/2\pi )I_2(\pi iG_2\mathrm{ln}sG_1N_c/2)`$. Note the cancelation of the $`\mathrm{ln}s`$ factor in the term proportional to $`G_2`$. The $`G_1`$ term consists of a color-octet exchange. It is not negligible compared to the first Born amplitude $`U_1`$ in the leading-log approximation where $`g^21`$ and $`g^2\mathrm{ln}s=O(1)`$. In fact it is the beginning of the reggeized gluon contribution. The color factor $`G_2=G_1\times G_1`$ consists of two color-octet exchanges, which can be decomposed into a color-octet and a color-singlet. Its coefficient in the impact-parameter representation is equal to the coefficient of $`G_1\times G_1`$ in $`\frac{1}{2}U_1^2`$, a factor of $`\mathrm{ln}s`$ down from the coefficient of the $`G_1`$ term. Its color-octet contribution is therefore negligible in the leading-log approximation, but its color-singlet component is leading and constitute the beginning of the Pomeron amplitude. For QED, the factor $`N_c`$ in front of $`G_1`$ is zero, and $`G_2`$ is replaced by 1. Hence $`U_2=\frac{1}{2}U_1^2`$ in the impact-parameter representation.

To summarize, a delicate cancelation of $`\mathrm{ln}s`$ occurs between Figs. 2(b) and 2(c) in QED, and in the color-singlet channel of QCD. This is reminiscent of the delicate cancellation of the volume factor in the grand partition function in statistical mechanics unless the linked-cluster expansion is used . As a result, the $`k`$ color-octet exchange amplitude in QCD is proportional to $`g^{2k}(g^2\mathrm{ln}s)^{2k}(k=1,2)`$, and it is this dependence that leads to the reggeon amplitudes in QCD. We shall now see that if we calculate $`U_2`$ via $`K_1^2`$ and $`K_2`$, then delicate cancelation is not needed, and the energy dependence needed for the reggeon structure becomes immediate.

Since $`K_1=U_1`$, it is given simply by Fig. 2(a). For $`K_2/2`$, it can be computed from Fig. 2(c), provided the product of the color matrices $`t_at_b`$ in the upper quark line is replaced by their commutator $`[t_a,t_b]`$ , thus making the color factor purely $`G_1`$. Otherwise, the result is identical to that of Fig. 2(c). In the impact-parameter representation where $`K_i`$ are diagonal, we obtain in this way the same result as before: $`K_2/2`$ is given by the $`G_1`$ coefficient of Fig. 2(c), or that of $`U_2`$, and $`K_1^2/2`$ is given by the $`G_2=G_1\times G_1`$ coefficient of $`U_2`$. The absence of the $`\mathrm{ln}s`$ factor in $`K_1^2/2`$ is not due to cancelation; it is simply a consequence that $`K_1`$ is independent of $`s`$. This is the same as in statistical mechanics when linked cluster expansion is used . The reggeized gluon contribution is now isolated in $`K_2/2`$ alone, without having to combine two Feynman diagrams to obtain it. Its octet origin can be traced back to the color commutator in the upper quark line, or the commutator structure of $`K_2`$ shown in (4). This gives $`K_2`$ a physical meaning as (part of) the one-reggeon amplitude.

The second example is pion-nucleon elastic scattering in the large-$`N_c`$ limit. The tree diagrams are shown in Fig. (d),(e), with pion-quark interaction given by a Hamiltonian of the form $`H_Iq^{}\stackrel{}{\mathrm{\Gamma }}\stackrel{}{\pi }q`$, where $`\stackrel{}{\mathrm{\Gamma }}`$ is a matrix that contains the isospin and spin information. The effective coupling at each vertex is proportional to $`H_I/\sqrt{N_c}`$, where $`\mathrm{}`$ represents the nucleon matrix element and $`1/\sqrt{N_c}`$ is the normalization factor needed for each external pion. Since a color-singlet nucleon contains $`N_c`$ quarks, $`H_IN_c`$, so the effective coupling at each vertex is of order $`\sqrt{N_c}`$, making each tree amplitude $`N_c`$. However, the term proportional to $`N_c`$ is canceled in the sum of the two diagrams, leaving a total amplitude $`U_2`$ of order unity in the large-$`N_c`$ limit. This can be understood  directly from (4) firstly because $`K_1=0`$ on account of energy-momentum conservation: an on-shell nucleon cannot absorb or emit a massive pion and remains on-shell. Thus $`U_2=\frac{1}{2}K_2`$. Since the commutator of two one-body operators is again a one-body operator, hence $`[H_1,H_2]q^{}[\stackrel{}{\mathrm{\Gamma }}_1\stackrel{}{\pi }_1,\stackrel{}{\mathrm{\Gamma }}_2\stackrel{}{\pi }_2]qN_c`$. Taking into account the normalization factor $`1/\sqrt{N_c}^2`$ for two pions, we conclude that $`U_2=\frac{1}{2}K_2N_c/N_c=1`$, as needed. Like the situation in Example 1, again there is no need to cancel the $`N_c`$ term explicitly, because the commutator structure of $`K_2`$ already provides for it. Physically, the one-body nature of $`K_2`$ tells us that the two pions must interact with the same quark in the nucleon.

With these two examples, it is clear that the commutator structure of $`K_2`$ is the key to the simplification. The question is whether $`K_n`$ still possesses such commutator structures for higher $`n`$. The answer is ‘yes’, though the detailed structure is increasingly more complicated. Let us define anti-hermitean operators $`C_n`$ to have the simple nested commutator structure,

$`C_n`$ $`=`$ $`(i)^n{\displaystyle __n}d^nt[H_1,[H_2,[\mathrm{},[H_{n1},H_n]\mathrm{}]]],\text{ }\text{cn}`$ (5)
where $`H_iH_I(t_i)`$. Then it can be shown  that

$`U_n={\displaystyle \underset{(m)}{}}{\displaystyle \underset{j=1}{\overset{k}{}}}{\displaystyle \frac{1}{_{i=j}^km_i}}C_{m_1}C_{m_2}\mathrm{}C_{m_k},\text{ }\text{xic}`$ (6)
where the sum is taken over all partitions $`(m)=(m_1m_2\mathrm{}m_k)`$ of the number $`n`$. In particular, $`U_1=C_1=K_1`$, $`U_2=\left(C_1^2+C_2\right)/2=\left(K_1^2+K_2\right)/2`$, $`U_3=\left(C_1^3+C_1C_2\right)/6+\left(C_2C_1+C_3\right)/3=K_1^3/6+\left(K_1K_2+K_2K_1\right)/4+K_3/3`$. From these relation, or more generally by comparing (1) with (6), one obtains the commutator structure of $`K_n`$ to be $`K_n=C_n+R_n`$, with $`R_1=R_2=0`$, $`R_3=[C_2,C_1]/4`$, and $`R_4=[C_3,C_1]/3`$. In general, $`R_n`$ is given by commutators of $`C_m`$, with increasingly complicated coefficients as $`n`$ increases.

We can use these relations between $`U_n`$ and $`C_m`$ to calculate quark-quark and electron-electron scatterings to higher orders, avoiding delicate cancelations, and obtaining directly the energy dependence necessary for reggeon structures in QCD, viz., an amplitude proportional to $`g^{2k}(g^2\mathrm{ln}s)^{nk}G_k`$ when $`k`$ color-octet objects are being exchange in the $`t`$-channel . They can also be used to obtain the correct amplitude  $`N_c^{1n/2}`$ for $`pi`$-nucleon inelastic scattering with $`n1`$ pions in the final state, without having to deal with delicate cancelation of $`n1`$ powers of $`N_c`$ encountered by summing Feynman tree diagrams . As mentioned at the beginnig, abelian and non-abelian eikonal approximations are special cases of these formulas. Infrared structure of QED can be obtained from these formulas, with $`K_1`$ containing the Bloch-Nordsieck result, $`K_2`$ giving rise to the Coulomb phase, and $`K_n=0`$ for $`n3`$ . The Landau-Pomeranchuk-Migdal effect  can also be cast into this form  with $`K_n=0`$ for $`n3`$. Other applications may require knowing the the diagrammatic rules for calculating $`K_n`$. For QED such rules are known . For tree diagrams in QCD involving an energetic particle such rules are also known , but the general case for QCD is still to be worked out.

So far we have emphasized situations when $`U`$ is the time-evolution operator with a known Hamiltonian. If the dynamics is unknown, we can still use $`K_n`$ in (1), or $`C_n`$ in (6), to parametrize the unitary dynamics. As long as $`K_n`$ is anti-hermitean, it is clear that $`U`$ is unitary. As long as $`C_n`$ is anti-hermitean, it can also be shown  that $`U`$ is unitary, though the proof is much more involved. A Wolfenstein-like parametrization of the CKM matrix to all orders can be deduced from these parametrizations. They can also be used to study the unitary matrix describing the overlap of the unperturbed and the perturb energy eigenfunctions .

In conclusion, we have given a number of examples to show that the unitary perturbation theory, in which $`K=\mathrm{ln}U`$ is expanded perturbatively, is in many ways better than the ordinary perturbation theory, where $`U`$ is directly expanded.

This research is supported in part by the Natural Sciences and Engineering Research Council of Canada, and the Fonds pour la formation de Chercheurs et l’Aide à la Recherche of Québec. I am indebted to Saad Ali, Hong-Mo Chan, Marc Grisaru, Greg Mahlon, Sheung-Tsun Tsou, Wu-ki Tung, and Kenneth Young for stimulating discussions.


text:
# Colour Reconnection in W Decays

## 1 Colour Reconnection

With some 450 pb-1 per experiment already recorded at $`\sqrt{s}=183202`$ $`\mathrm{GeV}`$, and more to come at higher energies, each of the four LEP experiments have selected up to now some 3200 WW $`qqqq`$ and some 2500 WW $`qq\mathrm{}\nu `$ candidate events. The mass and width of the W boson are measured from the kinematics of W decay products. Any energy-momentum exchange between W decay products not well simulated in Monte Carlo will affect the W mass and width measurement. Conventional MC’s treat the two $`q\overline{q}^{}`$ systems in a WW $`q\overline{q}^{}q\overline{q}^{}`$ event as independent. However, QCD interconnections, or colour reconnection (CR) can be expected . Perturbative CR effects are estimated to be small ; CR is a non-perturbative hadronization uncertainty that can only be studied in the context of various models.

CR models being used in these studies are those implemented in PYTHIA, ARIADNE and HERWIG. The models in PYTHIA, SK I, SK II and SK II’, are based on reconfiguration of overlapping or crossing strings . In the SK I model, the probability of reconnection is calculated as $`P_{\mathrm{reco}}=1e^{(k_iO)}`$, where $`O`$ is the overlap of two finite strings, and $`k_i`$ is a free parameter. In the SK II and II’ models, the string has no lateral dimension, and strings are reconnected when they cross. The ARIADNE models are based on rearrangement of colour dipoles if this reduces the string length . It should be noted that these models also affect LEP 1 data and are disfavoured from an OPAL study of the properties of quark- and gluon jets .

A reconfiguration of the colour flow is expected to change the (charged) particle multiplicity (typically decreasing it by 0.2 to 0.9 units), especially at low momentum, and more specifically between jets associated to the same W.

## 2 Multiplicities

### 2.1 Inclusive Charged Multiplicity

The charged multiplicity in WW events is measured by all four LEP experiments from charged tracks in the tracking system . The track multiplicity distribution is corrected to a charged particle multiplicity distribution by a matrix unfolding procedure. Alternatively, the multiplicity as a function of momentum (fragmentation function) or $`p_T`$ is determined and corrected bin-by-bin. The results are shown in Table 1.

The difference $`\mathrm{\Delta }<N_{ch}>=<N_{ch}^{4q}>2<N_{ch}^{2q}>`$ is also given in Table 1. DELPHI prefers to quote the ratio $`R=<N_{ch}^{4q}>/2<N_{ch}^{2q}>=0.990\pm 0.015\pm 0.011`$. Combining the results, it can be concluded that $`\mathrm{\Delta }<N_{ch}>`$ is consistent with 0 within an error of 0.3-0.4. A proper average is difficult due to differences in the definition, and the correlated systematics; the size of these correlated systematics (0.2-0.3), which is of the same size as the CR effects, limits the sensitivity of this method.

### 2.2 Fragmentation Function

By studying the particle multiplicity as a function of $`x_p=2p/\sqrt{s}`$, or $`\xi =\mathrm{log}(x_p)`$, one can study the low momentum region $`p<1`$ $`\mathrm{GeV}`$ where CR effects predominantly reside, but at the cost of reduced statistics. No significant effects at low $`x_p`$ are observed by any of the experiments.

### 2.3 Heavy Hadrons

Massive particles, like $`K^\pm `$ or (anti)protons, are more sensitive to CR effects than pions, by a factor 2 to 3. However, this is counterbalanced by the decreased statistics. DELPHI  and OPAL  have studied the production of heavy hadrons in $`qqqq`$ and $`qq\mathrm{}\nu `$ events and observe no significant differences.

## 3 Particle Flow

A more promising technique to study CR appears to be the study of the particle- or energy flow between jets from the same W and between different W’s, in analogy to studies of the string effect .

The construction of the particle flow is explained in Figure 1. A jet-finder is used to construct four jets in $`qqqq`$ events. Each pair of jets defines a plane onto which all reconstructed particles in the event are projected; for the energy flow weighted with the particle energy. In the preliminary studies submitted to this conference, L3  and ALEPH  use strong cuts on the angles between jets to select topologies with well separated jets and planar-like events, and obtain a selection efficiency of $`15`$%; OPAL  uses less restrictive cuts and a jet-pairing likelihood that gives a higher efficiency of $`42`$%, but selects also topologies with less clear separation between CR models. The flow is symmetrized with respect to the choice of jet-pairs, and particle angles between jets are rescaled between 0 and 1.

CR models indeed show a depletion of the particle flow between jets from the same W, and an increase between jets from different W’s, as naively expected. It is convenient to average the flow in the two regions between jets from the same W (regions j1-j2 and j3-j4), and to do the same for the flow in the two regions between jets from different W’s (regions j2-j3 and j4-j1). Subsequently, the ratio of these within-W/between-W flows is taken as a function of the rescaled particle angle.

L3 has studied 176 pb-1 of data taken at $`\sqrt{s}=189`$ $`\mathrm{GeV}`$. The ratio of the particle flow between jets from the same W and between jets from different W’s is shown in Figure 2, for data, PYTHIA without CR, and PYTHIA with SK I and GH, as a function of the rescaled angle $`\varphi _{\mathrm{resc}}`$. With this data sample only, a sensitivity of 3.5 $`\sigma `$ to SK I ($`k_i=1000`$) and 1.0 $`\sigma `$ to SK I ($`k_i=0.6`$) is reached. Varying the fraction of reconnected events, and calculating the $`\chi ^2`$ for the data-MC comparison, a fraction of $``$40% of reconnected events in the SK I model is favoured, and the No-CR scenario is disfavoured at 1.7 $`\sigma `$.

ALEPH has analyzed 347 pb-1 of data taken at $`\sqrt{s}=189200`$ $`\mathrm{GeV}`$. Their particle flow ratio is shown in Figure 3, for data, KORALW without CR, and KORALW with the SK I model for various values of $`k_i`$. Varying $`k_i`$, ALEPH finds the best data-MC agreement for $`k_i0.25`$, and puts a 1 $`\sigma `$ upper limit on $`k_i`$ of $`1.4`$ which corresponds to 45% of reconnected events at $`\sqrt{s}=189`$ $`\mathrm{GeV}`$.

OPAL has studied 183 pb-1 of data taken at $`\sqrt{s}=189`$ $`\mathrm{GeV}`$, and find sensitivities of 4.0 $`\sigma `$ for SK I ($`k_i=100`$), 1.1 $`\sigma `$ for SK I ($`k_i=0.9`$), 0.4 $`\sigma `$ for SK II and II’, and 0.5-1.8 $`\sigma `$ for AR2 and AR3. As a cross-check, OPAL uses the strong cuts like L3 and ALEPH, and observes slightly smaller sensitivities. The actual data is ambiguous, and prefers some reconnection in the default analysis, but no CR in the cross-check analysis. This, and in particular the role of the background, will be further studied.

DELPHI is also working on a similar analysis, but was not yet able to submit results to this conference.

With the full LEP 2 data sample, and combining all experiments, a further gain in sensitivity by a factor $`3.5`$ can be expected.

## 4 Effect on $`M_\mathrm{W}^{qqqq}`$

The estimates for $`\mathrm{\Delta }M_\mathrm{W}^{qqqq}`$ from the individual experiments calculated with their own Monte Carlo samples are summarized in Table 2 . A difference in reconnection probability in the SK II model can be expected from differences in the parton shower cutoff scale $`Q_0`$ , which ranges between 1.0 $`\mathrm{GeV}`$ (L3) and 1.9 $`\mathrm{GeV}`$ (OPAL).

Common samples of KORALW + JETSET Monte Carlo events, with (SK I) and without colour reconnection, have been generated for the four experiments; each experiment then passed these samples though detector simulation, event selection and analysis procedures. The resulting mass shifts found by the experiments were equal within errors, and a correlation between experiments of close to 100% was found. In view of this, further LEP collaboration will be needed to fully understand the differences in Table 2, especially in the ARIADNE estimates. For the LEP $`M_\mathrm{W}`$ combination, a CR systematic error of 50 $`\mathrm{MeV}`$ was used, fully correlated between experiments.

Estimates of the CR effect on the W width in the $`qqqq`$ channel range between +40 and +70 $`\mathrm{MeV}`$ in the SK models .

The studies of the particle flow between jets have proven to be sensitive to realistic CR model predictions. These studies will thus directly measure the amount of CR in the data. In models with a free parameter, such as SK I, this parameter can be measured from data; a calibration curve of mass shift versus $`k_i`$ can be used to estimate the CR systematic error. Already ALEPH, with a 1 $`\sigma `$ upper limit on $`k_i`$ of 1.4, puts a 1 $`\sigma `$ upper limit on $`\mathrm{\Delta }M_\mathrm{W}^{qqqq}`$ in the SK I model of 40 $`\mathrm{MeV}`$ . With the full LEP 2 data sample and combining all experiments, the CR systematic error on $`M_\mathrm{W}`$ is likely to be below that, and, most important, actually measured from data.


text:
# Stochastic resonance induced by random fields in ferroelectrics

## Abstract

We propose a mechanism of stochastic resonance induced by spatial noise (electric random fields) in ferroelectrics. The calculations demonstrate the characteristic of the stochastic resonance: certain random field can increase the susceptibility, the phonon occupancy, and the overcooled temperature. This mechanics of stochastic resonance may be general in other disordered systems.

Stochastic resonance (SR) is a nonlinear phenomenon in which the degree of order rather than disorder can be enhanced by the presence of optimized random noise. It was originally proposed to explain the periodicity in the recurrences of the Earth’s ice ages as the result of stochastic and weak periodic forces acting in concert on a bistable global climate model. Experimentally, SR was first demonstrated in a Schmitt trigger circuit and a bistable ring-laser system. Since then stochastic resonance has been reported in various systems from crayfish to SQUID and spin system. Variations and extensions of the classical definition of SR were also developed to include aperiodic (e.g., dc or wideband) signals. SR is particularly interesting for neurobiological systems because it provides a mechanism for such systems to detect and process weak signals. In applications, SR may be useful in physical, technological and biomedical contexts.

Stochastic resonance usually occurs in systems with multiple stable states. On the other hand, ferroelectrics are known as multistable systems since the polarization can emerge along more than one directions. In the so called relaxor ferroelectrics (relaxors) such as Pb(Mg1/3Nb2/3)O3 (PMN), there exists spatial random electric fields induced by point charge defects, nano-domain textures, and the direct interactions of polar ions, etc. Some specific characteristics of the relaxors were interpreted in terms of different models adopting random fields. Recently, it was shown that an optimal intensity of random fields can promote the spontaneous appearance of a first-order phase transition by increasing the overcooled temperature of ferroelectrics. Well then a question naturally arises: Can the random fields in ferroelectrics produce a certain stochastic resonance?

To resolve the above problem, in this paper we investigate a ferroelectric system with random fields, and show that a stochastic resonance does occur in this case.

The basic picture of stochastic resonance is generally illustrated using the following simple example: Imagine a particle subject to friction, moving in a double-well potential. A weak periodic signal force is applied on the particle (here, “weak” means that the signal force is far too small to excite the particle over the barrier between the wells). Meanwhile, random noise is applied to tilt the potential and could result in irregular switching of the particle between the wells. With the help of certain noise, the weak periodic signal can effectively make the particle hop over the barrier, and produce the maximal signal-to-noise ratio (SNR) at a nonzero level of input noise. This phenomenon is called stochastic resonance. In general, it is assumed that the signal period is longer than some characteristic intrawell relaxation time for the system. According to the statistical mechanics, the time average of a mechanical quantity of a system is equal to the ensemble average (ergodic hypothesis). So the above stochastic resonance can be constructed in another way: many particles moving in different double-well potentials, respectively. In this case, the applied noise does not change in time, but changes in space, i.e., double-well potentials tilt in random slopes. This is just the physical picture of the ferroelectric system with random fields where the tilt of potential is caused by random electric fields. Therefore, a stochastic resonance with “spatially extended” noise can be reasonably expected in ferroelectrics.

For the ferroelectrics with the second-order phase transition, it was shown that the critical temperature monotonously declines with increasing random fields. In this case, even if a phenomenon of stochastic resonance appears, it seems not to lead to interesting practical applications. Therefore, in this paper, we focus our study on the ferroelectrics with the first-order phase transition, i.e., with a triple-well potential (See Fig. 1). For this system, when the particle (polar ion) jumps from the middle well to the sideward wells, a paraelectrics-to-ferroelectrics phase transition takes place. Since optimal noise can enhance the hopping of particle between wells, a random-field-induced ferroelectric transition and nonmonotonic behaviors might be observed.

A microscopic Hamiltonian describing order-disorder ferroelectrics is introduced as:

$$H=\left[\begin{array}{cccccccc}Ep_0& \frac{\mathrm{\Omega }}{2}& 0& 0& 0& 0& 0& 0\\ \frac{\mathrm{\Omega }}{2}& Ep_0& 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{1}{3}Ep_0& \frac{\mathrm{\Omega }}{2}& 0& 0& 0& 0\\ 0& 0& \frac{\mathrm{\Omega }}{2}& \frac{1}{3}Ep_0& 0& 0& 0& 0\\ 0& 0& 0& 0& \frac{1}{3}Ep_0& \frac{\mathrm{\Omega }}{2}& 0& 0\\ 0& 0& 0& 0& \frac{\mathrm{\Omega }}{2}& \frac{1}{3}Ep_0& 0& 0\\ 0& 0& 0& 0& 0& 0& \frac{1}{3}Ep_0& \frac{\mathrm{\Omega }}{2}\\ 0& 0& 0& 0& 0& 0& \frac{\mathrm{\Omega }}{2}& \frac{1}{3}Ep_0\end{array}\right],$$
(1)
where $`E`$ is the equivalent field, and $`p_0`$ is the dipole moment of the ion in a local polar position. $`\mathrm{\Omega }`$ is the tunneling frequency between wells. Under a mean field approximation,

$$Ep_0=J\frac{p}{p_0}+E_{\mathrm{ext}}p_0+E_{\mathrm{rand}}p_0,$$
(2)
where $`p`$ is the thermal average of the dipole moment, and $`J`$ is the coupling energy between polar ions. $`E_{\mathrm{ext}}`$ is the applied external electric field, and $`E_{\mathrm{rand}}`$ is the internal random field in the system which is assumed to have a Gaussian distribution with a width $`\sigma _e`$. Hamiltonian in Eq. (1) produces a second-order phase transition for $`\mathrm{\Omega }<0.56J`$ and a first-order one for $`0.56J<\mathrm{\Omega }<2J`$ when there is no random field.

From Eq. (2) one can see that the uniform variable $`E_{\mathrm{ext}}`$ is the input signal, and the inhomogeneous $`E_{\mathrm{rand}}`$ is the input random noise. The average polarization $`p`$ is the output signal, and appears in the right hand of Eq. (2) as the feedback input through the ferroelectric interaction between polar ions. When there is no external input ($`E_{\mathrm{ext}}=0`$), the feedback mechanism can preserve a nonzero output at low temperatures and the system stays in a ferroelectric state. The over-cooled temperature $`T_{}`$ is used to measure the spontaneous phase transition in a cooling process. In Fig. 2, $`T_{}`$ is given as function of the intensity of noise ($`\sigma _e`$) when $`\mathrm{\Omega }=0.63J`$. It can be seen that $`T_{}`$ increases first and then declines with increasing random field, which is consistant with previous calculation. $`T_{}`$ reaches its maximum at $`\sigma _e=0.23J/p_0`$. This means that a nonzero random field can induce the spontaneous appearance of the ferroelectric state, i.e., enhance the output signal. Stochastic resonance really exists in the ferroelectrics.

When the temperature is higher than the maximum over-cooled temperature, the system stays in a paraelectric state for any noise. The effect of random fields can then be understood by investigating the response of the system to a weak external field. The ratio of the output and input signals is measured by the dielectric susceptibility,

$$\chi =\frac{P}{\epsilon _0E_{ext}}=\frac{N_0p}{\epsilon _0E_{ext}},$$
(3)
where $`P`$ is the polarization, and $`N_0`$ is the number of polar ions in unit volume. The calculation result is depicted in Fig. 3. It shows that appropriate noise can greatly enhances the output signal and results in a peak of the susceptibility, which is just the characteristic signature of stochastic resonance.

In the paraelectric state with no external field, the polarization does not keep zero but fluctuates due to the thermal movement. The polarization fluctuation can be described in terms of the ferroelectric soft mode and the phonon. The dependence of the soft-mode frequency and the phonon occupancy on the random fields is demonstrated in Fig. 4. Nonmonotonic behaviors are observed for both quantities. Stochastic resonance causes the decrease of the soft-mode frequency and the increase of the phonon occupancy, i.e., makes it easier for the thermal fluctuation of polarization to take place.

For ferroelectrics, the random fields are mainly determined by the impurity and the microstructure imhomogeneity in the system. If there are more charge impurities or the microstructure is more inhomogeneous, the random fields are stronger. The features of imhomogeneity include grain size, porosity, grain boundary, chemical imhomogeneity, and domain/twin structure, etc. The imhomogeneity can be controlled by thermal treatments in some degrees. So the stochastic resonance can be realized through the doping and thermal treating. Further experiments in this aspect are needed to testify the theory.

In summary, we have demonstrated a specific stochastic resonance in ferroelectric systems, which is not produced by the temporal noise, but by the spatial noise (random electric fields). The characteristics of stochastic resonance are revealed via the calculations of the over-cooled temperature, the dielectric susceptibility, the soft-mode frequency, and the phonon occupancy, which all show nonmonotonic dependence on the noise. The spatial stochastic resonance may be useful in materials applications as well as in theoretical interests. For example, the stability of paraelectric and ferroelectric phases can be adjusted by modulating the over-cooled temperature, and the dielectric susceptibility of materials can be increased by selecting appropriate amount of random fields. The mechanics of the spatial stochastic resonance revealed in this paper might also be applicable in a number of inhomogeneous systems whose properties are not well understood.

This work was supported by the Chinese National Science Foundation (Grant NO. 59995520) and State Key Program of Basic Research Development (Grant No. G2000067108).