# Annealing optimisation of Mixed ZX Phase Circuits

Stefano Gogioso Hashberg Ltd<EMAIL_ADDRESS>Quantinuum Ltd Richie Yeung
Quantinuum Ltd<EMAIL_ADDRESS>
###### Abstract

We present a topology-aware optimisation technique for circuits of mixed ZX
phase gadgets, based on conjugation by CX gates and simulated annealing.

##  1 Introduction

In this work, we build upon the Master’s thesis by one of the authors [19] and
present a topology-aware optimisation technique for circuits of mixed ZX phase
gadgets, based on conjugation by CX gates and simulated annealing. The basic
rules on CX conjugation of phase gadgets have previously appeared in the
literature [12, 3, 4, 8, 15]—which features other topology-aware
techniques—and are used by both the PyZX library [11, 10] and the
t$|$ket$\rangle$ compiler [17, 16]. We test the performance of our
optimisation technique on random circuits of mixed ZX phase gadgets. An open-
source implementation is made available as part of the Python library pauliopt
[7], which is also used to generate the circuit figures in this work.

With NISQ applications in mind, our optimisation target is the number
$\text{count}_{\text{CX}}\left(P\right)$ of nearest-neighbour (NN) CX gates
required to implement a mixed ZX phase circuit $P$ on a given topology, making
the following assumptions on compilation:

  * •

Each phase gadget will be compiled to a single-qubit rotation conjugated by
trees of CX gates.

  * •

Long-range CX gates in the topology will be compiled to double-ladders of NN
CX gates.

In these circumstances, the CX count of an individual phase gadget depends on
the mutual distances between its legs. Specifically, it is computed by finding
a minimum spanning tree with weights related to the distance of leg qubits in
the given topology.

Our optimisation technique is based on the observation that conjugating a
given mixed ZX phase circuit $P$ with an arbitrary CX circuit $C$ has the
effect of changing the legs of the individual gadgets, without otherwise
altering their angles or position within the circuit. The resulting circuit
$C^{\dagger}\left(P\right):=C^{\dagger}\circ P\circ C$ is again a mixed ZX
phase circuit, related to the original circuit by:

$P=C\circ C^{\dagger}\left(P\right)\circ C^{\dagger}$

However, the different gadget legs in $C^{\dagger}\left(P\right)$ might result
in a lower overall CX count:

$2\cdot\text{count}_{\text{CX}}\left(C\right)+\text{count}_{\text{CX}}\left(C^{\dagger}\left(P\right)\right)$

Finding a CX circuit $C_{opt}$ which (approximately) minimises the overall CX
count above results in an optimised implementation of our initial mixed ZX
phase circuit $P$, expressed as conjugation of another mixed ZX phase circuit
$C_{opt}^{\dagger}\left(P\right)$ by the CX circuit $C_{opt}$:

$P=C_{opt}\circ C_{opt}^{\dagger}\left(P\right)\circ C_{opt}^{\dagger}$

We use simulated annealing to explore the space of conjugating CX circuits and
find an approximately optimal $C_{opt}$, starting from the empty CX circuit.
We fix a number of layers and explore the space by progressively introducing
or removing NN CX gates within those layers, altering the overall CX count in
the process. In the early stages of optimisation, the annealing “temperature”
$t$ is high and we are free to explore the configuration space, even at the
expense of a—hopefully temporary—increase in overall CX count. As the
optimisation progresses, the temperature $t$ is lowered and the search becomes
progressively more greedy, accepting changes which increase the CX count by
$\Delta$ with progressively lower probability:

$\text{Prob}\left(\text{accept CX count increase
$\Delta$}\right)=\dfrac{1}{2^{\Delta/t}}$

Our technique is applicable to parametric circuits, such as the ansatzes used
in quantum machine learning, adiabatic quantum computation [6] and quantum
approximate optimisation [5]. Such ansatzes typically consist of a large
number of repeating layers, using the same phase gadgets with possibly
different parameters. In this context, a key observation about our technique
is that a single conjugating CX circuit is sufficient to optimise all layers
at once. The number $K$ of layers can then be absorbed into the overall CX
count calculation, with a single layer $L$ being optimised:

$2\cdot\text{count}_{\text{CX}}\left(C\right)+K\cdot\text{count}_{\text{CX}}\left(C^{\dagger}\left(L\right)\right)$

As a consequence, the optimisation cost scales with the size of $L$,
regardless of $K$, making our technique especially suited for application in
the aforementioned domains.

##  2 ZX Phase Gadgets

A _Z phase gadget_ on $n$ qubits is the exponential of an imaginary scalar
multiple of an element of the Pauli group $P_{n}$ on $n$ qubits involving only
the Z and I Pauli matrices (and with scalar factor $+1$):

$\exp\left(i\theta\bigotimes_{q=0}^{n-1}G_{q}\right)\hskip
14.22636pt\text{where}\hskip 14.22636pt\begin{cases}G_{q}=\text{Z}&\text{ if
qubit $q$ is a leg}\\\ G_{q}=\text{I}&\text{ otherwise}\end{cases}$

The element $\bigotimes_{q=0}^{n-1}G_{q}\in P_{n}$ is known as the _generator_
, the parameter $\theta\in U(1)$ is known as the _angle_ , and the qubits $q$
where $G_{q}=\text{Z}$ are known as the _legs_ of the gadget. An _X phase
gadget_ on $n$ qubits is analogously defined, using the Pauli matrix X in
place of Z. Below are two examples of phase gadgets on 3 qubits:

$\exp\left(i\frac{\pi}{4}\text{Z}\otimes\text{I}\otimes\text{Z}\right)\hskip
85.35826pt\exp\left(i\frac{3\pi}{4}\text{X}\otimes\text{X}\otimes\text{X}\right)$

Phase gadgets can be represented in the ZX calculus by attaching spiders for
the gadget basis to each leg and joining them at a “hub”, formed by spiders of
alternating colours and including the angle ($\theta=\frac{\pi}{4}$ and
$\theta=\frac{3\pi}{4}$ respectively) as a spider phase. Below are the same
two examples represented in the ZX calculus (Z phase gadget on the left, X
phase gadget on the right):

###  2.1 Conjugation by CX gates

Z and X phase gadgets are particularly well-behaved under conjugation by CX
gates: they either keep the same legs, they gain a leg, or they lose a leg,
but basis and angle are always left unchanged. Instead of thinking of CX gates
in terms of _control_ and _target_ , we think of them as having a Z qubit (the
control) and an X qubit (the target). Z phase gadgets are left unchanged under
conjugation if the X qubit is not a leg:

=

If the X qubit is a leg and the Z qubit is not a leg before conjugation, the Z
qubit becomes a leg after conjugation:

=

If the X qubit is a leg and the Z qubit is a leg before conjugation, the Z
qubit is no longer a leg after conjugation:

=

The behaviour of X phase gadgets is obtained by inverting Z and X everywhere.
X phase gadgets are left unchanged under conjugation if the Z qubit is not a
leg:

=

If the Z qubit is a leg and the X qubit is not a leg before conjugation, the X
qubit becomes a leg after conjugation:

=

If the Z qubit is a leg and the X qubit is a leg before conjugation, the X
qubit is no longer a leg after conjugation:

=

There are several different ways to obtain the conjugation rules above. A
simple one goes through the observation that conjugation commutes with matrix
exponentiation:

$A\left(\exp\left(i\theta\bigotimes_{q=0}^{n-1}G_{q}\right)\right)A^{\dagger}=\exp\left(i\theta\;A\left(\bigotimes_{q=0}^{n-1}G_{q}\right)A^{\dagger}\right)$

From this, it is enough to look at conjugation by CX for a handful of Pauli
group elements:

$\begin{array}[]{ccc}\text{cx}_{0,1}(\text{Z}\otimes\text{I})\text{cx}_{0,1}=\text{Z}\otimes\text{I}&\;\;\text{cx}_{0,1}(\text{I}\otimes\text{Z})\text{cx}_{0,1}=\text{Z}\otimes\text{Z}&\;\;\text{cx}_{0,1}(\text{Z}\otimes\text{Z})\text{cx}_{0,1}=\text{I}\otimes\text{Z}\\\
\text{cx}_{0,1}(\text{I}\otimes\text{X})\text{cx}_{0,1}=\text{I}\otimes\text{X}&\;\;\text{cx}_{0,1}(\text{X}\otimes\text{I})\text{cx}_{0,1}=\text{X}\otimes\text{X}&\;\;\text{cx}_{0,1}(\text{X}\otimes\text{X})\text{cx}_{0,1}=\text{X}\otimes\text{I}\\\
\end{array}$

###  2.2 Mixed ZX Phase circuits

A _mixed ZX phase circuit_ is, as the name suggests, a quantum circuit
consisting of a mix of Z phase gadgets and X phase gadgets. Examples of mixed
ZX phase circuits include the ansatzë used by adiabatic quantum computation
[6], the quantum approximate optimisation algorithm [5], as well as many
quantum machine learning techniques [19]. The pauliopt library allows these
circuits—called _phase circuits_ , for short—to be created directly from a
sequence of phase gadgets with given angles and legs:

⬇

circ = PhaseCircuit(3)

circ >>= Z(π/2)@{0,1}

circ >>= X(π)@{0}, Z(π/2)@{1}

circ >>= X(-π/4)@{1,2}

Utility methods are also made available to translate common gates into phase
gadgets:

⬇

circ = PhaseCircuit(3)

circ.ccz(0,1,2)

A generalisation of Euler’s decomposition to $n$ qubits (technically, a
consequence of Zassenhaus formula [14, 2]) implies that every unitary on $n$
qubits can be expressed using a bounded number of Z and X phase gadgets,
meaning that circuits of mixed ZX phase gadgets are universal.

##  3 Annealing optimisation

###  3.1 Topology-aware cost of phase gadgets

The optimisation method targets the CX count for phase gadget implementation
on a given topology, which we refer to as its _cost_. The cost of a gadget
relies on its implementation using a balanced tree of CX gates between its
legs [3], individually converted into double ladders of nearest-neighbour (NN)
CX gates (without further simplification). The balanced tree of minimum CX
count is obtained using Prim’s algorithm for minimum spanning trees, using the
following weights for any two distinct leg qubits $q_{i}$ and $q_{j}$ in the
phase gadget:

$w(q_{i},q_{j}):=4d(q_{i},q_{j})-2$

where $d(q_{i},q_{j})$ is the graph distance between qubits $q_{i}$ and
$q_{j}$ in the given topology. For example, consider the following 4-legged
gadget on a 3-by-3 grid topology:

The minimum spanning tree implementation for this gadget uses two (adjoint)
balanced trees of 3 CX gates each—between qubits $\\{0,3\\}$ (NN CX count: 1),
between qubits $\\{3,6\\}$ (NN CX count: 1), and between qubits $\\{3,5\\}$
(NN CX count: 3)—for a total implementation cost of 10 NN CX gates. In this
case, it is easy to see that every alternative would have been sub-optimal,
because the only common nearest neighbour between qubits 0 and 6 on the grid
is qubit 3: connecting both qubits to a qubit at distance 2 from both would
have incurred a NN CX count of at least $3+3=6$, while connecting one qubit to
a different nearest neighbour of the other (e.g. connecting 0 and 6 to 1)
would have incurred a NN CX count of at least $5+1=6$, both already exceeding
the optimal overall NN CX count of 5. The following snippet produces the phase
gadget in Qiskit [18] using the above minimum spanning tree:

⬇

from qiskit.circuit import QuantumCircuit, Parameter

cx_tree = QuantumCircuit(4)

cx_tree.cx(0, 2)

cx_tree.cx(1, 2)

cx_tree.cx(2, 3)

gadget = QuantumCircuit(8)

gadget.compose(cx_tree, qubits=(0,6,3,5), inplace=True)

gadget.rz(Parameter(”θ”), 5)

gadget.compose(cx_tree.inverse(), qubits=(0,6,3,5), inplace=True)

$\mapsto$

###  3.2 Conjugation by CX circuits

As discussed above, conjugation of mixed ZX phase circuits by CX circuits
preserves the order and angles of the gadgets, but modifies their legs. For
adequate choices of nearest-neighbour CX circuits, the increase in CX count
from the conjugating CX gates is counterbalanced by a decrease in
implementation cost for the phase gadgets on a given topology, due to leg
changes and reduction. The overall 2-qubit gate count can then be used as a
cost function for global optimisation methods (such as simulated annealing).
Here is example of a 10-gadget circuit on 4 qubits:

On a cycle topology, the circuit above has a CX count of 42. Below is an
equivalent circuit with a CX count of 14, obtained by conjugation from two
layers of nearest-neighbour CX gates:

Here is a second example, of a 20-gadget circuit on 9 qubits:

On a 3-by-3 grid topology, the circuit above has a CX count of 142. Below is
an equivalent circuit with a CX count of 126, obtained by conjugation from two
layers of nearest-neighbour CX gates:

The two layers of conjugating CX gates for each one of the two example above,
arranged on the corresponding qubit topologies, are displayed below:

###  3.3 Repeating layers

A remarkable feature of our optimisation technique concerns its use with
circuits composed of repeating layers. For such circuits, it is enough to
simplify an individual layer and then repeat it: the conjugating CX gates
between repeating layers cancel out in pairs, leaving a single pair of
conjugating CX block for the entire circuit.

This means that the per-layer CX cost reduction is multiplied by the number of
layers, while the CX cost for the conjugating blocks is fixed: as the number
of layers grows, the cost of larger conjugating CX blocks is offset by the
increased savings on the layers, at no additional computational expense for
the optimisation itself. As a concrete example, consider the following
4-gadget layer:

On a line topology with 5 qubits, the initial CX count for this layer is 22.
If optimisation is run for a single layer, the CX count is reduced by 18%,
from 22 to 18:

With 2 layers, the CX count is reduced by 36%, from 44 to 28:

With 5 layers, the CX count is reduced by 47%, from 110 to 58:

The 8 conjugating CX gates are independent of the number of layers, leaving 10
CX gates per simplified layer from an initial count of 22. In the limit of a
large number of layers, the relative CX count reduction for this conjugating
CX block approaches 12/22, or 54.5%.

This feature of our optimisation technique makes it particularly interesting
for applications in quantum machine learning (QML), adiabatic quantum
computation and the quantum approximate optimisation algorithm (QAOA), where
circuit ansatzë often feature the aforementioned repeating layers structure.
In particular, note that the technique is insensitive to the specific angles
used in the repeated layers: the angles can be different, and even parametric
(as indeed shown in the examples above). This means that an ansatz circuit
only need be optimised once, before assigning specific values to its
parameter. For a mixed ZX phase circuit formulation of several such ansatzë,
we refer the reader to [19].

###  3.4 Configuration space

The configuration space explored by our optimisation algorithm on a given
topology consists of all circuits of nearest-neighbour CX gates with a fixed
number of layers, which we refer to as the _CX blocks_. At any given point in
configuration space, i.e. at any such CX block $C$, the cost for a phase
circuit $P$ is given by sum of:

  * •

the CX count for $C^{\dagger}$;

  * •

the CX counts for the individual gadgets of $C^{\dagger}(P)$, computed on the
given topology;

  * •

the CX count for $C$;

where $C^{\dagger}(P):=C^{\dagger}\circ P\circ C$ is the phase circuit
obtained from $P$ through conjugation by $C^{\dagger}$.

For fixed topology and number of layers, we say that a CX block $C^{\prime}$
is obtained from another CX block $C$ by _flipping_ a CX gate
$\text{cx}_{i,j}$ on a layer $l$ if either:

  * •

the gate $\text{cx}_{i,j}$ is present at layer $l$ in $C$, and $C^{\prime}$ is
obtained from $C$ exactly by removing $\text{cx}_{i,j}$ from layer $l$;

  * •

there is no gate at layer $l$ in $C$ which is incident to either qubit $i$ or
qubit $j$ (or both), and $C^{\prime}$ is obtained from $C$ exactly by adding
$\text{cx}_{i,j}$ to layer $l$.

We write $C\stackrel{{\scriptstyle l;i,j}}{{\longrightarrow}}C^{\prime}$ to
denote this fact; because CX gates are self-inverse, this is equivalent to
$C\stackrel{{\scriptstyle l;i,j}}{{\longleftarrow}}C^{\prime}$. We consider
two CX blocks to be nearest neighbours in configuration space if one can be
obtained from the other by a single CX gate flip: this endows the
configuration space with the structure of an undirected graph, where edges are
labelled by the triples $(l;i,j)$ describing the CX gate flips. Our
optimisation algorithm performs a random walk on this graph, starting from the
empty CX block and attempting to (globally) minimise the cost for a given
phase circuit.

###  3.5 Simulated Annealing

Simulated annealing is a global optimisation method which explores a discrete
problem-dependent configuration space $\mathcal{X}$, in an attempt to find an
approximate optimum for a an arbitrary cost function $f$. At each step of the
optimisation, the algorithm sits at some configuration $x\in\mathcal{X}$ and
selects a uniformly random “direction” to explore; that is, it samples a
uniformly random $x^{\prime}\in N(x)$, where
$N(x)\subseteq\mathcal{X}\backslash\\{x\\}$ a bounded set of nearest
neighbours of $x$ in configuration space. The algorithm the computes
$f(x^{\prime})$ and compares it to $f(x)$: if $f(x^{\prime})\leq f(x)$, then
the algorithm moves to $x^{\prime}$ deterministically; otherwise, it moves to
$x^{\prime}$ with a probability which decreases exponentially in the
difference $\Delta:=f(x^{\prime})-f(x)$:

$\text{Prob}\left(x\rightarrow
x^{\prime}\right)=\exp\left(-\frac{\Delta\log(2)}{t}\right)$

The parameter $t>0$ is known as the _temperature_ : it decreases from an
initial value $t_{0}$ to a final value $t_{N}\approx 0$ over $N$ iterations,
in a way specified by the “temperature schedule”. The temperature schedule for
the annealing determines the probability of accepting changes which increase
the cost function: an increase of $t\log_{2}(1/p)$ has probability $p$ of
being accepted. At the start of the annealing, when the temperature is high,
the random walk is allowed to explore the configuration space, even at the
cost of an increase in cost. As iterations progress and the temperature is
lowered, the random walk progressively favours those steps which decrease the
cost. Towards the end, the probability of accepting steps which increase the
cost becomes vanishingly small, and the algorithm is reduced to a greedy
optimisation.

###  3.6 Simulated Annealing for Mixed ZX Phase Circuit Optimisation

For our problem, the simulated annealing configuration space consists of a
fixed number of layers of NN CX gates on the given topology, with nearest
neighbours in configuration space determined by CX gate flips in any one of
the layers. Given a mixed ZX phase circuit, a topology and a fixed number of
layers for the CX blocks, we construct an optimized circuit container;
initially, this consists of the phase circuit conjugated by empty CX blocks.
To minimise the cost of over the space of CX blocks, we perform a random walk
using simulated annealing [9], for a given number of steps/iterations and
following a given temperature schedule. At the end of the annealing, the
optimized circuit contains a simplified phase circuit conjugated by
(typically) non-empty CX blocks. Every iteration of our simulated annealing
algorithm proceeds as follows:

  1. 1.

We obtain the current temperature $t$ from the temperature schedule using the
current iteration number (e.g. linearly interpolating from initial temperature
to final temperature).

  2. 2.

We select a random CX gate flip for the current block, i.e. select a random
neighbour.

  3. 3.

We flip the selected CX gate (see Subsection 3.7 below).

  4. 4.

We calculate the CX count $c_{new}$ for the new circuit and compare it to the
previous CX count $c_{old}$, by computing $\Delta:=c_{new}-c_{old}$ and
sampling a uniformly random number $r\in[0,1]$:

     * •

if $r\leq\exp\left(-\frac{\Delta\log(2)}{t}\right)$, and in particular if
$c_{new}\leq c_{old}$, the move is accepted and we go to the next iteration
from the new CX block;

     * •

otherwise, the move is undone—by flipping the selected CX gate again—and we go
the next iteration from the old CX block.

Our selection of temperature schedule depends on the shape of the temperature
curve (linear, geometric, reciprocal or logarithmic) and two parameters: the
amount of time we want to spend exploring the configuration space at the start
of the annealing and the amount of time we want to spend performing greedy
optimisation at the end of the annealing. As a convention, we take the
temperature value $t=2m+2$ (where $m$ is the number of phase gadgets) to
denote the end of the exploration phase and the temperature value $t=2$ to
denote the start of the greedy optimisation phase. Because the CX gates are
nearest-neighbour, $2m+2$ is the maximum possible increase in CX count from a
single CX gate flip: 2 CX gates added plus a leg added to all gadgets, at a
worst-case CX count increase of 2 per gadget. On the other hand, $2$ is the
minimum possible (positive) increase: CX count always changes by multiples of
2. This defines the exploration phase as the time where every CX gate flip has
a probability higher than $50\%$ of being accepted, and the greedy
optimisation phase as the time where every increase in CX count has a
probability lower than $50\%$ of being accepted. As an example, below is a
sketch of a linear temperature schedule over $N$ iterations:

###  3.7 Flipping CX gates

Because simulated annealing requires thousands of iterations, optimised
execution of CX gate flips is paramount. Naively, one could simply maintain
the current CX block $C$ and compute the conjugated circuit from scratch at
every iteration: this requires conjugating every gadget in the phase circuit
by every CX gate in the CX block, even though the CX block itself has only
changed by one gate flip. In this work, we adopt a different strategy, where
we maintain both the current CX block $C$ and the current conjugated phase
circuit $C^{\dagger}(P)$. Imagine we wish to perform a gate flip
$C\stackrel{{\scriptstyle l;i,j}}{{\longrightarrow}}C^{\prime}$, e.g.
performing the flip $C\stackrel{{\scriptstyle
3;5,4}}{{\longrightarrow}}C^{\prime}$ boxed in dark blue below.

To perform the flip, we go through the following 4 steps:

  1. 1.

Starting from layer $l-1$ and moving inwards towards layer $0$, we create the
(partially ordered) set $(l;i,j)\\!\\!\downarrow$ of all CX gates that are
strictly in the “past” of gate $\text{cx}_{i,j}$ at layer $l$. These are boxed
and highlighted below.

  2. 2.

We undo all the CX gates in $(l;i,j)\\!\\!\downarrow$, conjugating the phase
circuit by each one, working our way outwards from layer $0$ (closest to the
phase circuit) up to layer $l-1$ (just below the gate we’ll flip).

  3. 3.

We flip gate $\text{cx}_{i,j}$ at layer $l$, conjugating the phase circuit by
$\text{cx}_{i,j}$ (adding $\text{cx}_{i,j}$ to layer $l$ if it wasn’t there,
removing it from layer $l$ if it was). Below, we remove $\text{cx}_{5,4}$ from
layer $3$.

  4. 4.

We redo all the CX gates in $(l;i,j)\\!\\!\downarrow$, conjugating the phase
circuit by each one, this time working our way inwards from layer $l-1$ down
to layer $0$.

To define the poset $(l;i,j)\\!\\!\downarrow$ of past gates formally, consider
the partial order on the gates $(l;x,y)$ of the CX block obtained by
transitive-reflexive closure of the following relation:

$(l_{1};x_{1},y_{1})<(l_{2};x_{2},y_{2})\Leftrightarrow l_{1}<l_{2}\text{ and
}\\{x_{1},y_{1}\\}\cap\\{x_{2},y_{2}\\}\neq\emptyset$

The poset $(l;i,j)\\!\\!\downarrow$ consists of all gates strictly less than
$(l;i,j)$ in this partial order.

##  4 Performance Evaluation

###  4.1 Random mixed ZX phase circuits

To evaluate the performance of our technique, we test it on random mixed ZX
phase circuits with a varying number 2- and 3-legged gadgets, on 16, 25 and 36
qubits arranged in a square grid.

Our random model samples a fixed number of gadgets num_gadgets on a fixed
number of qubits num_qubits, selecting legs in a independent, identically
distributed way for each gadget. Legs for a gadget are sampled by first
selecting the number of legs uniformly randomly in a given range (from
min_legs to max_legs, both included), and then selecting a uniformly random
subset of qubits with the desired cardinality.

⬇

PhaseCircuit.random(num_qubits, num_gadgets, *, min_legs, max_legs)

Below is an example of a random circuit of 10 gadgets on 6 qubits, with
between 1 and 3 legs:

We can also make our random circuits purely parametric, but this doesn’t
impact our benchmarking.

###  4.2 Internal representation of mixed ZX phase circuits

The internal representation of phase circuits in our technique is optimised
for execution of the the annealing algorithm. A circuit ($n$ qubits, $m_{z}$ Z
gadgets and $m_{x}$ X gadgets, in any order) is broken down into the following
components:

  * •

a pair of binary matrices (an $n\times m_{z}$ matrix $L_{z}$ for Z gadgets and
an $n\times m_{x}$ matrix $L_{x}$ for X gadgets) encoding the positions of the
legs;

  * •

a pair of lists (length $m_{z}$ for Z gadgets and length $m_{x}$ for X
gadgets) mapping the columns of each matrix to the position of the
corresponding gadget in the circuit;

  * •

a list (length $m_{z}+m_{x}$) of angles for the gadgets (concrete or
parametric).

The annealing algorithm operates on gadget legs only—i.e. on the the two
matrices—without any alteration to the original gadget order and angles. The
algorithm further acts on the columns of each matrix independently, leaving
ample scope for caching (implemented) and GPU parallelisation (not yet
implemented). As an example, consider the following circuit of $m_{z}+m_{x}=5$
gadgets on $n=3$ qubits:

The matrices for this circuit are as follows:

$L_{z}=\left(\begin{array}[]{cc}1&1\\\ 1&0\\\ 0&1\end{array}\right)\hskip
85.35826ptL_{x}=\left(\begin{array}[]{ccc}1&0&1\\\ 0&1&1\\\
1&1&0\end{array}\right)$

The lists mapping columns to gadgets are $(0,3)$ for $L_{z}$ and $(1,2,4)$ for
$L_{x}$. The list of angles for the gadgets is
$(\pi/2,\pi,7\pi/4,\pi/4,\pi/2)$.

###  4.3 Benchmarking Results

In all benchmarks below, we proceed as follows:

  1. 1.

We generate 50 random mixed ZX phase circuits, for various combinations of
hyper-parameters—number of qubits, number of gadgets per layer, initial
temperature, schedule type—randomly sampling between 2-legged and 3-legged
gadgets. The number of CX block layers is fixed to 3.

  2. 2.

We run annealing optimisation with varying number of iterations, from 100 to
5000, in increments of 100. The number of circuit layer repetitions is fixed
to 5.

  3. 3.

We look at the reduction in NN CX count—in average or in distribution—from
each original random circuit to the corresponding optimised circuit, relative
to the NN CX count of the original circuit. The NN CX count is computed using
Prim’s algorithm for minimum spanning trees.

To start with, we look at the performance for four different initial
temperatures $t_{0}=1,5,10,20$ and two annealing schedules (linear and
geometric). Below is the average CX count reduction over all runs for the
linear annealing schedule at the four initial temperatures.

Below is the average CX count reduction over all runs for the geometric
annealing schedule at the four initial temperatures. Because of the inferior
average performance of other combinations, we restrict our attention to
linearly annealed runs with $t_{0}=10$ for the remainder of this section.

For our optimisation, we always use 5 layer repetitions for the phase circuit
and 3 layers for the conjugating CX blocks. Below is the average CX count
reduction—as a percentage of the initial CX count—for varying number $m$ of
phase gadgets per circuit layer, ranging from 10 to 35 in increments of 5.

The average CX count reduction degrades as the number of gadgets/layer
increases, decreasing by approximately 0.5% per additional gadget on 36
qubits.

This behaviour is roughly linear when the number of gadgets/layer is small
compared to the number of qubits, but approaches a horizontal asymptote as the
number of gadgets/layer increases. The following graphs showcase this
behaviour over smaller grid topologies ($3\times 3$, $3\times 4$, $4\times 4$
and $4\times 5$ respectively).

The optimised flipping of CX gates allows us to avoid recomputing the whole
conjugated circuit every time. We also implement some basic caching techniques
for gadget legs, and use a modified version of Prim’s algorithm to efficiently
compute the topologically-aware CX count for each gadget. Our Python
implementation is reasonably efficient, taking between $400\mu s$ and $800\mu
s$ per iteration.

Over the full range of iterations we considered, from 100 to 5000, our
experiments never took more than $4s$, with the vast majority under $3s$.

To evaluate the performance of our approach under strict time constraints, we
restrict our attention to the runs that finished under $1s$, using almost all
the time they were allowed (in the 90th to 100th percentile within runs under
$1s$ for the same number of qubits and gadget/layers). We observe that the
average CX count reduction decreases as the number of gadgets/layer increases,
and that the distribution narrows.

To further evaluate the effect of time on our optimisation, we restrict our
attention to runs at 30 gadgets/layer, partitioned into time bins with $0.5s$
resolution. The average CX count reduction improves only slightly as more time
is allowed, indicating diminishing returns for progressively higher number of
iterations. However, we observe a marked improvement of worst-case CX count
reduction from anneals under $0.5s$ to anneals up to $1s$, with only marginal
improvement after that point.

##  5 Future Work

Our first avenue for future work concerns the annealing itself: rather than
using simulated annealing on classical hardware, we have started working on a
formulation of the problem suitable for execution on D-Wave quantum annealers.
This would allow our exploration of the space of configurations to efficiently
scale to much larger examples, involving hundreds of gadgets and qubits. It
would also make our technique a “quantum optimisation of quantum circuits”,
which is fun to say. Thanks to QPL 2022 Anonymous Reviewer #1 for pointing us
to [13], which provides an Ising formulation for the minimum spanning tree
objective.

Our second avenue for future work concerns the scope and method of
optimisation: instead of limiting ourselves to mixed phase gadgets and random
CX circuits, we have have generalised our observation to arbitrary layers of
Pauli gadgets using random Clifford circuits. We are in the process of
implementing this technique and will report on its performance in future work.
We also plan to incorporate commutation between phase gadgets in different
bases to further improve performance.

We will benchmark our techniques explicitly against previous literature, using
realistic circuit families from quantum chemistry, adiabatic quantum
computation, quantum approximate optimisation and quantum machine learning.

## Acknowledgements

The authors would like to acknowledge financial support by Hashberg Ltd for
the execution of numerical experiments. The authors would like to thank the
three anonymous QPL 2022 reviewers for their very helpful comments and
suggestions.

## References

  * [1]
  * [2] Fernando Casas, Ander Murua & Mladen Nadinic (2012): _Efficient computation of the Zassenhaus formula_.  Computer Physics Communications 183(11), pp. 2386–2391, 10.1016/j.cpc.2012.06.006.
  * [3] Alexander Cowtan, Silas Dilkes, Ross Duncan, Will Simmons & Seyon Sivarajah (2020): _Phase gadget synthesis for shallow circuits_.  EPTCS 318, 10.4204/EPTCS.318.13.
  * [4] Alexander Cowtan, Will Simmons & Ross Duncan (2020): _A Generic Compilation Strategy for the Unitary Coupled Cluster Ansatz_.  arXiv preprint, 10.48550/arXiv.2007.10515.
  * [5] Edward Farhi, Jeffrey Goldstone & Sam Gutmann (2014): _A quantum approximate optimization algorithm_.  arXiv preprint, 10.48550/arXiv.1411.4028.
  * [6] Edward Farhi, Jeffrey Goldstone, Sam Gutmann & Michael Sipser (2000): _Quantum computation by adiabatic evolution_.  arXiv preprint quant-ph/0001106, 10.48550/arXiv.quant-ph/0001106.
  * [7] Stefano Gogioso & Richie Yeung (2021): _pauliopt Library_.  Available at https://github.com/sg495/pauliopt.
  * [8] Arianne Meijer-van de Griend & Ross Duncan (2020): _Architecture-aware synthesis of phase polynomials for NISQ devices_.  arXiv preprint, 10.48550/arXiv.2004.06052.
  * [9] S. Kirkpatrick, C. D. Gelatt & M. P. Vecchi (1983): _Optimization by Simulated Annealing_.  Science 220(4598), pp. 671–680, 10.1126/science.220.4598.671.
  * [10] Aleks Kissinger & John van de Wetering (2019): _pyzx library_.  Available at https://github.com/Quantomatic/pyzx.
  * [11] Aleks Kissinger & John van de Wetering (2020): _Pyzx: Large scale automated diagrammatic reasoning_.  EPTCS 318, 10.4204/EPTCS.318.14.
  * [12] Aleks Kissinger & John van de Wetering (2020): _Reducing the number of non-Clifford gates in quantum circuits_.  Physical Review A 102(2), p. 022406, 10.1103/PhysRevA.102.022406.
  * [13] Andrew Lucas (2014): _Ising formulations of many NP problems_.  Frontiers in Physics 2, 10.3389/fphy.2014.00005.
  * [14] Wilhelm Magnus (1954): _On the exponential solution of differential equations for a linear operator_.  Communications on Pure and Applied Mathematics 7(4), pp. 649–673, 10.1002/cpa.3160070404.
  * [15] Beatrice Nash, Vlad Gheorghiu & Michele Mosca (2020): _Quantum circuit optimizations for NISQ architectures_.  Quantum Science and Technology 5(2), p. 025010, 10.1088/2058-9565/ab79b1.
  * [16] Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington & Ross Duncan (2019): _pytket library_.  Available at https://github.com/CQCL/pytket.
  * [17] Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington & Ross Duncan (2020): _t $|$ket$\rangle$: a retargetable compiler for NISQ devices_.  Quantum Science and Technology 6(1), p. 014003, 10.1088/2058-9565/ab8e92.
  * [18] The IBM Qiskit Team & Contributors (2021): _Qiskit: An Open-source Framework for Quantum Computing_ , 10.5281/zenodo.2573505.
  * [19] Richie Yeung (2020): _Diagrammatic Design and Study of Ansätze for Quantum Machine Learning_.  10.48550/arXiv.2011.11073.


# PHYSICAL PARAMETERS OF THE RELATIVISTIC SHOCK WAVES IN GRBs: THE CASE OF 30
GRBs

SAŠA SIMIĆ  Faculty of Science, Department of Physics, University of
Kragujevac
Radoja Domanovića 12, Kragujevac, Serbia, 34000
<EMAIL_ADDRESS>LUKA Č. POPOVIĆ  Astronomical Observatory, Volgina 7
Belgrade, Serbia, 11000
<EMAIL_ADDRESS>
(18.01.2012 for publication in IJMPD)

###### Abstract

Using the modified internal shock wave model we fit the gamma ray burst (GRB)
light and spectral curves of 30 GRBs observed with BATSE. From the best
fitting we obtain basic parameters of the relativistic shells which are in
good agreement with predictions given earlier. We compare measured GRB
parameters with those obtained from the model and discuss connections between
them in the frame of the physical processes laying behind GRB events.

###### keywords:

Gamma rays: bursts; Gamma rays: theory; Shock waves

##  1 Introduction

The quest on resolving a Gamma Ray Burst (GRB) event consists of finding an
explanation for several parts e.g. spatial distribution of the event,
afterglow, spectral and light curve, collimation, etc. In the past decade the
GRB phenomenon has been thoroughly investigated both observationally see e.g.
[Bloom98]-[Ostlin08] and theoretically see e.g. [Meszaros97]-[Zitouni08].
However, a mystery of GRB phenomena lies in its heart, where the central
engine ejects material with relativistic energies and velocities. Due to the
high optical depth of the expanding material in the first phase of the GRB
event, structural observations of the central engine, that is located near the
core of the progenitor, are not possible. The only information that one can
obtain from the observations, in the first minute of a GRB event, is the
temporal variability of the $\gamma$-ray light curve. This usually shows
strong, but short fluctuations of the energy output with a typical time-scale
of the order of milliseconds to a couple of seconds [Norris96]. The numerical
simulations of [Kobayashi97] and [Nakar02] for the early phase of the
explosion, revealed that $\gamma$-ray light curve pulses replicate temporal
activity of the ’inner engine’. This can give information about the connection
between observations and physical processes occurring in the GRB core.

In order to make conclusion about the physical processes in the GRB core one
should assume the mechanism of GRB origin. The mechanism of pulse creation in
a GRB light curve is proposed to be connected with mutual interactions
(collisions) between faster and slower spreading shells see e.g. [Piran05], in
the so called internal shock model. Here we accepted slightly modified version
of the model with the additional assumption of the non-zero density
environment and also with different treatment of the slower material (shell)
accumulated at some distance from the GRB engine. In a previous paper
[Simic07], where we considered interaction of an incoming faster shell with
dense barrier, we demonstrated that the model is able to reproduce and fit
well the observed GRB light curves.

In this paper we apply the model to a sample of 30 GRB light curve pulses
observed with Burst and Transient Source Experiment (BATSE) in order to
discuss the physical parameters of relativistic shells.

The paper is organized as follow: in §2 we give a brief description of the
model and considered assumptions; in §3 we describe the observational sample
and fitting procedure; in §4 we discuss our results and finally in §5 we
outline our conclusions.

##  2 The model and assumptions

Here we describe the physical scenario and some approximations used in the
model of GRB light and spectral curves. First, we assume that the GRB engine
ejects an amount of relativistically expanding material, that spreads
isotropically from the center of explosion. The material is subsequently
ejected from time to time depending on the central engine activity. Here we
will not consider the nature of GRB progenitors, i.e. for the model it is not
relevant whether a GRB event is originated in the collapse of a massive star
[Woosley93] or in the process of merging of two compact objects – neutron
stars or black holes [Norris00], [Fryer99]. The most important for the model
is the assumption that the ejected material can mutually interact or interact
with the surrounding environment.The ejected material is probably irregular in
nature, with different initial parameters (mass $M_{\rm ej}$ and Lorentz
factor $\Gamma_{0}$).

During the expansion of the ejected material (closer to the GRB core), a slow
moving material is followed by a fast moving one, thus the faster moving
material will overtake the slower material and plunge into it. This
interaction produces the relativistic shells and shock waves that can
accelerate particles to very high speeds. The new formed relativistic shells
are probably with different velocities and could collide mutually producing
the observed GRB light curve pulses. The observed GRB light curves show vast
variety of pulses, ranging from very intense ones to those almost equal to
noise, and from symmetrical to highly asymmetrical ones. To model different
observational light curve pulse profiles, we consider here the modified
standard internal model, assuming that the accumulated matter is able to form
a slow moving barrier. This, allowed us to generate more diversity in light
curve pulse profiles.

In the shell interaction processes, the mass of a newly created shell is
approximately the sum of the two colliding shells and its Lorentz factor is
smaller than it was in the faster shell. It could be additionally reduced by
surrounding medium or further shell interactions, thus producing the
accumulation effect. Depending on the density and width of such created dense
barrier, further collision with it produces wider/thinner or higher/lower
pulses. Also, considering the energy density of expanding shells and barriers,
one can get more symmetrically shaped light curve pulses in case of the shells
with high energy and low energy barriers, and high asymmetrical pulse shapes
in the opposite case. This is the main mechanism that we use to explain high
temporal variability of the observed $\gamma$-ray light curves. From the
scenario described above, it follows that the selection of shell parameters
can have a random distribution in a given interval of values.

Similar as in [Kobayashi97], we consider that the ejected GRB material is
organized in an ultra-relativistic flow of well defined and collimated shells
with random initial energy. In contrary to [Piran05] and [Kobayashi97], where
null density hypothesis was used, we assume here that the surroundings regions
around the GRB central engine consist of at least small number of particles
with densities $n_{0}>1\ \rm cm^{-3}$. This allows us to analyze the density
of the moving shell and to model the density distribution and shape of the
barrier in a specific way (here we use Gaussian function, see Eq. 5 further in
the text). The assumption of $n_{0}\neq 0$ around the central engine seems to
be generally valid. If there are some kind of ejections from the central
engine, one can expect an amount of scattered material distributed to the
surrounding region, as e.g. in the collapsar model where a central star is
Wolf-Rayet type. Such star ejects a huge amount of material during its final
stage, therefore one can expect non zero density in its environment.

To describe the evolution of a relativistic shell, we adopt a phenomenological
model based on the [Huang98] that presents a system of the first order of
differential equations (where the distance $R$, Lorentz factor $\Gamma$ and
mass $m_{s}$ of the shell are included, see [Simic07]):

${\frac{{dR}}{{dt}}}=c\sqrt{\Gamma^{2}-1}{\left[{\Gamma+\sqrt{\Gamma^{2}-1}}\right]},$
(1)

${\frac{{d\Gamma}}{{dm_{s}}}}=-{\frac{{\Gamma^{2}-1}}{{M_{\rm
ej}+2(1-\xi)\Gamma m_{s}+\xi m_{s}}}},$ (2)

${\frac{{dm_{s}}}{{dt}}}=2\pi
nm_{p}(1-\cos\theta){\frac{{R^{2}}}{{\Gamma^{3}}}}\left({3\Gamma{\frac{{dR}}{{dt}}}-2R{\frac{{d\Gamma}}{{dt}}}}\right),$
(3)

where the parameter $\xi$ takes values from 0 in case of adiabatic expansion
to 1 in fully radiative case. $M_{\rm ej}$ and $\theta$ are the initial mass
and collimation angle of the shell, $n$ is the number density and $m_{p}$ is
the proton mass. Eqs. (1) - (3) are derived for an observer reference frame,
and they have to be solved simultaneously, together with the density equation.
The initial values of parameters and variables are highly dependent on
physical properties of the shocks.

The density of the barrier cerated from accumulating decelerated shells
(emitted by the central engine) could be described with (see [Blandford76]):

$\frac{n_{2}}{n_{1}}=\frac{\kappa_{2}\gamma_{2}+1}{\kappa_{2}-1}$ (4)

where $n_{2}$ and $n_{1}$ are number densities after and in front of the
shock, $\kappa_{2}$ is the ratio of the specific heat for the shocked fluid
and $\gamma_{2}$ is the Lorentz factor of the shocked fluid.

This equation gives a connection between density of perturbed and unperturbed
material. In the case of the ultra-relativistic expansion, the ratio of the
specific heats has a constant value of $\kappa=4/3$, then Eq. 4 can be reduced
to $\frac{n_{2}}{n_{1}}=4\gamma_{2}+3$. Also, in the relativistic regime, the
Lorentz factor of the shell $\Gamma$ is directly proportional to the Lorentz
factor of the shocked particles $\gamma_{2}$ [Blandford76].

In case of the collision of relativistic shells, a slower interacting shell
(or barrier in our model) presents a density perturbation in surrounding media
for the incoming faster shell. The barrier (or accumulated material) at a
distance $R_{c}$ from the central engine should have a distribution of the
density along the path of penetration of the faster shell. This can be
included in calculation, and here we assume a Gaussian density distribution
as:

$n=n_{0}\left({{\frac{{R_{0}}}{{R}}}}\right)^{s}(4\Gamma+3)\left({1+a\cdot\exp{\left[{-\left({{\frac{{R-R_{c}}}{{b}}}}\right)^{2}}\right]}}\right)$
(5)

where $a$ and $b$ describe the Gaussian intensity and width of the barrier and
$n_{0}$ is the density of the surrounding region. $R_{0}$ is the initial
position of the faster shell, $\Gamma$ is the Lorentz factor of the shell as
used in Eqs. (1-3). Note here that in Eq. 5 (as it was mentioned above)
$n_{0}>0$, and it is a crucial difference between our and the standard
internal shock (IS) model. In the standard IS model, the density of the
interacting shells is fixed by the mass loss rate from the central engine, the
Lorentz factor of the shells and the distance from the source where the
collisions take place. It is not directly related to the density of the
external medium. As a result, even in the absence of an external medium,
prompt emission will result from shocks taking place in the material ejected
from the source. Indeed there is some evidence of bursts occurring in very low
density environments which have a prompt emission but no detectable afterglow,
but as we mentioned above, one could expect that $n_{0}>0$ (especially close
to the accumulated material) in the central engine surroundings.

Generally, one can expect that the ejected material in form of the shell can
collide with the ISM which has a certain density distribution. In such a
highly relativistic physical system the relative motion of charged particles
of the ISM can generate an intense magnetic field in the reference frame of
the moving fluid. We calculate the magnetic field in a similar way as in
[Huang00], by assuming that the energy of the magnetic field is a certain
fraction, $\xi_{b}$, of the total energy of the relativistic shell. In the
comoving reference frame the magnetic field is calculated as:

$B^{{}^{\prime}}=\sqrt{8\pi\xi_{b}n_{0}\Gamma
m_{p}c^{2}(4\Gamma+3)\left({{\frac{{R_{0}}}{{R}}}}\right)^{s}\left({1+a\cdot
e^{-\left({{\frac{{R-R_{c}}}{{b}}}}\right)^{2}}}\right)},$ (6)

where the variables used in this equation are same as in the Eqs. 1-5.

The emission mechanism of shock waves is mainly based on the synchrotron
radiation, but for higher energy bands additional flux may be gained by the
Inverse Compton (IC) radiation [Piran05]. In the first approximation we
neglect the IC radiation.

We calculate the intensity of the radiation emitted by particles in the
relativistic shell using the formulae given by [Rybicki79]. Then the total
emitted flux can be calculated as e.g. in [Huang00a]. Note here that an
expanding shell contains relativistic electrons and baryons which contribute
to the synchrotron radiation. However, taking into account the difference in
velocities of these constituents, one can neglect the contribution of baryons
to the total emitted flux. Then in the comoving reference frame the total flux
is given as:

$P_{\nu}^{{}^{\prime}}=A\cdot\int_{\gamma_{emin}}^{\gamma_{emax}}\gamma_{e}^{-(p+1)}F(\nu^{{}^{\prime}}/\nu_{syn}^{{}^{\prime}})d\gamma_{e}$
(7)

where $A$, $\gamma_{emin}$ and $\gamma_{emax}$ are:

$A=\frac{\sqrt{3}e^{3}B^{{}^{\prime}}}{m_{e}c^{2}}\frac{m_{s}}{m_{p}};\gamma_{emin}=\xi_{e}\Gamma\frac{(p-2)}{(p-1)}\frac{mp}{me};\gamma_{emax}\rightarrow\infty;$
(8)

and

$F(x)=\int_{x}^{\infty}K_{5/3}(x)dx,$

where, $K_{5/3}$ is the Bessel function of the second order and
$\nu_{syn}^{{}^{\prime}}=3\gamma_{e}^{2}eB^{{}^{\prime}}/4\pi m_{e}c$ is the
critical frequency of the synchrotron radiation.

In Eq. 7 we neglect the effects of the surface curvature of an emitting shell,
since it has a small influence on the pulse shape.

In the case of ultra relativistic shells the cooling time is much shorter then
the dynamical (time of expansion) (see [Molinari07]). This is particulary
interesting in the gamma phase of explosion where shells are interacting with
each other.

With the model described above we are able to simulate collisions of
relativistic shells in the first phase of a GRB event, which produce the peaks
in the light curve [Simic07]. This model has been used to fit the light and
energy curves of the sample of 30 GRBs taken from the BATSE database.

##  3 Model vs. observations

In [Simic07] we demonstrated that the model is able to reproduce (simulate)
the observed light curves of GRBs. Moreover, the model can properly fit the
observed light curves. Here we selected 30 GRBs from the BATSE database and
fit them with the model in order to find the basic parameters of interacting
shells. In this section we describe the selected sample and the fitting
procedure.

###  3.1 The sample

From a large BATSE database (3rd channel, $E=100-300$ keV, for the light
curve) we select a sample of 30 GRBs using following criteria: (i) GRBs have
isolated light curves with the clear peak maximum. For proper application of
the model we avoid the complex pulse profiles. 111The process of pulse
creation is stochastic in nature and may result in a complex pulse shape, e.g.
it is often observed that two or more pulses are superposed; (ii) we avoid
small pulses because of their low temporal resolution; (iii) we include in the
sample as much as possible different GRBs (long and short lasting, with
different profiles, different intensity and different profile asymmetries).

In Table 3.1 we give a list and basic parameters of selected GRBs. The
parameters in the table are (from the first to the last column respectively):
the Full Width at Half Maximum of the intensity of light curve pulse (FWHM),
the time of peak intensity for the observed pulse ($t_{\rm peak}$), the total
duration of the pulse from the beginning to the end of its lower tail ($\Delta
t$), the maximal intensity ($J_{m}$) of the pulse measured in ${\rm
erg/cm^{2}sHz}$, and the asymmetry indicator ($w$) calculated as a ratio of
half-halfwidths before and after the maximum.

Parameters of the selected GRB light curves. From the first to the last column
the following parameters are given: Full width at half maximum (FWHM), time of
peak intensity $t_{\rm peak}$, duration of the pulse $\Delta t$, pulse
intensity ($J_{m}\rm[erg/cm^{2}sHz]$), asymmetry indicator $w_{s}$.  GRB FWHM
[s] $\Delta t_{\rm peak}[s]$ $\Delta t[s]$ $J_{m}[\times 10^{-27}]$ $w_{s}$
GRB910629 0.55 0.4 1 13 0.6 GRB911104 0.85 0.5 1.3 13 0.5 GRB920715 0.5 0.45 1
3.2 0.4 GRB920720 0.75 0.6 1.4 40 0.7 GRB920808 0.5 0.25 0.7 11.3 0.4
GRB920811 0.2 0.25 0.3 3.5 0.7 GRB920830 2.75 1.9 6 7.8 0.6 GRB920912 0.4 0.58
1.2 5.4 0.7 GRB920924 0.4 0.4 0.8 4. 0.3 GRB921021 1.8 1 4 5.1 0.4 GRB921207
1.4 1.4 3 60. 0.6 GRB921208 1 0.85 2 2.25 0.8 GRB921222 1 0.85 2 2.25 0.8
GRB950129B 0.6 0.65 1.22 3.7 0.7 GRB950211B 1.5 0.95 3.5 6.5 0.4 GRB960111 1.2
0.5 2.6 9 0.5 GRB960207 0.33 0.25 0.6 9.2 0.8 GRB960229 0.55 0.35 1 5.4 0.8
GRB960311 0.4 0.4 0.7 6 0.6 GRB960409 2.4 2 5 10.5 0.7 GRB960418 0.9 0.8 1.5
4.7 0.4 GRB960524 0.55 0.4 1.2 6.2 0.8 GRB960528 2.2 0.55 4.2 4.2 0.6
GRB960530 5.5 3 11 5.7 0.6 GRB960613 1.6 1.8 4 6 0.3 GRB960617B 0.7 0.55 1.1
3.3 0.8 GRB970424 0.4 0.3 0.7 3.3 0.5 GRB991105 0.65 0.45 1.6 6 0.4 GRB991213
1.4 1 2.4 4 0.4 GRB000107 0.5 0.45 0.7 4.7 1.

Fig. 1 shows statistics of the measured parameters given in Table 3.1. As it
can be seen in Fig. 1 the values of observed parameters do not follow the
Gaussian distribution. Due to the relatively small number of GRBs in the
sample, here we are not able to discuss the power law indices of parameters
(for detailed studies of these GRB properties see [Norris96], [Hakkila08]).

Figure 1: Histograms of basic GRB characteristics for the sample of GRBs
presented in Table 3.1. On the Figures (a) to (d) we place distributions of
(a) the Full Width at Half Maximum of the intensity of light curve pulse
(FWHM), (b) the time of peak intensity for the observed pulse ($t_{\rm
peak}$), (c) the maximal intensity ($J_{m}$) of the pulse and (d) the
asymmetry indicator ($w$)

.

###  3.2 Fitting procedure

In order to follow the relativistic shell evolution and collision, we consider
the case of only one shell expanding from the GRB core. It propagates through
the surrounding media, which can contain a barrier with the mentioned Gaussian
profile. We fit with our model light and energy curves of GRBs given in Table
3.1. In order to find the best fitting we vary the parameters of the faster
shock and barrier. We consider the following parameters as free: the Lorentz
factor $\Gamma_{0}$, the total initial ejected mass of the shell $M_{ej}$, the
density of surrounding media $n_{0}$, the opening angle of the jet
$\theta_{m}$, the distance of collision $R_{c}$ and the parameters $a$ and $b$
which describe the shape (height and width) of the density barrier.
Additionally, we assumed that the barrier can move, thus we also put as a free
parameter the Lorentz factor of the barrier $\Gamma_{b}$. In Fig. 2 (left
panels) we show the best fit of three isolated pulses with different shapes:
GRB000508, GRB911104 and GRB911117. The light curves of these GRBs do not have
a standard form, i.e. the shape of pulses does not always follow the FRED
(Fast Rise Exponential Decay) behavior. As one can see from Fig. 2 the shapes
of the light and energy curves can be very well fitted with the model. In
right panels we show the best fit of the averaged spectral energy distribution
(ASED) taken from all four BATSE channels. Although the data for measured
counts in the energy channels are with large uncertainties, the fit of ASED
can be used for the confirmation of the validity of the GRB light curve fit,
since the same parameters are used to fit both curves.

The coefficients $\xi,\xi_{e}$ and $\xi_{b}$222$\xi$ describes fraction of
total shell energy that has been converted into radiation, $\xi_{e}$ fraction
of the total shell energy devoted to the electron plasma component and
$\xi_{b}$ fraction of total shell energy contained in the magnetic field.,
determine the radiation efficiency for expanding relativistic shell. The
synchrotron radiation is directly proportional to intensity of the magnetic
field, as well as the energy of the electron component of plasma.
Consequently, if one puts high values for those three components it will
increase the intensity of radiation and also the intensity of light curve
pulses. However, the conversion from the kinetic to radiative energy is the
process with relatively low efficiency (see e.g. [Eichler05]), therefore the
above coefficient must not exceed 10 to 20% [Piran05]. They are expected to be
smaller than 0.2. In order to reduce the number of free parameters in the
fitting procedure we fixed them to $\xi=0.1$, $\xi_{e}=0.2$, $\xi_{b}=0.2$
(see also [Zhang07]). Moreover, the different values of these parameters will
mainly affect the intensity of GRB light curve, but not the shape which is
mainly considered in the fitting procedure.

Also, the distribution of electrons in the shell follows a power law function
where the index of the electron distribution $p$ usually takes a value between
2 and 3 [Gallant99]. In the fitting we fix it to $p=2.5$ which well reproduce
the obtained values from the fits of the GRB light and energy curves. In the
fitting procedure we used $\chi^{2}$ minimization.

##  4 Results and discussion

Figure 2: The light (dashed line - left panels) curves and the averaged
spectral energy distribution in four BATSE channels I ch:(20-50)keV, II
ch:(50-100)keV, III ch:(100-300)keV, IV ch:$>$300keV (full circles - right
panels) for GRB911104, GRB911117 and GRB000508 (from top to bottom,
respectively) fitted with the model (solid lines - left panels and dashed
lines right panels).

The parameters of the shells and barriers for 30 GRBs obtained from the best
fittings are given in Table 4. The distribution of parameters are presented in
Fig. 3. As it can be seen from Table 4 (also see Fig. 3) the most of GRB light
curve pulses in our sample indicate an initial Lorentz factor around 95, with
the mean value of about 93. Cases with $\Gamma_{0}<90$ could be explained if
we consider the physical mechanism of the collisions. In the first place, the
initial value of the Lorentz factor directly determines the initial energy and
velocity of the incoming shell. So, for those shells (higher $\Gamma_{0}$)
produced light curve pulses will be more intense and short lasting. On the
contrary, smaller $\Gamma_{0}$ will produce low intensity and long lasting
pulses. This is in a good agreement with the observations of the pulse width -
luminosity correlation found in [Hakkila08].

The internal shell and barrier parameters obtained from the best fitting of
light and energy curves for the sample of 30 GRBs. In the last two rows we put
the mean and appropriate standard deviation of parameters.  parameter
$\Gamma_{0}$ $M_{\rm ej}$ $\Gamma_{b}$ $n_{0}$ $\theta_{m}$ $R_{c}$ $\Delta R$
$n_{b}$ units - $\cdot 10^{-11}[M_{\rm sun}]$ - $[\rm cm^{-3}]$ $[\rm rad]$
$\cdot 10^{14}\rm[cm]$ $\cdot 10^{13}\rm[cm]$ $\cdot 10^{4}\rm[cm^{-3}]$
GRB910629 109 5.5 50 43 0.05 2.0 5.1 43. GRB911104 111 10. 43 63 0.06 2.5 5.1
56.7 GRB920715 103 1.3 51 33 0.06 2.2 5.8 33. GRB920720 90 25. 49 70 0.06 1.5
5.8 105. GRB920808 110 8. 50 33 0.07 1.9 7.8 8.3 GRB920811 111 1.1 65 65 0.06
1.2 2.2 58.5 GRB920830 71 13.5 45 37 0.1 3.5 25. 1.1 GRB920912 100 10. 74 70
0.05 2.0 3.5 4.2 GRB920924 125 2. 50 50 0.04 2.5 3.8 150. GRB921021 95 5. 65
20 0.08 3.5 21.7 0.6 GRB921207 73 53. 63 75 0.1 2.5 10. 150. GRB921208 71 3.8
49 70 0.07 1.9 7.3 4.9 GRB921222 87 2.5 70 20 0.06 2.1 7.8 6. GRB950129B 84 4.
70 20 0.05 1.8 7.5 100. GRB950211B 73 18. 56 90 0.05 1.7 8.3 4.5 GRB960111 95
10. 68 45 0.08 1.9 7.7 1.8 GRB960207 109 6.8 63 57 0.06 1.5 3.2 39.9 GRB960229
98 5.2 61 75 0.057 1.7 4.8 9.8 GRB960311 110 13. 77 19 0.04 1.9 5.2 9.5
GRB960409 75 18. 52 110 0.1 5.5 30. 5.5 GRB960418 75 3. 48 50 0.09 2.7 10.7
50. GRB960524 87 2.7 65 50 0.09 2.1 8.8 335. GRB960528 91 5.5 56 43 0.06 2.0
13.3 2.2 GRB960530 62 40. 46 31 0.055 3.8 30.1 0.6 GRB960613 87 13. 63 65 0.04
2.4 10.2 3.9 GRB960617B 103 3. 87 50 0.07 2.3 8.3 150. GRB970424 99 0.9 53 35
0.061 1.5 3. 105. GRB991105 99 6. 57 69 0.07 2.3 6.8 4.8 GRB991213 75 5. 55 10
0.05 2.1 9.3 9. GRB000107 115 4.3 90 10 0.05 1.5 3.3 10. Mean 93 10 60 50
0.064 2.26 9.4 48 Deviation 13.5 7.6 9.6 19.4 0.014 0.57 5.0 52

Figure 3: Histograms for obtained parameters from the best fit: $\Gamma_{0}$,
$M_{\rm ej}$, $\Gamma_{b}$, $n_{o}$, $\theta_{m}$, $R_{c}$ and $\Delta R$ (a-g
panels, respectively).

The total released isotropic energy of a GRB lies in the range from $8\times
10^{47}$ erg (GRB 980425 associated with supernova SN1998bw, see e.g. Ref
[Galama98], [Pian00]), to $2\times 10^{54}$ erg (GRB 990123 Ref.
[Andersen99]). Note here that using the conical jet model [Frail01] found that
the gamma-ray energies are clustered in a narrow interval around $5\times
10^{50}$ erg. The mass of the ejected shell in the most GRB light curve pulses
is in the interval of $\sim 10^{-11}-10^{-10}M_{\rm sun}$, with a mean value
of $10^{-10}M_{\rm sun}$. These could be used to calculate the initial energy
of a particular shell in the moment of the ejection, giving values in the
range from $10^{45}$ to $10^{47}$ergs. The sum of energies of all ejected
shells during the GRB event can indicate the total energy released by the
central engine in the particular GRB.

Additionally we assume that the barrier is moving, with the Lorentz factor
$\Gamma_{b}$. In the first approximation we take that its velocity is constant
until the collision. This parameter has multiple influence mainly on the
intensity and width of the created light curve pulse, as well as on the FRED
pulse shape. For the faster barriers interaction is long-lasting, but with a
low intensity. On the other hand, if the barrier is moving significantly
slower than the incoming shell, it will cause a more intense interaction
followed by the short-lasting and very strong pulses. The FRED shape is more
dominant in the former case.

The density of the ISM in the region around the central engine is assumed to
be homogeneous $(s=0)$, with a density approximately an order of magnitude
higher than at the distances where afterglow starts (see for example
[Kobayashi04]). We obtained values from several to few tens of particles per
cm-3 (see Table4) that is in the expected range $n_{0}\simeq(10^{0}-10^{3})$
cm-3. The distribution of $n_{0}$ (Fig. 3d) for our sample has a bell-like
shape with the maximal value around 60 cm-3. This values is more appropriate
to hypernova then to a merger scenario, because in the merger scenario one can
expect significantly smaller densities [Hakkila08].

Similar distribution is obtained for opening angle of relativistic shells
$\theta_{m}$, which is considered to be constant during the evolution. The
most probable values are around 0.05 radians ($\approx 3^{o}$) (see [Frail01]
and reference therein). Existence of this particular value is determined by
physical processes in the vicinity of the GRB central engine. Namely, if one
take a higher value of $\theta_{m}$, the resulting pulse is broader and a slow
decay feature is much more visible than in the case of smaller angles, where
the pulse is thinner and has a symmetrical shape. That is in a good agrement
with engaged physical processes during the shell expansion.

We suppose that the shell interaction occurs mainly close to the GRB engine,
(distance $R_{c}\sim 10^{14}\ \rm cm$), as it was proposed in the internal
shock scenario (see [Piran05]). This parameter has very small influence on the
shape of GRB pulses, but has an influence on the intensity of pulses. Pulses
produced in a collision closer to the GRB engine are more intense than ones
originating at larger distances.

The parameters of the Gaussian of the density distribution which describe the
structure of the barrier, the width at the half maximum $b$ and the density in
the central part $a$, both influence the shape of GRB light curve pulses. They
are translated into appropriate variables, number density $n_{b}$ and width
$\Delta R$ of the barrier, respectively, using the connections given by Eq. 9
and 10:

$n_{b}=n_{0}(1+a)$ (9)

$\Delta R=2b\sqrt{\ln\frac{2a}{a-1}}$ (10)

Note here that in the case of dense barrier one can expect $n_{b}/n_{0}>>1$
(see also Table 4), consequently $a>>1$ and Eq. 10 can be rewritten as $\Delta
R=2b\sqrt{\ln{2}}$. The influence of the barrier parameters on pulse profiles
is following: a barrier with narrow width and high number density will have a
strong interaction with fast shell, producing symmetrical and intense light
curve pulses. On the contrary, a barrier with larger width and lower number
density will cause small intensity pulses with high asymmetries.

In general, comparing the obtained values of parameters (Table 4) for
different GRBs, one can conclude that there is no significant difference
between them even when the shapes and durations of GRB pulses are different.
This suggests that the nature of GRBs is similar and that there should be no
big difference between the physical conditions of GRB progenitors (see
[Ghirlanda11]). On the other hand, the barrier density distribution can differ
from the Gaussian, that is assumed here, and it may reflect the values of
basic parameters. But in any case, one can expect that the density
distribution of the barrier has to be taken into account in the shock model.

###  4.1 Connection between the shell parameters

In order to find physical meaning of the obtained parameters, we explore
correlations between them. In Table 4 the parameters are divided in two
groups: one that describes a shell and other connected with the barrier. We
expect that parameters from those two groups are not in correlation, since
they are independent. The results are presented in Figs. 4 and 5, where we
separately denoted the long (t$>$2s) and short (t$<$2s) GRBs with open and
full circles, respectively. As one can see in Figs. 4 and 5 there are no
strong correlations between these quantities, but only in some cases, there is
a slight connections between different parameters, as e.g. $\theta_{m}$ vs.
$\Gamma_{0}$ (Fig.4a), $M_{\rm ej}$ vs. $\Gamma_{0}$ (Fig. 4b). In the case of
the barrier parameters a slight correlation can be found in $\Delta R$ vs.
$n_{b}$ (Fig. 5a). Additionally, we examine the correlation which may be
established between Lorentz factors of the incoming fast shock wave
$\Gamma_{0}$ and the moving barrier $\Gamma_{b}$, presented in Figure 5b.

Figure 4: Parameter dependance: $\theta_{m}$ vs. $\Gamma_{0}$ (a) and $M_{\rm
ej}$ vs. $\Gamma_{0}$ (b). The long (t$>2$s) GRBs are denoted with full circle
and short with open circle. Dashed line presents the border above which there
is no parameter values for this sample of GRBs.

Weak correlations or non-correlations can be noticed in each case particulary.
For example, with larger opening angle of the shock wave $\theta_{m}$, the
volume of the shell increases, that causes smaller initial Lorentz factors of
the shell and vice versa (see Fig. 4a). For the certain shell energy, the
increase of the space angle $\theta_{m}$ causes the increase of the ejected
mass $M_{\rm ej}$, and that causes the decrease of the initial $\Gamma_{0}$.
This is in a agrement with the trends presented in Figure 4b, where one can
see that for a higher value of $M_{\rm ej}$, the initial Lorentz factor tends
to be smaller. Of course, one can not expect to see obvious confirmation for
this conclusion with such small sample of examined GRBs, but rather signs of
trends for mentioned dependence.

Figure 5: Parameter dependance: $\Delta R$ vs. $n_{b}$ (a) and $\Gamma_{0}$
vs. $\Gamma_{b}$ (b), with $r$ designating the correlation coefficient. The
notation for short and long pulses is same as in Fig. 4.

In the case of the barrier parameters, there is an indication of the
connection between the particle number density $n_{b}$ and width of barrier
$\Delta R$ as it is shown in Fig. 5a, as $\Delta R\sim 1/n_{b}$. Physically
one can expect such situation, i.e. for a broader barrier, the density of
barrier tends to be smaller and vice versa.

In Fig. 5b the $\Gamma_{0}$ vs. $\Gamma_{b}$ is presented. There is no global
correlation between these Lorentz factors, but taking into account only short
pulses there is some indication that for a faster barrier the $\Gamma_{0}$ of
the shell is higher. The dashed and solid line in Fig. 5b present the linear
dependance of the given parameters, indicating that the incoming shell must
have higher $\Gamma_{0}$ than the slower moving barrier, providing just enough
necessary conditions for event of collision to happen.

###  4.2 Connection between shell parameters and observed pulse parameters

Additionally, we explore possible correlations between parameters obtained
from fitting the light curves and measured ones (given in §3.1). In Fig. 6 we
plot measured values against $\Gamma_{0}\cdot M_{\rm ej}$, which is
proportional to the energy of the incoming shell.

It is interesting that the pulse intensity for small energies
($\Gamma_{0}\cdot M_{\rm ej}<$0.2) shows nearly linear trend with the energy
(see Fig. 6c), but for $\Gamma_{0}\cdot M_{\rm ej}>0.2$ this trend is not
present.

Figure 6: Parameter dependance: FWHM, $\Delta t$, $J_{m}$ and $w$ vs.
$\Gamma_{0}\cdot M_{\rm ej}$ (scaled to the Max[$\Gamma_{0}\cdot M_{\rm
ej}$]), panels a - d, respectively. The notation for short and long pulses is
same as in Fig. 4 and 5.

In Figs 7a and b we present $\Delta t$ vs. $\Gamma_{0}$, $\Gamma_{b}$,
respectively and in Fig. 7c FWHM vs $\Delta R$. It is obvious that for a
faster barrier and shell the interaction will be shorter. Also, for a broader
barrier we obtain long lasting (wider) pulses as it is shown in Fig. 7c. There
is a correlation between FWHM and $\Delta R$, with correlation coefficients r
= 0.61 for short and r = 0.84 for long lasting GRBs. A linear relationship
between FWHM and $\Delta R$ is present as FWHM=a$\cdot\Delta R$ \+ b (where a
= 0.54 and b = 0.22 for short and a = 1.2 and b = 0.024 for long lasting GRBs,
see Fig. 7c).

Figure 7: Correlation of pulse and shell parameters. Case of $\Delta t$ vs.
$\Gamma_{0}$, $\Gamma_{b}$ (panels a, b), and FWHM vs. $\Delta R$, (panel c).
The notation for short and long pulses is same as in Fig. 4 and 5.

We found some connections between $n_{b}$ vs. $\Delta T$ and $n_{b}$ vs.
$J_{m}$, as presented in Fig. 8. This is expected, since for a barrier with
higher density there are more intense pulses (Fig. 8a). The smaller values of
density produce broader light curve pulses (Fig. 8b). It can be seen from Fig.
8a that short pulses have mainly smaller densities. This can be easily
understood analyzing the physical processes in the moment of collision. When
the particle number density is higher more particles are taking part in the
interaction, producing a higher gamma-ray flux. But if the $n_{b}$ decreases
the interaction is prolonged, caused by the stretching of the barrier
material, so long lasting pulses are produced.

Figure 8: Correlation of pulse and shell parameters. Case of $n_{b}$ on
$\Delta T$ (a) and $n_{b}$ on $J_{m}$ (b).

##  5 Conclusions

In this paper we extracted the basic parameters of internal shock waves during
the first phase of a GRB event by fitting 30 observed GRBs (from BATSE
database). To fit the observed GRB light curves we used the modified internal
shock model given in [Simic07], assuming the collision of a fast shell with a
slow moving (with respect to the velocity of the shell) barrier. We analyzed
the obtained parameters in order to find physical processes behind the GRB
origin, and we came to the following conclusions:

(i) Relativistic shell parameters obtained from the fitting of GRB light
curves are in a good agreement with expected ones and also with estimations
given earlier by other authors mentioned throughout the text.

(ii) The obtained values of internal shell physical parameters for GRBs with
different light curves are in the short interval, showing that the physical
processes behind the GRB creation are similar, i.e. there should be the
ejected mass that collides with surrounding regions - or accumulated slow
moving material.

Also, we analyzed possible connections between parameters obtained from the
best fitting of GRB light curves with measured ones. From this analysis we can
conclude:

(i) There is no correlation between parameters obtained from the best fitting
(Figs 4 and 5), only some indication that long GRBs have higher values of
Lorentz factor, and we found a slight trend between Lorentz factor of the
shell and moving barrier for short pulses.

(ii) There is a correlation between the intensity of pulses and the energy
density of the shell only for a low energy pulses ($\Gamma_{0}\cdot M_{\rm
ej}<0.2$, see Fig. 6c).

(iii) The FWHM of GRB light curve pulses is in the correlation with the width
of the barrier. Using this we give a relation between FWHM (that can be
measured from observed light curves) and $\Delta R$ that is a parameter of the
model (see Fig. 7b).

Finally, we can conclude that the modified internal shock model that assumes
also the barrier density distribution can well describe the first phase of the
GRB origin. Moreover, the obtained parameters for the internal shocks and
barriers well fit the physics of the GRBs.

## Acknowledgments

The work was supported by the Ministry of Education and Science of R. Serbia
through the project ”Astrophysical Spectroscopy of Extragalactic Objects”. We
would like to thank to the anonymous referee for very useful comments and
suggestions and D. Ilić for careful reading of this manuscript and comments.

## References

  * [1] J. S. Bloom, D. A. Frail, S. R. Kulkarni, et al., ApJ 508 (1998) L21.
  * [2] E. Costa, F. Frontera, J. Heise, et al., Nature 387 (1997) 783.
  * [3] G. Ghirlanda, G. Ghisellini, L. Nava, MNRAS in publication.
  * [4] P. J. Groot, T. J. Galama, R. A. M. Vreeswijk, et al., ApJ 502 (1998) L123
  * [5] S. R. Kulkarni, S. G. Djorgoski, A. N. Ramaprakash, et al., Nature 393 (1998) 35.
  * [6] S. R. Kulkarni, S. G. Djorgovski, S. C. Odewahn, et al., Nature 398 (1999) 389.
  * [7] M. I. Andersen, A. J. Castro-Tirado, J. Hjorth, et al., Science 283 (1999) 2075.
  * [8] F. A. Harrison, J. S. Bloom, D. A. Frail, et al., ApJ 523 (1999) L121.
  * [9] T. J. Galama, M. S. Briggs, R. Wijers, et al., Nature 398 (1999) 394.
  * [10] P. A. Price, F. A. Harrison, T. J. Galama et al., ApJ 549L (2001) 7P.
  * [11] S. A. Yost, D. A. Frail, F. A. Harrison, et al., ApJ 577 (2002) 155Y.
  * [12] E. S. Rykoff, D. A. Smith, P. A. Price, et al., ApJ 601 (2004) 101.
  * [13] S. Campana, V. Mangano, A. J. Blustin, P. Brown, D. N. Burrows, G. Chincarini, Nature 442 (2006) 1008
  * [14] G. Östlin, E. Zackrisson, J. Sollerman, S. Mattila, M. Hayes, MNRAS 387 (2008) 1227.
  * [15] P. Mészáros and M. J. Rees, ApJ 476 (1997) 238.
  * [16] M. Vietri, ApJ 488 (1997) L105.
  * [17] E. Waxman, ApJ 485 L5.
  * [18] R. Sari, T. Piran, R. Narayan, ApJ 497 (1998) L17.
  * [19] Y. F. Huang, Z. G. Dai, T. Lu, A&A 336 (1998) L69.
  * [20] Y. F. Huang, Z. G. Dai, T. Lu, MNRAS 309 (1999) 513.
  * [21] E. Nakar and T. Piran, ApJ 598 (2003) 400N.
  * [22] T. Piran, Rev. Modern Phys. 76 (2005) 1143.
  * [23] H. Zitouni, F. Daigne, R. Mochkovich, T. H. Zerguini, MNRAS 386 (2008) 1597.
  * [24] J. P. Norris, R.J., Bonnell, J. T. Nemiroff, AIPC 384 (1996a) 77N.
  * [25] S. Kobayashi, T. Piran and R. Sari, ApJ 490 (1997) 92.
  * [26] E. Nakar, T. Piran, ApJ 572L (2002) 139N.
  * [27] S. Simić, L. Č. Popović, M. I. Andersen, and L. Christensen, A&SS 309 (2007) 173.
  * [28] S. E. Woosley, ApJ 405 (1993) 273.
  * [29] J. P. Norris, G. F. Marani, J. Bonnell, ApJ 534 (2000) 248.
  * [30] C. L. Fryer, S. E. Woosley, M. Herant, M. B. Davies, ApJ 520 (1999) 650.
  * [31] R. D. Blandford and C. F. McKee, Physics of Fluids 19 (1976) 1130.
  * [32] Y. F. Huang, Z. G. Dai, T. Lu, it A&A 355L (2000a) 43H.
  * [33] G. B. Rybicki and A. P. Lightman, Radiative Processes in Astrophysics, (Wiley, New York, 1979).
  * [34] Y. F. Huang, L. J. Gou, Z. G. Dai, T. Lu, ApJ 543 (2000b) 90H.
  * [35] E. Molinari, S. D. Vergani, D. Malesani, et al., A&A 469 (2007) L13.
  * [36] J. Hakkila, T. W. Giblin, J. P. Norris, et al., ApJ 677 (2008) 81.
  * [37] D. Eichler and E. Waxman, ApJ 627 (2005) 861.
  * [38] B. Zhang, E. Liang, K. L. Page, et al., ApJ 655 (2007) 989.
  * [39] Y. A. Gallant, A. Achterberg, MNRAS 305 (1999) L6-10.
  * [40] T. J. Galama, P. J. Groot, J. Vanparadijs, Nature 387 (1998) 479G.
  * [41] E. Pian, L. Amati, L. A. Antonelli, ApJ 536 (2000) 778P.
  * [42] D. A. Frail, S. R. Kulkarni, R. Sari, et al., ApJ 562 (2001) 5.
  * [43] S. Kobayashi, P. Meszaros and B. Zhang, ApJ 601 (2004) 13.