format: text

text:
# Nuclear composition of gamma-ray burst fireballs

## 1 Introduction

Cosmological gamma-ray bursts (GRBs) are powerful explosions in distant galaxies. A physical picture of these phenomena has emerged in the past decade (see Mészáros 2002 for a review): the GRBs are generated by compact, dense, and energetic engines, and they are likely related to black hole formations. The typical mass of the central engine is believed to be a few $`M_{}`$, its size is $`10^610^7`$ cm, and its temperature is $`kT=110`$ MeV. The engine ejects a hot outflow (“fireball”) made of free nucleons, $`e^\pm `$ pairs, trapped blackbody radiation, and magnetic fields.

The initial nuclear composition of the fireball and its evolution during expansion turn out crucial for the explosion physics. In particular, neutrons were shown to be an important component in the explosion (Derishev, Kocharovsky, & Kocharovsky 1999) and an extreme case of fireballs with a large neutron excess was studied (Fuller, Pruet, & Abazajian 2000). Particular attention was given to the relative motion of the neutron and ion components of the fireball, which can lead to observable multi-GeV neutrino emission (Bahcall & Mészáros 2000; Mészáros & Rees 2000). The two components can be ejected with substantially different velocities, and this can affect the observed afterglow emission of GRBs (Pruet & Dalal 2002; Bulik, Sikora, & Moderski 2002). In an accompanying paper (Beloborodov 2002, hereafter Paper 2) we show that the presence of neutrons crucially changes the fireball interaction with an external medium and implies a new mechanism of the afterglow production.

The present paper focuses on processes that shape the nuclear composition of GRB fireballs. In § 2 we study the production of neutrons by the central engine and the resulting neutron-to-proton ratio in the ejected material (scale $`Rr_0<10^7`$ cm). In § 3, we calculate nuclear reactions in the expanding fireball and find abundances of survived free nucleons and synthesized helium ($`R=10^710^9`$ cm). Similar nucleosynthesis calculations, with different codes, have been done recently by Lemoine (2002) and Pruet, Guiles, & Fuller (2002). In § 4, we study spallation of helium at later stages when internal motions develop in the fireball ($`R=10^910^{12}`$ cm). The subsequent dynamics of GRB blast waves with the survived neutron component ($`R=10^{15}10^{17}`$ cm) is investigated in Paper 2.

## 2 Neutronization

There are various models for the central engines of GRBs. It may be a young neutron star whose rotational energy is emitted in a magnetized wind, a neutron star merger, or a massive star collapse. The latter two scenarios proceed via the formation of a black hole of mass $`MM_{}`$ and subsequent disk-like accretion of a comparable mass. The baryonic component of the ejected fireball is then picked up from the accretion disk.

The central engines are sufficiently dense for the electron capture reaction. We will calculate at what densities and temperatures this process creates a neutron excess (neutron-to-proton ratio above unity), and then show that GRB engines are likely to satisfy these conditions. We will illustrate with the accretion-type models of GRBs, where the matter density can be relatively low and neutronization is most questionable.

### 2.1 The equilibrium $`Y_e`$

Consider a dense, $`\rho >10^7`$ g cm-3, and hot, $`kT>m_ec^2`$, matter. The rates of photon emission and absorption are huge and the matter is filled with Planckian radiation. Also, the rates of $`e^\pm `$ pair creation and annihilation ($`\gamma +\gamma e^{}+e^+`$) are huge, and the pairs are in perfect thermodynamic equilibrium with the baryonic matter and radiation. The $`e^\pm `$ number densities are

$$n_\pm =\frac{(m_ec)^3}{\pi ^2\mathrm{}^3}_0^{\mathrm{}}f_\pm (\sqrt{p^2+1})p^2dp.$$
(1)
Here $`p`$ is particle momentum in units of $`m_ec`$ and $`f_\pm `$ is the Fermi-Dirac occupation function,

$$f_\pm (x)=\frac{1}{\mathrm{exp}[(x\mu _\pm )/\theta ]+1}=\frac{1}{\mathrm{exp}[(x\pm \mu )/\theta ]+1},$$
(2)
where $`\theta =kT/m_ec^2`$ and $`\mu \mu _{}=\mu _+`$ is the electron chemical potential in units of $`m_ec^2`$ (thermodynamic equilibrium of $`e^\pm `$ with radiation implies $`\mu _++\mu _{}=0`$). In addition to equation (1), we have the charge-neutrality condition,

$$n_{}n_+=Y_e\frac{\rho }{m_p},$$
(3)
where $`Y_e`$ is the proton-to-nucleon ratio, which would equal the electron-to-nucleon ratio in the absence of $`e^\pm `$ pairs. Equations (1) and (3) determine $`\mu `$ and $`n_\pm `$ for given $`T`$, $`\rho `$, and $`Y_e`$. The electrons become degenerate if $`\mu `$ exceeds $`\theta `$, which happens below the characteristic degeneracy temperature,

$$\theta _{\mathrm{deg}}=\frac{\mathrm{}}{m_ec}\left(\frac{\rho }{m_p}\right)^{1/3},kT_{\mathrm{deg}}=7.7\left(\frac{\rho }{10^{11}\mathrm{g}\mathrm{cm}^3}\right)^{1/3}\mathrm{MeV}.$$
(4)
Degeneracy exponentially suppresses the positron density, $`n_+/n_{}\mathrm{exp}(\mu /\theta )`$, because then $`e^\pm `$ are created only at energies $`E/m_ec^2\mu >\theta `$, in the exponential tail of the thermal distribution.

At temperatures and densities under consideration, the baryonic matter is in nuclear statistical equilibrium, and it is dominated by free nucleons in the unshadowed region of Figure 1. The boundary of this region has been calculated with the Lattimer-Swesty Equation of State code (Lattimer & Swesty 1991). The free protons and neutrons can capture $`e^{}`$ and $`e^+`$ via charged current reactions,

$$e^{}+pn+\nu ,e^++np+\overline{\nu }.$$
(5)
These reactions can rapidly convert protons into neutrons and neutrons back into protons, and establish an equilibrium $`Y_e=n_p/(n_n+n_p)`$, where $`n_p`$ and $`n_n`$ are number densities of protons and neutrons, respectively.111The neutron decay $`np+e^{}+\overline{\nu }`$ is a slow process on GRB timescales and it is neglected here. We now calculate the equilibrium $`Y_e(T,\rho )`$, and in § 2.2 we will show how it applies to GRB central engines.

The exact equilibrium $`Y_e`$ depends on whether the opposite reactions — reabsorption of the emitted $`\nu `$ and $`\overline{\nu }`$ — are also significant. We first consider the $`\nu `$-transparent case, where reabsorption can be neglected, and then address the $`\nu `$-opaque case.

#### 2.1.1 Neutrino-transparent matter

The rates of $`e^{}`$ and $`e^+`$ capture can be derived from the standard electro-weak theory (e.g., Shapiro & Teukolsky 1983, Bruenn 1985). We assume the nucleons to be non-degenerate, and then the rates are

$$\dot{n}_{e^{}p}=Kn_p_0^{\mathrm{}}f_{}(\omega +Q)(\omega +Q)^2\left[1\frac{1}{(\omega +Q)^2}\right]^{1/2}\omega ^2d\omega ,$$
(6)
$$\dot{n}_{e^+n}=Kn_n_{Q+1}^{\mathrm{}}f_+(\omega Q)(\omega Q)^2\left[1\frac{1}{(\omega Q)^2}\right]^{1/2}\omega ^2d\omega ,$$
(7)
where $`\omega `$ is neutrino energy in units of $`m_ec^2`$, $`Q=(m_nm_p)/m_e=2.531`$, and $`K6.5\times 10^4`$ s-1. The constant $`K`$ can be expressed in terms of the mean lifetime of neutrons with respect to $`\beta `$-decay, $`\tau _\beta 900`$ s, as $`K(1.7\tau _\beta )^1`$ (Shapiro & Teukolsky 1983, pp. 316, 524).

An equilibrium $`Y_e`$ is established when the rates of $`e^{}`$ and $`e^+`$ captures are equal,222In the $`\nu `$-transparent regime, neutrinos are sometimes prescribed a zero chemical potential, and the balance $`\mu +\mu _p=\mu _n+\mu _\nu `$ is used with $`\mu _\nu =0`$ to determine the equilibrium $`Y_e`$. In fact, the balance of chemical potentials does not hold because the neutrinos are out of thermodynamic equilibrium. This balance would be valid only in the cold limit $`T0`$ (Landau & Lifshitz 1980).

$$\dot{n}_{e^+n}=\dot{n}_{e^{}p}.$$
(8)
Equations (1), (3), and (8) determine $`Y_e`$ for given $`\rho `$ and $`T`$. Contours of the function $`Y_e(T,\rho )`$ on the $`T\rho `$ plane are shown in Figure 1.

The transition from a proton excess ($`Y_e>0.5`$) to a neutron excess ($`Y_e<0.5`$) takes place in the region of mild degeneracy where $`\mu <\theta `$. We now focus on this region and derive the equilibrium $`Y_e`$ analytically. At $`\mu <\theta `$ and $`\theta >Q+1`$ equations (6) and (7) simplify,

$$\dot{n}_{e^{}p}=Kn_p\theta ^5\left[\frac{45}{2}\zeta (5)+\frac{7\pi ^4}{60}\frac{(2\mu Q)}{\theta }\right],$$
(9)
$$\dot{n}_{e^+n}=Kn_n\theta ^5\left[\frac{45}{2}\zeta (5)\frac{7\pi ^4}{60}\frac{(2\mu Q)}{\theta }\right],$$
(10)
where $`\zeta (5)=1.037`$ is Riemann $`\zeta `$-function. We neglected here next-order terms $`O(Q^2/\theta ^2)`$, $`O(\mu ^2/\theta ^2)`$, and $`O[(Q+1)^5/\theta ^5]`$, and used the formula $`_0^{\mathrm{}}(e^x+1)^1x^ndx=(12^n)\mathrm{\Gamma }(n+1)\zeta (n+1)`$ with $`\mathrm{\Gamma }(n+1)=n!`$ for integer $`n`$. Equating the two rates, we get the equilibrium $`Y_e=n_p/(n_n+n_p)`$,

$$Y_e=\frac{1}{2}+\frac{7\pi ^4}{1350\zeta (5)}\frac{(Q/2\mu )}{\theta }=\frac{1}{2}+0.487\frac{(Q/2\mu )}{\theta }.$$
(11)
In the non-degenerate limit, $`\mu /\theta 0`$, this gives $`Y_e=0.5+0.616/\theta >0.5`$ and implies a proton excess, which is due to the positive difference $`Q`$ between the neutron and proton masses. A very mild degeneracy $`\mu =Q/2<\theta `$ is sufficient to drive $`Y_e`$ below 0.5. This happens because the $`e^+`$ density is reduced by the degeneracy effects and the $`e^{}`$ capture becomes preferential.

It is instructive to write the $`e^{}`$ and $`e^+`$ densities using the linear expansion of equation (1) in $`\mu /\theta `$,

$$n_\pm =\frac{1}{\pi ^2\lambda \text{}^3}\left[\frac{3}{2}\zeta (3)\theta ^3\frac{\pi ^2}{6}\mu \theta ^2\right],\mu <\theta ,$$
(12)
where $`\lambda \text{}=\mathrm{}/m_ec=3.862\times 10^{11}`$ cm and $`\zeta (3)=1.202`$. Then equation (3) yields

$$\mu =3Y_e\frac{\lambda \text{}^3\rho }{m_p\theta ^2}.$$
(13)
The condition $`\mu >Q/2`$ that defines the neutron-excess region on the $`T\rho `$ plane ($`Y_e<0.5`$) can now be written as $`\theta <\theta _n(\rho )`$, where

$$\theta _n=\left(\frac{3\lambda \text{}^3\rho }{Qm_p}\right)^{1/2},kT_n=33\rho _{11}^{1/2}\mathrm{MeV}.$$
(14)
This simple formula perfectly coincides with the numerical results (Fig. 1). That a very mild electron degeneracy is sufficient to create a neutron excess can also be seen by comparing $`T_n`$ with $`T_{\mathrm{deg}}`$,

$$\frac{T_n}{T_{\mathrm{deg}}}4.3\rho _{11}^{1/6}.$$
(15)
A useful explicit formula for the equilibrium $`Y_e(T,\rho )`$ in $`\nu `$-transparent matter with mild degeneracy is derived from equations (11) and (13),

$$Y_e(\theta ,\rho )=\frac{1}{2}\left(1+\frac{0.487Q}{\theta }\right)\left(1+1.46\frac{\lambda \text{}^3\rho }{m_p\theta ^3}\right)^1.$$
(16)
It agrees with the numerical calculations shown in Figure 1 with a high accuracy, $`\delta Y_e/Y_e<1`$% at $`Y_e>0.35`$. In a more degenerate region, where the equilibrium $`Y_e<0.35`$, the formula can still be used, though its error increases to 10% at $`Y_e=0.2`$ and 30% at $`Y_e=0.1`$.

#### 2.1.2 Neutrino-opaque matter

If the matter is opaque to the emitted neutrinos, a complete thermodynamic equilibrium is established. A detailed balance now holds, $`e^{}+pn+\nu `$ and $`e^++np+\overline{\nu }`$, and the equilibrium $`Y_e`$ is determined by the condition

$$\mu +\mu _p=\mu _n+\mu _\nu ,$$
(17)
where $`\mu `$, $`\mu _p`$, $`\mu _n`$, and $`\mu _\nu `$ are chemical potentials (in units of $`m_ec^2`$) of the electrons, protons, neutrons, and neutrinos, respectively, all including the particle rest-mass energy. The antineutrinos have chemical potential $`\mu _{\overline{\nu }}=\mu _\nu `$, so that $`\mu _++\mu _n=\mu _p+\mu _{\overline{\nu }}`$ is also satisfied. The neutrons and protons have Maxwellian distributions, which gives $`n_nn_p^1=\mathrm{exp}[(\mu _nm_ec^2m_nc^2)/kT]\mathrm{exp}[(\mu _pm_ec^2+m_pc^2)/kT]`$ and

$$\mu _n\mu _p=\theta \mathrm{ln}(n_n/n_p)+Q.$$
(18)
The thermalized $`\nu `$ and $`\overline{\nu }`$ obey Fermi-Dirac statistics, and they are described in the same way as $`e^\pm `$ (see eqs.  and ); the only difference is that the statistical weight of energy states is one for neutrinos and two for electrons. The neutrino chemical potential vanishes if $`\nu `$ and $`\overline{\nu }`$ have equal densities, $`n_\nu =n_{\overline{\nu }}`$. If, however, the matter emits non-equal numbers of $`\nu `$ and $`\overline{\nu }`$ (its $`Y_e`$ is changing) then $`n_\nu n_{\overline{\nu }}`$ and $`\mu _\nu 0`$.

Suppose $`N_\nu `$ neutrinos and $`N_{\overline{\nu }}`$ anti-neutrinos have been emitted per nucleon. This causes a change of $`Y_e`$,

$$Y_eY_e^0=N_{\overline{\nu }}N_\nu ,$$
(19)
where $`Y_e^0`$ is an initial value that the matter had before the neutrino emission. If all the emitted neutrinos remain trapped in the matter then $`(n_\nu n_{\overline{\nu }})/n_b=Y_eY_e^0`$, where $`n_b=n_n+n_p`$ is the total nucleon density. If a fraction $`x`$ of the emitted neutrino diffused out of the matter then $`(n_\nu n_{\overline{\nu }})/n_b=(1x)(Y_eY_e^0)`$ and

$$\frac{n_\nu n_{\overline{\nu }}}{n_{}n_+}=(1x)\left(1\frac{Y_e^0}{Y_e}\right).$$
(20)
Let us first consider the case $`x1`$ (efficient neutrino cooling). Then $`|n_\nu n_{\overline{\nu }}|n_{}n_+`$ and the neutrino chemical potential can be neglected compared to that of the electrons. Thus, the chemical equilibrium in $`\nu `$-cooled, $`\nu `$-opaque matter reads

$$\mu =\theta \mathrm{log}(n_n/n_p)+Q.$$
(21)
Equation (21) combined with equations (1) and (3) determines an equilibrium $`Y_e`$ as a function of $`T`$ and $`\rho `$. Note however that equation (21) assumes free nucleons, which is invalid in the shadowed region of the $`T\rho `$ plane shown in Figure 1 (the boundary of this region is also plotted in Fig. 2). The chemical balance can be extended to this region if $`\mu _n`$ and $`\mu _p`$ are corrected for the heavy nuclei formation, which we do using the Lattimer-Swesty Equation of State code. Then we find the equilibrium $`Y_e`$ that satisfies $`\mu =\mu _n\mu _p`$. The results are shown in Figure 2.

The neutron excess boundary, $`Y_e=0.5`$, lies in the region where all nucleons are free and the electrons are mildly degenerate, $`\mu <\theta `$. Equation (21) applies here and shows that this boundary is defined by $`\mu =Q`$. Using the linear expansion in $`\mu /\theta `$ (eq. 13), we get a simple equation for $`Y_e(\theta ,\rho )`$,

$$\frac{3Y_e}{\theta ^2}\frac{\rho \lambda \text{}^3}{m_p}=\theta \mathrm{log}\left(\frac{1Y_e}{Y_e}\right)+Q.$$
(22)
In particular, one sees that $`Y_e=0.5`$ corresponds to

$$\theta _n=\left(\frac{3\lambda \text{}^3\rho }{2Qm_p}\right)^{1/2},kT_n=23.1\rho _{11}^{1/2}\mathrm{MeV}.$$
(23)
It coincides exactly with the contour $`Y_e=0.5`$ calculated with the Lattimer-Swesty code (Fig. 2). Note also that $`T_n=3T_{\mathrm{deg}}\rho _{11}^{1/6}`$.

Finally, we address the regime where the neutrinos are not only thermalized but also trapped in the matter ($`x0`$ in eq. ). It can happen in GRB accretion flows with high accretion rates, $`\dot{M}>1M_{}`$/s (e.g. Di Matteo, Perna, & Narayan 2002). Then $`\mu _\nu `$ should not be neglected in the chemical balance. Like equation (13) we derive for neutrinos

$$\mu _\nu =\frac{6\lambda \text{}^3(n_\nu \nu _{\overline{\nu }})}{\theta ^2},$$
(24)
and substitute $`n_\nu \nu _{\overline{\nu }}=(Y_eY_e^0)(\rho /m_p)`$. The chemical balance $`\mu \mu _\nu =\mu _n\mu _p`$ now reads

$$\frac{3(2Y_e^0Y_e)}{\theta ^2}\frac{\rho \lambda \text{}^3}{m_p}=\theta \mathrm{log}\left(\frac{1Y_e}{Y_e}\right)+Q,$$
(25)
and $`Y_e=0.5`$ corresponds to

$$\theta _n=\left[\frac{3(2Y_e^00.5)\lambda \text{}^3\rho }{Qm_p}\right]^{1/2}.$$
(26)
At $`Y_e^0=0.5`$ it coincides with equation (23) as it should: $`Y_e=0.5=Y_e^0`$ requires that $`\nu `$ and $`\overline{\nu }`$ are emitted in equal numbers and $`\mu _\nu =0`$.

#### 2.1.3 Effects of magnetic fields

Magnetic fields have been neglected in the above calculations, which is a valid approximation as long as the field does not affect the particle distribution functions. Magnetic fields can be generated in the GRB central engines by dynamo, and their energy is likely a fraction $`ϵ_B<1`$ of the total energy density $`w`$, so that

$$B=\sqrt{8\pi ϵ_Bw}.$$
(27)
Here $`w`$ includes the energy of baryons, radiation, and $`e^\pm `$; it also includes the neutrino energy if the disk is $`\nu `$-opaque. The field has the strongest effect on light charged particles — $`e^\pm `$. It introduces the discrete energy levels (Landau & Lifshitz 1980),

$$\frac{E_j}{m_ec^2}=\left[1+p_z^2+2j\frac{\mathrm{}\mathrm{\Omega }_B}{m_ec^2}\right]^{1/2},j=0,1,\mathrm{},$$
(28)
where $`\mathrm{}<p_z<\mathrm{}`$ is the component of the electron momentum parallel to the field and $`\mathrm{\Omega }_B=eB/m_ec`$. The magnetic field also changes the phase-space factor in equation (1): $`p^2\mathrm{d}p`$ is replaced by $`(\mathrm{}\mathrm{\Omega }_B/m_ec^2)(\mathrm{d}p_z/2)`$. Both effects are important if $`(\mathrm{}\mathrm{\Omega }_B/m_ec^2)>p_z^2`$. The mean parallel momentum of the relativistic $`e^\pm `$ equals $`\sqrt{3}kT/c`$, and the condition for the field effects to be important reads

$$\frac{\mathrm{}\mathrm{\Omega }_B}{m_ec^2}>3\theta ^2.$$
(29)
For any plausible $`ϵ_B<1`$, the magnetic field is important only where the electrons are degenerate, and the energy density $`w`$ is dominated by either baryons, $`w_b=(3/2)kT\rho /m_p`$, or degenerate electrons, $`w_{\mathrm{deg}}=(9/4)(\pi /3)^{4/3}\mathrm{}c(Y_e\rho /m_p)^{4/3}`$. The condition (29) can be written as

$$\theta <\left(\frac{\rho \lambda \text{}^3}{m_p}\right)^{1/3}\times \{\begin{array}{cc}(4\pi \alpha _fϵ_B/3)^{1/3}\hfill & ϵ_B>(\pi ^3Y_e^4/32\alpha _f),\hfill \\ (\pi /3)^{7/12}(6\alpha _fϵ_B)^{1/4}Y_e^{1/3}\hfill & ϵ_B<(\pi ^3Y_e^4/32\alpha _f).\hfill \end{array}$$
(30)
where $`\alpha _f`$ is the fine structure constant. The upper line corresponds to $`w_b>w_{\mathrm{deg}}`$, and the lower line to $`w_{\mathrm{deg}}>w_b`$.

### 2.2 GRB central engines

There exists a general constraint on the electron degeneracy in the engine. Its derivation makes use of two facts: (1) The engine is gravitationally bound — otherwise it would explode, and a high baryon contamination of the fireball would be inevitable. GRB models normally envision a quasi-static engine that liberates a fraction of its gravitational binding energy and passes it to a tiny amount of mass outside the engine, thus creating a highly relativistic outflow. (2) The engine is compact (size $`r<10^7`$ cm), and it has a relativistic blackbody temperature $`\theta =kT/m_ec^2>1`$.

The sound speed in a gravitationally bound object of size $`r`$ and mass $`M`$ must be smaller than $`(GM/r)^{1/2}`$. It gives the constraint

$$\frac{P}{\rho }<0.1c^2\left(\frac{r}{3r_g}\right)^1,$$
(31)
where $`r_g=2GM/c^2`$ is the gravitational radius of the object. Pressure $`P=P_\gamma +P_\pm +P_\nu +P_b`$ includes contributions from radiation, $`e^\pm `$, neutrinos, and baryons. We have $`P_\gamma +P_\pm =(11/12)aT^4`$ if the $`e^\pm `$ are weakly degenerate, and a maximum $`P_\nu =(7/24)aT^4`$ for each neutrino species if it is thermalized, where $`a=(\pi ^2k^4/15\mathrm{}^3c^3)=7.56\times 10^{15}`$ erg cm-3 K-4 is the radiation constant. Approximately,

$$PaT^4+\frac{\rho }{m_p}kT.$$
(32)
When comparing the two terms in equation (32), it easy to see that

$$\frac{aT^4}{(\rho /m_p)kT}=\frac{\pi ^2}{15}\frac{\theta ^3}{\theta _{\mathrm{deg}}^3}.$$
(33)
Hence, the baryonic pressure gets dominant at $`\theta \theta _{\mathrm{deg}}`$ \[and at even lower temperatures, $`\theta <(\pi /4)Y_e^{4/3}\theta _{\mathrm{deg}}`$, the pressure of degenerate electrons $`P_{\mathrm{deg}}=(3/4)(\pi /3)^{4/3}\mathrm{}c(Y_e\rho /m_p)^{4/3}`$ takes over\]. Our goal is to derive an upper bound on $`\theta /\theta _{\mathrm{deg}}`$ and therefore we consider $`\theta >\theta _{\mathrm{deg}}`$ with $`PaT^4`$. Constraint (31) then reads

$$\frac{\theta }{\theta _{\mathrm{deg}}}<2\rho _{11}^{1/12}\left(\frac{r}{3r_g}\right)^{1/4},$$
(34)
which also gives a lower bound on the electron chemical potential (eq. 13),

$$\frac{\mu }{Q}>5Y_e\rho _{11}^{1/2}\left(\frac{r}{3r_g}\right)^{1/2}.$$
(35)
We know from § 2.1 that the condition $`Y_e<0.5`$ reads $`\mu >Q/2`$ for $`\nu `$-transparent matter and $`\mu >Q`$ for $`\nu `$-opaque matter. We now see that any engine of size $`r`$ and density $`\rho >10^{10}(r/3r_g)^1`$ g/cm3 satisfies this condition and tends to an equilibrium $`Y_e<0.5`$. This is evidently the case in models of neutron star mergers (Ruffert & Janka 1999), as well as magnetized neutron stars. The collapsar scenario (MacFadyen & Woosley 1999) invokes a relatively low-density accretion flow, and here neutronization is questionable. We therefore shall study accretion flows in more detail. We also need to check whether the equilibrium $`Y_e`$ is achieved on the accretion timescale.

#### 2.2.1 GRB accretion flows

All accretion models of GRBs invoke rotation that creates a funnel along which the fireball can escape. The accretion flow is viewed as a rotating disk maintained in hydrostatic balance in the vertical direction, which gradually spirals to the central black hole. A standard model assumes turbulent viscosity in the disk with a stress tensor $`W_{r\varphi }=\alpha c_s^2`$, where $`c_s=(P/\rho )^{1/2}`$ is the isothermal sound speed and $`\alpha =0.010.1`$ (Balbus & Hawley 1998). Same stress tensor can be described in terms of a viscosity coefficient $`\nu =(2/3)\alpha c_sH`$, where $`H`$ is the half-thickness of the disk.

The standard disk theory gives the velocity of accretion at radius $`r`$

$$u^r=\frac{W_{r\varphi }}{\mathrm{\Omega }_\mathrm{K}rS},S(r)=1\left(\frac{r_{\mathrm{in}}}{r}\right)^{1/2},$$
(36)
where $`\mathrm{\Omega }_\mathrm{K}(r)=(GM/r^3)^{1/2}`$ is the angular velocity of Keplerian rotation, $`M`$ is the black hole mass, and $`r_{\mathrm{in}}`$ is the inner boundary of the disk (the marginally stable orbit); $`r_{\mathrm{in}}=3r_g`$ for non-rotating and $`r_{\mathrm{in}}=r_g/2`$ for extremely rotating black holes.333We give here simple Newtonian estimates and trace the hole spin only through its effect on $`r_{\mathrm{in}}`$; other relativistic corrections are modest and weakly affect the results. $`S(r)`$ varies from $`0`$ to $`1`$, and it equals $`0.5`$ at a characteristic radius $`r=4r_{\mathrm{in}}`$ where the integrated dissipation rate peaks ($`4\pi Hr^2\dot{q}^+S/r`$, see eq.  below). One expects rotating black holes in GRBs, and $`r3r_g`$ is a reasonable characteristic radius.

The vertical hydrostatic balance reads $`c_s=H\mathrm{\Omega }_\mathrm{K}`$, which we use to write the accretion time as

$$t_a(r)=\frac{r}{u^r}=\frac{S}{\alpha \mathrm{\Omega }_\mathrm{K}}\left(\frac{H}{r}\right)^2=2.9\times 10^3S\left(\frac{2H}{r}\right)^2\left(\frac{\alpha }{0.1}\right)^1\left(\frac{r}{3r_g}\right)^{3/2}\left(\frac{M}{M_{}}\right)\mathrm{s}.$$
(37)
We neglect here the disk gravity, which is a valid approximation as long as $`\dot{M}t_aM`$ where $`\dot{M}`$ is the accretion rate. The typical $`t_a`$ is much shorter than the burst duration, and accretion is viewed as a quasi-steady process. It can power a relativistic outflow (“fireball”) with luminosity

$$L=ϵ_f\dot{M}c^210^{51}\left(\frac{ϵ_f}{0.01}\right)\left(\frac{\dot{M}}{10^{32}\mathrm{g}\mathrm{s}^1}\right)\mathrm{erg}\mathrm{s}^1,$$
(38)
where $`ϵ_f`$ is the efficiency of $`\dot{M}c^2`$ conversion into a fireball, which is below the net efficiency of accretion ($`ϵ0.1`$). Most of the accretion energy is either carried away by neutrinos or advected by the accretion flow (in case of small neutrino losses).

The disk baryonic density is given by

$$\rho =\frac{\dot{M}t_a}{4\pi r^2H}6.6\times 10^{10}S\dot{M}_{32}\left(\frac{2H}{r}\right)^3\left(\frac{r}{3r_g}\right)^{3/2}\left(\frac{\alpha }{0.1}\right)^1\left(\frac{M}{M_{}}\right)^2\mathrm{g}\mathrm{cm}^3,$$
(39)
and it is heated viscously with rate

$$\dot{q}^+=\frac{3\dot{M}\mathrm{\Omega }_\mathrm{K}^2}{8\pi H}S5.1\times 10^{33}S\dot{M}_{32}\left(\frac{2H}{r}\right)^1\left(\frac{r}{3r_g}\right)^4\left(\frac{M}{M_{}}\right)^3\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1.$$
(40)
The flow has a huge optical depth for radiation and the radiation is trapped — its diffusion is negligible on the accretion timescale. The only cooling mechanism of the flow is neutrino emission which becomes efficient at $`\dot{M}>10^{31}`$ g s-1 and small $`r`$ (e.g. Popham, Woosley, & Fryer 1999; Narayan, Piran, & Kumar 2001; Kohri & Mineshige 2002). An upper bound on the temperature is derived from the assumption that the neutrino cooling is absent. Then the flow does not loose the dissipated energy and instead traps it and advects. Its energy density can be estimated as $`w_{\mathrm{max}}\dot{q}^+t_a(3/8)(r_g/r)\rho c^2`$ (we use eqs.  and  with $`S=0.5`$). The internal energy of such a hot advective flow is dominated by radiation and $`e^\pm `$, so

$$w_{\mathrm{max}}=\frac{11}{4}aT_{\mathrm{max}}^4\frac{3r_g}{8r}\rho c^2,$$
(41)
where factor $`11/4`$ accounts for the contribution of relativistic weakly degenerate $`e^\pm `$ (neutrino contribution can further increase this factor if the neutrinos are reabsorbed, and then $`T_{\mathrm{max}}`$ will be slightly lower). Equation (41) yields

$$kT_{\mathrm{max}}13\rho _{11}^{1/4}\left(\frac{r}{3r_g}\right)^{1/4}\mathrm{MeV}.$$
(42)
The actual temperature can be significantly lower if the neutrino cooling is significant. The main cooling process is $`e^\pm `$ capture on nucleons: $`e^{}+pn+\nu `$ and $`e^++np+\overline{\nu }`$, which also shapes $`Y_e`$ as we discussed in § 2.1. We now evaluate the characteristic accretion rate $`\dot{M}_{\mathrm{eq}}`$ above which the $`e^\pm `$ capture is rapid enough to establish an equilibrium $`Y_e`$.

The equilibrium $`Y_e`$ is achieved when the flow has emitted one neutrino per nucleon. It is easy to see that disks with efficient neutrino cooling always reach the equilibrium. Indeed, the mean energy of the emitted neutrinos, $`\overline{E_\nu }<5kT`$, is below the liberated accretion energy per nucleon, $`E_n=100300`$ MeV, and an efficient cooling implies that more than one neutrino per nucleon is produced. Thus, $`\dot{M}_{\mathrm{eq}}`$ should be looked for in the inefficient (advective) regime with $`TT_{\mathrm{max}}`$. Such a flow is only mildly degenerate, and the rates of $`e^\pm `$ capture read (in zero order in $`\mu /\theta `$)

$`\dot{n}_{e^{}p}1.5\times 10^2n_p\theta ^5\mathrm{cm}^3\mathrm{s}^1,`$ (43)
$`\dot{n}_{e^+n}1.5\times 10^2n_n\theta ^5\mathrm{cm}^3\mathrm{s}^1.`$ (44)
The neutronization timescale is $`t_n=n_p/\dot{n}_{e^{}p}70\theta ^5`$ s, and it should be compared with the accretion timescale $`t_a`$ (eq. 37),

$$\frac{t_n}{t_a}70\frac{\alpha \mathrm{\Omega }_\mathrm{K}}{S\theta ^5}\left(\frac{H}{r}\right)^2.$$
(45)
We substitute $`T=T_{\mathrm{max}}`$, take into account the hydrostatic balance

$$\frac{H}{r}=\left(\frac{w_{\mathrm{max}}}{3\rho }\right)^{1/2}\frac{1}{\mathrm{\Omega }_\mathrm{K}r}=\frac{1}{2},$$
(46)
and use equation (39) with $`S=0.5`$ to get

$`{\displaystyle \frac{t_n}{t_a}}1.7\times 10^2\dot{M}_{32}^{5/4}\left({\displaystyle \frac{r}{3r_g}}\right)^{13/8}\left({\displaystyle \frac{\alpha }{0.1}}\right)^{9/4}\left({\displaystyle \frac{M}{M_{}}}\right)^{3/2}.`$ (47)
We conclude that disks with

$$\dot{M}>\dot{M}_{\mathrm{eq}}=3.8\times 10^{30}\left(\frac{r}{3r_g}\right)^{13/10}\left(\frac{\alpha }{0.1}\right)^{9/5}\left(\frac{M}{M_{}}\right)^{6/5}\mathrm{g}\mathrm{s}^1$$
(48)
have $`t_n<t_a`$ and hence approach the equilibrium $`Y_e`$. Radius $`r`$ entered this expression in the 13/10 power, which implies that the characteristic $`\dot{M}_{\mathrm{eq}}`$ depends significantly on the black hole spin: the characteristic radius (where viscous dissipation peaks) decreases from $`r10r_g`$ to $`rr_g`$ as the hole spin increases from zero to a maximum value.

The equilibrium $`Y_e`$ is below $`0.5`$ if the flow temperature is below $`T_n`$ that was calculated in § 2.1 (eqs.  and ). It is instructive to compare $`T_n`$ with the maximum accretion temperature $`T_{\mathrm{max}}`$ (eq. 42),

$$\frac{T_n}{T_{\mathrm{max}}}=\kappa \rho _{11}^{1/4}\left(\frac{r}{3r_g}\right)^{1/4},$$
(49)
where $`\kappa =2.5`$ for $`\nu `$-transparent and $`\kappa =1.8`$ for $`\nu `$-opaque flows. The condition for a neutron-excess, $`T<T_n`$, is most difficult to satisfy in low-$`\dot{M}`$ flows where the neutrino cooling is inefficient and $`TT_{\mathrm{max}}`$. Such flows are $`\nu `$-transparent, so that $`\kappa =2.5`$ should be used in equation (49). Substituting equation (39), $`T=T_{\mathrm{max}}`$, and $`H/r=1/2`$, we find that $`T_{\mathrm{max}}<T_n`$ if

$$\dot{M}>\dot{M}_n=7.6\times 10^{30}\left(\frac{r}{3r_g}\right)^{1/2}\left(\frac{\alpha }{0.1}\right)\left(\frac{M}{M_{}}\right)^2\mathrm{g}\mathrm{s}^1.$$
(50)
The “neutronization” accretion rate $`\dot{M}_n`$ is comparable to $`\dot{M}_{\mathrm{eq}}`$. Plausible $`\dot{M}`$ in GRB accretion flows are $`10^{32}`$ g/s and higher, and they should have a neutron excess.

### 2.3 De-neutronization in the fireball?

Are outflows from neutron-rich engines also neutron-rich? The fireball picks up baryons from the surface of the central engine, and the surface density is relatively low. $`Y_e`$ might change while the matter is escaping into a fireball.

Let us consider fireballs produced by accretion disks. The disk is turbulent (the turbulence is the source of viscosity that is responsible for accretion), and it is mixed in the vertical direction on the sound-crossing timescale $`t_{\mathrm{mix}}H/c_s`$ (the turbulent velocity is somewhat smaller than $`c_s`$, however, the thickness of surface layers is also smaller than $`H`$, and $`t_{\mathrm{mix}}H/c_s`$ is about right). Given the hydrostatic balance, $`H=c_s/\mathrm{\Omega }_\mathrm{K}`$, one gets

$$t_{\mathrm{mix}}\mathrm{\Omega }_\mathrm{K}^1.$$
(51)
The turbulent material circulates rapidly up to the surface and back to the interior of the disk, and a small portion of it can (also rapidly) escape in each circulation. The escape timescale of the fireball is comparable to $`\mathrm{\Omega }_\mathrm{K}^1`$. As an element of matter elevates to the surface its $`Y_e`$ would increase if it adjusted instantaneously to a new equilibrium value. However, the time $`\mathrm{\Omega }_\mathrm{K}^1`$ the element has before it sinks back into the disk (or escapes) can be too short for the adjustment. Then $`Y_e`$ of the escaping matter corresponds to $`\rho `$ and $`T`$ inside the disk, where it has spent almost all the time before the sudden escape.

To check this picture, let us evaluate the timescale of “de-neutronization” of an initially neutron rich material that has suddenly expanded into a low-density, hot fireball. Neutrons tend to convert back into protons via two charged current reactions: $`e^+`$ capture and $`\nu `$ absorption ($`\beta `$-decay is slow and negligible). The fireball has temperature $`\theta >1`$ and non-degenerate electrons, $`\mu \theta `$. Equation (44) then gives the $`e^+`$ capture timescale

$$t_+=\frac{n_n}{\dot{n}_{e^+n}}\frac{70}{\theta ^5}\mathrm{s},$$
(52)
$$t_+\mathrm{\Omega }_\mathrm{K}=\left(\frac{kT}{8\mathrm{MeV}}\right)^5\left(\frac{r}{3r_g}\right)^{3/2}\left(\frac{M}{M_{}}\right)^1.$$
(53)
For fireball temperatures up to 8 MeV the $`e^+`$-capture is slow compared to $`\mathrm{\Omega }_\mathrm{K}^1`$ and does not affect $`Y_e`$ of the escaping material.

The $`\nu `$ absorption by neutrons is potentially more important for de-neutronization. The cross section of this reaction is

$$\sigma _a(\omega )=2.4\times 10^{44}\left[1f_{}(\omega +Q)\right](\omega +Q)^2\left[1\frac{1}{(\omega +Q)^2}\right]^{1/2}\mathrm{cm}^22.4\times 10^{44}\omega ^2\mathrm{cm}^2,$$
(54)
and the corresponding rate of de-neutronization is

$$\dot{n}_n=cn_n_0^{\mathrm{}}n_\omega \sigma _ad\omega ,$$
(55)
where $`n_\omega =\mathrm{d}n_\nu /\mathrm{d}\omega `$ is the spectrum of the neutrino number density. This gives the timescale of de-neutronization,

$$t_\nu =\frac{n_n}{|\dot{n}_n|}=\frac{1.4\times 10^{33}}{n_\nu \overline{\omega ^2}}\mathrm{s},$$
(56)
where bar signifies an average over a neutrino spectrum. An upper bound on the neutrino density outside the disk is given by

$$n_\nu <\frac{3\dot{M}\mathrm{\Omega }_\mathrm{K}^2S}{8\overline{\omega }m_ec^3},$$
(57)
which states that the energy flux of the (electron) neutrino from the two faces of the disk, $`2F_\nu 2n_\nu \overline{\omega }m_ec^3/\pi `$, is smaller than the accretion energy released per unit time per unit area of the disk, $`F=(3/4\pi )GM\dot{M}S/r^3`$. Thus, we get (with $`S0.5`$)

$$t_\nu \mathrm{\Omega }_\mathrm{K}>10^2\frac{\overline{\omega }}{\overline{\omega ^2}}\dot{M}_{32}^1\left(\frac{M}{M_{}}\right)\left(\frac{r}{3r_g}\right)^{3/2}.$$
(58)
The mean $`\omega `$ and $`\omega ^2`$ are determined by the state of material that emits the neutrinos. If the disk is $`\nu `$-transparent, the neutrinos have the spectrum (see eq. )

$$n_\omega (\omega +Q)^2\left[\mathrm{exp}\left(\frac{\omega +Q\mu }{\theta }\right)+1\right]^1\left[1\frac{1}{(\omega +Q)^2}\right]^{1/2}\omega ^2\omega ^4\left[\mathrm{exp}\left(\frac{\omega }{\theta }\right)+1\right]^1,$$
(59)
where $`\theta `$ is temperature inside the disk, and $`Q`$ and $`\mu `$ have been neglected compared to the typical $`\omega =\overline{\omega }`$. This gives $`\overline{\omega }=(31/6)[\zeta (6)/\zeta (5)]\theta 5.07\theta `$ and $`\overline{\omega ^2}=(63/2)[\zeta (7)/\zeta (5)]\theta ^230.6\theta ^2`$.

The optical depth of the disk for neutrino absorption by neutrons is444Absorption makes a major contribution to the disk opacity for neutrinos.

$$\tau _a\dot{M}_{32}^{3/2}\left(\frac{\omega }{5\theta }\right)^2\left(\frac{T}{T_{\mathrm{max}}}\right)^5\left(\frac{r}{3r_g}\right)^{7/4}\left(\frac{\alpha }{0.1}\right)^{3/2}\left(\frac{M}{M_{}}\right)^2.$$
(60)
Here we assumed $`T>T_{\mathrm{deg}}`$, so that $`Y_e`$ is not much below 0.5 and $`2H/r=(T/T_{\mathrm{max}})^2`$ ($`PT^4`$, see eqs.  and ). From equations (58) and (60) one can see that $`\nu `$-transparent disks ($`\tau _\nu <1`$) have $`t_\nu \mathrm{\Omega }_\mathrm{K}>1`$ for $`\alpha >0.01`$, and hence their fireballs do not experience any significant absorption of neutrinos.

High-$`\dot{M}`$ disks with neutrino luminosity $`L_\nu >10^{53}`$ erg/s create a sufficiently dense bath of neutrinos that can impact $`Y_e`$ of the ejected fireball. High-$`\dot{M}`$ disks are $`\nu `$-opaque and the emitted neutrinos can be approximated by a blackbody spectrum with $`\overline{\omega }(7/2)[\zeta (4)/\zeta (3)]\theta _\nu 3.15\theta _\nu `$ and $`\overline{\omega ^2}15[\zeta (5)/\zeta (3)]\theta _\nu ^212.9\theta _\nu ^2`$, where $`\theta _\nu `$ is the temperature of the $`\nu `$-photosphere in the disk. The disk also emits $`\overline{\nu }`$ from a corresponding $`\overline{\nu }`$-photosphere. Both $`\nu `$ absorption by neutrons and $`\overline{\nu }`$ absorption by protons occur at the base of the fireball and set a new equilibrium $`Y_e`$ such that the rates of $`\nu `$ and $`\overline{\nu }`$ absorptions are equal. One can expect the new $`Y_e`$ to be below $`0.5`$ (see Qian et al. 1993, where a similar problem is discussed in the context of supernova engines, and Pruet, Fuller, & Cardall 2001). This expectation is based on two facts: (1) the number densities of $`\nu `$ and $`\overline{\nu }`$ are approximately equal, $`|n_\nu n_{\overline{\nu }}|n_\nu `$, — otherwise they carry away too large a leptonic number from the disk, and (2) the absorption cross section is proportional to neutrino energy squared (in the limit $`\omega Q`$, which is valid for the hot high-$`\dot{M}`$ disks). It gives the equilibrium $`n_n/n_p\theta _{\overline{\nu }}^2/\theta _\nu ^2`$ where $`\theta _\nu `$ and $`\theta _{\overline{\nu }}`$ are temperatures of the $`\nu `$ and $`\overline{\nu }`$ photospheres. One expects $`\theta _\nu <\theta _{\overline{\nu }}`$ because the $`\nu `$-photosphere is likely closer to the disk surface (the neutronized disk is more opaque for $`\nu `$) and then $`n_n/n_p>1`$. This argument assumes, however, that the disk temperature decreases toward its surface. It may not hold if the dissipation rate peaks at the surface and makes it hotter than the interior of the disk.

## 3 Nucleosynthesis

The ejected fireball has a low baryon loading and a high temperature, and its nucleons are initially in the free state. As expansion proceeds, the fireball cools adiabatically and, when its temperature decreases to $`kT_{}100`$ keV, fusion reactions shape the nuclear composition like they do during the primordial nucleosynthesis in the Universe (Wagoner, Fowler, & Hoyle 1967). In both cases, we deal with an expanding blackbody fireball with initially free nucleons, and the following three parameters control the outcome of nucleosynthesis: (1) photon-to-baryon ratio $`\varphi =n_\gamma /n_b`$, (2) expansion timescale $`\tau _{}`$ at the time of nucleosynthesis, and (3) n/p-ratio prior to the onset of nucleosynthesis. In the Universe, $`\varphi 3\times 10^9`$, $`\tau _{}10^2`$ s, and $`n_n/n_p1/7`$. Below we formulate the nucleosynthesis problem in GRBs and give a qualitative comparison with the big bang. Then we make detailed calculations of nuclear reactions in GRB fireballs.

### 3.1 Fireball model

A customary GRB model envisions a central engine that ejects baryonic matter at a rate $`\dot{M}_b`$ \[g/s\], thermal energy at a rate $`L_{\mathrm{th}}\dot{M}_bc^2`$, and magnetic energy at a rate $`L_\mathrm{P}`$. The produced fireball expands and accelerates as its internal energy is converted into bulk kinetic energy. The fireball is likely to carry strong magnetic fields and the Poynting luminosity $`L_\mathrm{P}`$ may be higher than the thermalized luminosity. In this case, the magnetic fields prolong the acceleration after all the thermal energy has been converted into bulk expansion. A maximum Lorentz factor of the fireball can be estimated as $`\mathrm{\Gamma }_{\mathrm{max}}=(L_{\mathrm{th}}+L_\mathrm{P})/\dot{M}_bc^2`$, and it is at least $`10^2`$ for GRB explosions.

The nucleosynthesis occurs when the fireball temperature drops to about $`kT_{}100`$ keV at a radius $`R_{}10^710^9`$ cm. The timescale of expansion to this radius, $`R_{}/c`$, is shorter than the duration of the engine activity. Therefore, in the nucleosynthesis calculations, the fireball can be modeled as a quasi-steady outflow.

Let the outflow expand in an axisymmetric funnel with a cross section

$$S(R)=S_0\left(\frac{R}{r_0}\right)^\psi .$$
(61)
For example, $`\psi =2`$ for a radial funnel (and also for a spherically symmetric explosion), and $`\psi =1`$ for a parabolic funnel which may develop in a collapsing progenitor of the GRB (Mészáros & Rees 2001). The outflow is a relativistic ideal fluid with baryon density $`\rho `$, pressure $`P`$, and energy density $`w=3P\rho c^2`$; all these magnitudes are measured in the fluid rest frame. In spherical coordinates $`x^i=(t,R,\theta ,\varphi )`$ the outflow has 4-velocity $`u^i=\mathrm{d}x^i/\mathrm{d}\tau =(u^t,u^R,u^\theta ,u^\varphi )`$, where $`\tau `$ is proper time. We assume $`u^\varphi =0`$ and $`Ru^\theta u^t,u^R`$. The latter assumption is satisfied at all radii for a radial explosion ($`u^\theta =0`$) and at $`Rr_0`$ for a collimated explosion ($`\psi <2`$).

The outflow dynamics is governed by the conservation laws $`_i(\rho u^i)=0`$ and $`_i(T_k^i)=0`$, where $`T_k^i=u^iu_k(w+P)+P\delta _k^i`$ is the stress-energy tensor. The electromagnetic tensor is not included here, which greatly simplifies the problem (we are interested in the early hot stage when the expansion is likely driven by the thermal pressure even in the presence of strong fields). Then the baryon and energy conservation laws read

$`S\rho u^R=\dot{M}_b,S(w+P)u^tu^R=L_{\mathrm{th}}.`$ (62)
The high Lorentz factor of the expansion $`\mathrm{\Gamma }=u^t1`$ implies $`u^R/c\mathrm{\Gamma }`$. Equations (62) then yield

$`\rho ={\displaystyle \frac{\dot{M}_b}{S\mathrm{\Gamma }c}},{\displaystyle \frac{w}{\rho }}={\displaystyle \frac{3}{4}}{\displaystyle \frac{L_{\mathrm{th}}}{\dot{M}_b\mathrm{\Gamma }}}.`$ (63)
We assume that the outflow does not exchange mass or energy with the surroundings, i.e., $`\dot{M}_b(R)=\mathrm{const}`$ and $`L_{\mathrm{th}}(R)=\mathrm{const}`$. Equation (63) then gives $`w\rho /\mathrm{\Gamma }`$ while the first law of thermodynamics $`\mathrm{d}(w/\rho )=P\mathrm{d}(1/\rho )`$ gives $`w\rho ^{4/3}`$. Excluding $`\rho `$ from these two relations one gets $`w^{1/4}1/\mathrm{\Gamma }`$ and

$$T\mathrm{\Gamma }=T_0,$$
(64)
where $`T_0`$ is a constant. Equation (64) is strictly valid at temperatures $`kTm_ec^2=511`$ keV where the $`e^\pm `$ energy density can be neglected. At small $`\mathrm{\Gamma }`$ (where $`kT>m_ec^2`$), the inclusion of $`e^\pm `$ reduces $`T`$ by a modest factor $`(11/4)^{1/4}`$.

Equation (62) implies $`w=aT^4(\mathrm{\Gamma }^2S)^1`$ and, combining with equation (64), one gets

$$\mathrm{\Gamma }=\frac{T_0}{T}=\left(\frac{S}{S_0}\right)^{1/2}=\left(\frac{R}{r_0}\right)^{\psi /2}.$$
(65)
The outflow may be collimated already at its base into a solid angle $`\mathrm{\Omega }_0=S_0/r_0^2<4\pi `$ — its transonic dynamics near the central engine is complicated and unknown. With reasonable accuracy, the simple estimate $`L_{\mathrm{th}}S_0aT_0^4c`$ gives

$$kT_01.2\left(\frac{L_{\mathrm{th}}}{10^{51}\mathrm{erg}\mathrm{s}^1}\right)^{1/4}\left(\frac{r_0}{3\times 10^6\mathrm{cm}}\right)^{1/2}\left(\frac{\mathrm{\Omega }_0}{4\pi }\right)^{1/4}\mathrm{MeV}.$$
(66)
The value of $`T_0`$ is most sensitive to the engine size $`r_0`$. It can be as small as $`3\times 10^5`$ cm (if the outflow is powered by a Kerr black hole via the Blandford-Znajek process) or as large as $`10^7`$ cm (if the outflow is powered by an accretion disk).

A major parameter of the nucleosynthesis problem is the ratio of photon density,

$$n_\gamma =\frac{30\zeta (3)w}{\pi ^4kT}\frac{w}{2.70kT},$$
(67)
to baryon density $`n_b`$. This ratio is constant during the expansion,555It is proportional to entropy per baryon: $`s/k=(2\pi ^4/45\zeta [3])n_\gamma /n_b3.602n_\gamma /n_b`$ where $`k`$ is Boltzmann constant, and it must be constant in an adiabatic expansion. and its value is (cf. eqs.  and )

$$\varphi =\frac{n_\gamma }{n_b}=\frac{3}{4}\frac{m_p}{\dot{M}_b}\frac{L_{\mathrm{th}}}{2.7kT_0}=7.7\times 10^4\left(\frac{\eta _{\mathrm{th}}}{300}\right)\left(\frac{kT_0}{1\mathrm{MeV}}\right)^1,\eta _{\mathrm{th}}=\frac{L_{\mathrm{th}}}{\dot{M}_bc^2}.$$
(68)
The typical $`\varphi `$ is 5 orders of magnitude smaller compared to that of the Universe.

The nucleosynthesis must occur during the acceleration stage of the fireball because the acceleration ends at a low temperature $`T=T_0/\mathrm{\Gamma }_{\mathrm{max}}T_{}`$. The accelerated expansion is described in the comoving time by $`\mathrm{d}\tau =\mathrm{d}t/\mathrm{\Gamma }\mathrm{d}R/c\mathrm{\Gamma }`$, which gives (we use eq. )

$$\tau (T)=\{\begin{array}{cc}\frac{r_0}{c}\mathrm{log}\left(\frac{T_0}{T}\right)\hfill & \psi =2,\hfill \\ \frac{r_0}{c}\frac{2}{(2\psi )}\left[\left(\frac{T}{T_0}\right)^{(\psi 2)/\psi }1\right]\hfill & \psi <2.\hfill \end{array}$$
(69)
The time of nucleosynthesis depends on the shape of the funnel: $`\tau (R_{})3r_0/c`$ for $`\psi =2`$ and $`\tau (R_{})2(T_0/T_{})(r_0/c)`$ for $`\psi =1`$. The timescale of the density fall off is 3 times shorter than that of the temperature because $`\rho T^3`$. Therefore a reasonable choice of the characteristic expansion timescale during nucleosynthesis is $`\tau _{}r_0/c`$ for $`\psi =2`$ and $`\tau _{}(2r_0/3c)(T_0/T_{})`$ for $`\psi =1`$.

Note that $`L_{\mathrm{th}}`$, $`\dot{M}_b`$, and $`T_0`$ enter the nucleosynthesis problem only in combinations that determine $`\varphi `$ and $`\tau (T)`$, and play no other role. For example, the initial temperature is not important as long as $`kT_0>200`$ keV.

### 3.2 Simple estimates and comparison to the big bang

First we evaluate the temperature $`T_{}`$ at which we expect the nucleosynthesis to begin. Let us remind why $`kT_{}80`$ keV in the big bang. If the nuclear statistical equilibrium (NSE) were maintained, recombination of nucleons into $`\alpha `$-particles would occur at $`kT200`$ keV (e.g. Meyer 1994). The reason of the nucleosynthesis delay until $`kT80`$ keV is what is sometimes called “deuterium bottleneck”. Before fusing into $`\alpha `$-particles, the nucleons have to form lighter nuclei — deuterium, tritium, or 3He. At $`kT>100`$ keV, these elements have very low equilibrium abundances, which implies a very long timescale for their fusion, and helium is not formed even though it is favored by the NSE. In particular, the deuterium abundance is suppressed by the very fast photodisintegration $`\gamma +dn+p`$ that balances the opposite reaction $`n+pd+\gamma `$ at the equilibrium value

$$Y_d7\times 10^6\varphi ^1Y_nY_pT_9^{3/2}\mathrm{exp}\left(\frac{25.8}{T_9}\right),$$
(70)
(e.g. Esmailzadeh, Starkman, & Dimopoulos 1991). Here $`Y_i=n_i/n_b`$ is abundance of the $`i`$-th element and $`T_9=T/10^9\mathrm{K}=kT/86.17\mathrm{keV}`$. Nucleosynthesis starts when the exponential wins the pre-exponential factor to give a noticeable $`Y_d10^3`$. Substituting $`\varphi =3\times 10^9`$, one finds that it happens at $`kT_{}80`$ keV.

In GRB outflows, $`\varphi `$ is 5 orders of magnitude smaller, and the NSE favors the recombination of nucleons into $`\alpha `$-particles at a temperature as high as $`500`$ keV. For the same reason as in the big bang, the nucleosynthesis is delayed to a lower $`T_{}`$. Using equation (70) with $`\varphi =10^5`$, we find $`kT_{}140`$ keV, which is almost twice as high as the big bang $`T_{}`$ — it depends logarithmically on $`\varphi `$.

At a given $`T`$, nuclear reaction rates scale as $`\dot{Y}\rho \varphi ^1`$ (here $`Y`$ is an element abundance). The ratio of the expansion timescale $`\tau _{}`$ to a reaction timescale, $`\tau _{\mathrm{reac}}=Y/\dot{Y}`$, behaves as $`\tau _{}/\varphi `$. This combination is $`330`$ times smaller in GRBs compared to the big bang, which may seem not a crucial difference. In fact, the big bang nucleosynthesis was dangerously close to a freezeout of the neutron capture reaction $`n+pd+\gamma `$, without which nucleosynthesis cannot start. Therefore the difference by a factor of 10 can be crucial and it is instructive to compare accurately the capture timescale with the expansion timescale in GRBs.

The neutron capture rate varies slowly with temperature (see Fig. 11 in Smith, Kawano, & Malaney 1993), and near 100 keV it is approximately $`\dot{Y}_c2.5\times 10^4\rho Y_pY_n`$ s-1. The reaction timescale is

$$\tau _c=\frac{\mathrm{min}\{Y_p,Y_n\}}{\dot{Y}_c}=\frac{1.2\times 10^9\varphi }{T_9^3\mathrm{max}\{Y_n,Y_p\}}\mathrm{s}.$$
(71)
It should be compared with the timescale of the density fall off, $`\tau _{}`$. For the big bang, $`\tau _{}/\tau _c30`$. For GRBs, we get

$$\tau _c10^4\left(\frac{\varphi }{10^5}\right)\left(\frac{Y_n}{0.5}\right)^1\mathrm{s},$$
(72)
$$\frac{\tau _{}}{\tau _c}\{\begin{array}{cc}\left(\frac{\varphi }{10^5}\right)^1\left(\frac{r_0}{3\times 10^6}\right)\left(\frac{Y_n}{0.5}\right)\hfill & \psi =2,\hfill \\ \left(\frac{\varphi }{10^5}\right)^1\left(\frac{r_0}{3\times 10^6}\right)\left(\frac{Y_n}{0.5}\right)\left(\frac{T_0}{T_{}}\right)\hfill & \psi =1.\hfill \end{array}$$
(73)
One can see that, in a radial explosion ($`\psi =2`$), the neutron capture rate is marginal for a successful nucleosynthesis. In a collimated explosion, the ratio $`\tau _{}/\tau _c`$ is higher by a factor of $`T_0/T_{}10`$, and nucleosynthesis is efficient.

An important difference between the GRBs and the big bang is the n/p-ratio. In the big bang, $`n_n/n_p=1/7`$ ($`Y_e=7/8`$), which leads to 25% mass fraction of helium after the n-p recombination, while 75% of mass remains in protons (and a tiny amount of other nuclei). In GRBs, $`n_n/n_p>`$1 ($`Y_e<0.5`$), and there are leftover neutrons even if all protons are consumed by helium production. The minimum mass fraction of leftover neutrons is

$$X_n=12Y_e.$$
(74)
### 3.3 Detailed calculation

#### 3.3.1 The code

We will keep track of elements with mass numbers less than 5 (like the big bang, the abundances of heavier nuclei are very small). The six elements under consideration are neutrons, protons, 2H, 3H, 3He, and 4He; they are denoted by $`n,p,d,t,3,\alpha `$, respectively, and the photons are denoted by $`\gamma `$. Their abundances are measured by $`Y_i=n_i/n_b`$ or mass fraction $`X_i=A_iY_i`$, where $`n_i`$ and $`A_i`$ are the number density and the mass number of the i-th species. The photon and matter densities are known functions of temperature,

$$n_\gamma =\frac{2}{\pi ^2}\zeta (3)\left(\frac{kT}{\mathrm{}c}\right)^3=2.02\times 10^{28}T_9^3\mathrm{cm}^3,\rho =3.39\times 10^4\varphi ^1T_9^3\mathrm{g}\mathrm{cm}^3.$$
(75)
The evolution of nuclear composition is described by the set of equations

$`\dot{Y_i}`$ $`=`$ $`{\displaystyle Y_kY_l[klij]}{\displaystyle Y_iY_j[ijkl]},`$ (76)
$`\dot{Y}_n`$ $`=`$ $`\dot{Y_d}2\dot{Y_t}\dot{Y_3}2\dot{Y_\alpha }{\displaystyle \frac{Y_n}{\tau _\beta }},`$ (77)
$`\dot{Y}_p`$ $`=`$ $`\dot{Y_d}\dot{Y_t}2\dot{Y_3}2\dot{Y_\alpha }+{\displaystyle \frac{Y_n}{\tau _\beta }}.`$ (78)
Here dot signifies a derivative with respect to proper time $`\tau `$, $`[ijkl]=n_b\overline{\sigma v}_{ijkl}`$ is the rate of reaction $`i+jk+l`$, and all quantities are measured in the rest frame of the outflow. It is sufficient to calculate the reaction rates for deuterium, tritium, and helium isotopes 3He and 4He ($`i=d,t,3,\alpha `$ in eq. 76). Then $`\dot{Y}_n`$ and $`\dot{Y}_p`$ are found from the neutron and proton conservation laws (eqs.  and ). The latter include the $`Y_n/\tau _\beta `$ term — the conversion of neutrons into protons via $`\beta `$-decay with the mean lifetime of neutrons $`\tau _\beta =900`$ s. The $`\beta `$-decay is negligible in GRBs and we keep it for the code tests on big bang nucleosynthesis.

The sums in equation (76) are taken over all possible reactions with participation of the i-th nuclei. Not all reactions are important (Smith et al. 1993). For example, reactions that destroy $`\alpha `$-particles can be neglected as $`Y_\alpha `$ is far below its equilibrium value. We include the following reactions in the calculations: $`n+pd+\gamma `$, $`n+3p+t`$, $`n+3\alpha +\gamma `$, $`p+d3+\gamma `$, $`p+t\alpha +\gamma `$, $`d+dp+t`$, $`d+dn+3`$, $`d+tn+\alpha `$, $`d+3p+\alpha `$, and take their rates from Smith et al. (1993) and Esmailzadeh et al. (1991). We then have

$`\dot{Y_d}`$ $`=`$ $`Y_nY_p[npd\gamma ]Y_dY_\gamma [d\gamma np]Y_pY_d[pd3\gamma ]2Y_d^2[ddn3]Y_dY_3[d3p\alpha ]`$ (79)
$`Y_dY_t[dtn\alpha ]2Y_d^2[ddpt],`$
$`\dot{Y_t}`$ $`=`$ $`Y_d^2[ddpt]+Y_nY_3[n3pt]Y_dY_t[dtn\alpha ]Y_pY_t[ptn3]Y_pY_t[pt\alpha \gamma ],`$ (80)
$`\dot{Y_3}`$ $`=`$ $`Y_pY_d[pd3\gamma ]+Yd^2[ddn3]+Y_pY_t[ptn3]Y_dY_3[d3p\alpha ]Y_nY_3([n3pt]+[n3\alpha \gamma ]),`$ (81)
$`\dot{Y_\alpha }`$ $`=`$ $`Y_dY_t[dtn\alpha ]+Y_dY_3[d3p\alpha ]+Y_nY_3[n3\alpha \gamma ]+Y_pY_t[pt\alpha \gamma ].`$ (82)
The reaction rates are functions of $`T`$ and $`\rho `$, and the set of equations is closed by equation (69) that relates $`T`$ and $`\tau `$. We solve equations (77-82) numerically with initial conditions $`Y_p^0=Y_e`$, $`Y_n^0=1Y_e`$, $`Y_\alpha =0`$, and $`\dot{Y}_d=\dot{Y}_t=\dot{Y}_3=0`$ at an initial temperature $`kT_01`$ MeV. At high temperatures, the abundances of all elements except $`n`$ and $`p`$ are negligibly small. The abundances of deuterium, tritium, and 3He are close to quasi-steady equilibrium at $`kT>150`$ keV, i.e., the rates of their production $`\dot{Y}^+`$ and sink $`\dot{Y}^{}`$ are much higher than the expansion rate $`Y/\tau `$, and hence $`\dot{Y}^+`$ and $`\dot{Y}^{}`$ almost balance each other. Where an element abundance $`Y`$ approaches the equilibrium value, it is calculated by setting the net $`\dot{Y}=\dot{Y}^+\dot{Y}^{}=0`$. (Thus, one avoids numerical integration where $`\dot{Y}`$ is a small difference of big numbers.) The code has been tested with the big bang nucleosynthesis. It successfully reproduced the standard evolution of the cosmological nuclear composition.

#### 3.3.2 Results

Figure 3 compares the nucleosynthesis in GRBs and the big bang (BB). Parameters of the GRB fireball in this example are chosen as $`\varphi =10^5`$, $`r_0=3\times 10^6`$ cm, $`\psi =2`$, and $`Y_e=0.5`$. As expected, GRB nucleosynthesis occurs at higher temperatures. The freezeout mass fractions of helium and deuterium are $`X_\alpha 0.16`$ and $`X_d0.03`$, and about 81% of mass remains in free nucleons. Interestingly, the deuterium evolution is qualitatively different in GRBs: $`X_d`$ increases monotonically and freezes out at the 3% level which is almost 3 orders of magnitude higher than in the big bang (see also Pruet et al. 2002). This difference is due to a high abundance of neutrons and a negligible $`\beta `$-decay. As a result, the deuterium production by the neutron capture reaction outweighes its burning rate and $`X_d`$ grows. Also, the freezeout happens quickly in the radial expansion ($`T\mathrm{exp}[\tau /\tau _0]`$ where $`\tau _0=r_0/c`$) and $`X_d`$ could not decrease much below 1% even if the neutron capture reaction were completely switched off at, e.g., $`kT<80`$ keV.

The post-nucleosynthesis fireball is predominantly composed of free nucleons and helium, and it would be useful to know how their freezeout ratio

$$f=\frac{X_\alpha }{X_n+X_p}$$
(83)
depends on the four parameters of the explosion $`r_0`$, $`\varphi `$, $`Y_e`$, and $`\psi `$. Figure 4 shows $`f(r_0,\varphi )`$ for radial explosions with $`Y_e=0.5`$. Contours $`f=\mathrm{const}`$ on the $`r_0\varphi `$ plane are perfect straight lines $`\varphi r_0^{0.85}`$ (or $`r_0\varphi ^{1.18}`$). This implies that $`f`$ depends on combination $`r_0\varphi ^{1.18}`$. A similar situation takes place for collimated explosions (Fig. 5). The slope of all contours in Figure 5 is $`0.90`$ and $`f`$ can be viewed as a function of $`r_0\varphi ^{1.11}`$. Calculations with $`Y_e0.5`$ give different values of $`f`$, however, contours $`f(r_0,\varphi )=\mathrm{const}`$ have exactly the same slope as in the case of $`Y_e=0.5`$.

The results of nucleosynthesis calculations can be represented in a compact form if we define a new variable

$$\xi =\{\begin{array}{cc}0.2\left(\frac{\varphi }{10^5}\right)^{1.18}\left(\frac{r_0}{3\times 10^6}\right)\hfill & \psi =2,\hfill \\ 7.6\left(\frac{\varphi }{10^5}\right)^{1.11}\left(\frac{r_0}{3\times 10^6}\right)\hfill & \psi =1.\hfill \end{array}$$
(84)
Figure 6 shows the outcome of nucleosynthesis $`f(\xi )`$ for different $`Y_e`$. Strikingly, $`f(\xi )`$ is practically identical for radial and collimated explosions. Figure 6 gives a full description of helium production in relativistic explosions: from this figure and equation (84) one finds $`f`$ for given $`\varphi `$, $`r_0`$, $`Y_e`$, and $`\psi =1,2`$.

One can see that nucleosynthesis is most efficient in $`Y_e=0.5`$ models, where $`X_n^0=X_p^0`$. In the most interesting region of parameters, where $`f`$ changes from $`f1`$ to $`f1`$, one observes a simple behavior

$$f=\xi ,0.03<f<30,$$
(85)
for both radial and collimated explosions (this explains our choice of constants 0.2 and 7.6 in eq. ).

Fireballs with $`Y_e0.5`$ have an upper bound on the freezeout abundance of helium $`\mathrm{max}X_\alpha =2\mathrm{min}\{Y_e,1Y_e\}`$. For explosions with a neutron excess, $`Y_e<0.5`$, this gives a maximum $`f`$,

$$f_{\mathrm{max}}=\frac{2Y_e}{12Y_e}.$$
(86)
With increasing $`\xi `$, the dependence $`f\xi `$ switches to $`f=f_{\mathrm{max}}=\mathrm{const}`$ (Fig. 6).

The post-nucleosynthesis fireball is dominated by $`\alpha `$-particles if $`f>1`$, which requires $`Y_e>0.25`$ and $`\xi >1`$. For example, in a radial explosion with $`Y_e=0.5`$ it requires $`\varphi <2.5\times 10^4(r_0/3\times 10^6)^{0.85}`$, and a similar condition for parabolic explosions reads $`\varphi <6.2\times 10^4(r_0/3\times 10^6)^{0.9}`$. The typical GRB parameters happen to be just marginal for a successful nucleosynthesis: $`f`$ varies from $`10^2`$ to $`f_{\mathrm{max}}`$ in the expected range of $`\varphi `$ and $`r_0`$.

## 4 Spallation

Synthesized helium may be destroyed during the subsequent evolution of the explosion. Spallation reactions can occur when an $`\alpha `$-particle collides with another particle with a relative energy exceeding the nuclear binding energy. The fireball temperature at the acceleration stage is too low for such reactions, however, the collisions may become energetic if (1) there appears a substantial relative bulk velocity between the neutron and ion components or (2) internal shocks occur and heat the ion fireball to a high temperature, much above the blackbody value.

### 4.1 Neutron-ion collisions during the acceleration stage

Near the central engine, the neutron component of the fireball is well coupled to the ion component by elastic collisions with a small relative velocity $`\stackrel{~}{\beta }=\stackrel{~}{v}/c`$ (Derishev et al. 1999, Bahcall & Mészáros 2000). The collision cross section can be approximated as $`\sigma _i=\sigma _0/\stackrel{~}{\beta }`$, where $`\sigma _03\times 10^{26}`$ cm2 if the ions are protons (and $`\sigma _0`$ approaches $`10^{25}`$ cm2 for $`\alpha `$-particles). In the fluid frame, the mean collisional time for neutrons is

$$\tau _{\mathrm{coll}}=\frac{1}{n_i\sigma _i\stackrel{~}{\beta }c}=\frac{1}{n_i\sigma _0c},$$
(87)
where $`n_i`$ is the ion density. The ion density is a fraction of the total nucleon density $`n_b`$ which behaves as $`n_b=\mathrm{const}/\mathrm{\Gamma }^3`$ during the acceleration stage (see eqs.  and ). The constant of proportionality can be expressed in terms of $`\eta _{\mathrm{th}}=L_{\mathrm{th}}/\dot{M}_bc^2`$ and $`\varphi =n_\gamma /n_b`$,

$$n_b=\frac{45^3\zeta ^4(3)}{4\pi ^{14}}\left(\frac{m_pc}{\mathrm{}}\right)^3\frac{\eta _{\mathrm{th}}^3}{\varphi ^4\mathrm{\Gamma }^3}=5.61\times 10^{38}\frac{\eta _{\mathrm{th}}^3}{\varphi ^4\mathrm{\Gamma }^3}\mathrm{cm}^3.$$
(88)
In the fixed lab frame, the ion fluid is accelerated by radiative or magnetic pressure, and its Lorentz factor is doubled on timescale $`t=R/c`$. Neutrons “miss” the ion acceleration by $`\mathrm{\Delta }\mathrm{\Gamma }/\mathrm{\Gamma }t_{\mathrm{coll}}/t=\mathrm{\Gamma }c\tau _{\mathrm{coll}}/R<1`$ and have a smaller Lorentz factor $`\mathrm{\Gamma }_n=\mathrm{\Gamma }\mathrm{\Delta }\mathrm{\Gamma }`$. The relative velocity of the neutron and ion components is $`\stackrel{~}{\beta }=(\mathrm{\Gamma }^2\mathrm{\Gamma }_n^2)/(\mathrm{\Gamma }^2+\mathrm{\Gamma }_n^2)(\mathrm{\Gamma }\mathrm{\Gamma }_n)/\mathrm{\Gamma }`$, which is

$$\stackrel{~}{\beta }=\frac{t_{\mathrm{coll}}}{t}=\frac{\mathrm{\Gamma }}{Rn_i\sigma _0}\frac{\mathrm{\Gamma }^4}{R}.$$
(89)
The energy of neutron-ion collisions becomes sufficient for spallation reactions if $`\stackrel{~}{\beta }`$ exceeds $`\stackrel{~}{\beta }_{\mathrm{sp}}0.1`$. This can happen at late stages of the fireball acceleration, at high $`\mathrm{\Gamma }`$ but not exceeding $`\mathrm{\Gamma }_{\mathrm{max}}=L/\dot{M}_bc^2=\eta _{\mathrm{th}}+\eta _\mathrm{P}`$ ($`L=L_{\mathrm{th}}+L_\mathrm{P}`$ is the total luminosity of the fireball that includes the Poynting flux). In the case of $`\mathrm{\Gamma }(R)R/r_0`$ (radial explosion), $`\stackrel{~}{\beta }`$ reaches $`\stackrel{~}{\beta }_{\mathrm{sp}}`$ when $`\mathrm{\Gamma }`$ reaches

$$\mathrm{\Gamma }_{\mathrm{sp}}=0.3\mathrm{\Gamma }_{\mathrm{max}}\left(\frac{L_{\mathrm{th}}}{L}\right)\left(\frac{\varphi }{10^5}\right)^{4/3}\left(\frac{\stackrel{~}{\beta }_{\mathrm{sp}}}{0.1}\right)^{1/3}\left(\frac{r_0}{3\times 10^6}\right)^{1/3}\left(\frac{4n_i}{n_b}\right)^{1/3}.$$
(90)
Spallation takes place if $`\mathrm{\Gamma }_{\mathrm{sp}}<\mathrm{\Gamma }_{\mathrm{max}}`$, which requires $`\varphi >\varphi _{\mathrm{sp}}`$,

$$\varphi _{\mathrm{sp}}4\times 10^4\left(\frac{L_{\mathrm{th}}}{L}\right)^{3/4}.$$
(91)
We keep here $`L_{\mathrm{th}}/L<1`$ as the most uncertain parameter which may be much below unity. Note also, that in the case of a non-radial explosion, $`r_0`$ should be replaced by $`R/\mathrm{\Gamma }`$ in equation (90). Using equation (68) it is easy to show that the condition $`\varphi >\varphi _{\mathrm{sp}}`$ is equivalent to

$$\mathrm{\Gamma }_{\mathrm{max}}>160\left(\frac{kT_0}{\mathrm{MeV}}\right)\left(\frac{L}{L_{\mathrm{th}}}\right)^{1/4}.$$
(92)
Let $`Y_n`$ be abundance of neutrons that survived the nucleosynthesis, and suppose that $`\varphi >\varphi _{\mathrm{sp}}`$. The lifetime of $`\alpha `$-particles bombarded by neutrons with $`\stackrel{~}{\beta }\stackrel{~}{\beta }_{\mathrm{sp}}`$ is $`t_{\mathrm{life}}(t_{\mathrm{coll}}/Y_n)(\sigma _0/\sigma _{\mathrm{sp}})(\stackrel{~}{\beta }_{\mathrm{sp}}/Y_n)t`$, where $`\sigma _{\mathrm{sp}}\sigma _0`$ is the spallation cross section (it is roughly equal to the geometrical size of the nucleus $`\pi r_N^25.3\times 10^{26}A^{2/3}\mathrm{cm}^2`$ with $`A=4`$ for $`\alpha `$-particles). Hence a modest $`Y_n0.1`$ should be sufficient for a significant spallation.

In case of a very low baryon loading (high $`\mathrm{\Gamma }_{\mathrm{max}}`$) the neutrons decouple from the ions before the end of the acceleration stage and their Lorentz factor $`\mathrm{\Gamma }_n`$ saturates at $`\mathrm{\Gamma }_{\mathrm{dec}}<\mathrm{\Gamma }_{\mathrm{max}}`$. The relative velocity $`\stackrel{~}{\beta }`$ approaches unity when the decoupling happens, and the last neutron-ion collisions are energetic enough for pion production which leads to multi-GeV neutrino emission (see Derishev et al. 1999, Bahcall & Mészáros 2000). The fireball Lorentz factor at this moment, $`\mathrm{\Gamma }_{\mathrm{dec}}`$, is given by equation (90) with $`\stackrel{~}{\beta }_{\mathrm{sp}}`$ replaced by unity, i.e. $`\mathrm{\Gamma }_{\mathrm{dec}}=10^{1/3}\mathrm{\Gamma }_{\mathrm{sp}}`$. The decoupling takes place if $`\mathrm{\Gamma }_{\mathrm{dec}}<\mathrm{\Gamma }_{\mathrm{max}}`$, which requires $`\varphi `$ greater than

$$\varphi _{\mathrm{dec}}7\times 10^4\left(\frac{L_{\mathrm{th}}}{L}\right)^{3/4}.$$
(93)
The condition $`\varphi >\varphi _{\mathrm{dec}}`$ is equivalent to $`\mathrm{\Gamma }_{\mathrm{max}}>300(kT_0/\mathrm{MeV})(L/L_{\mathrm{th}})^{1/4}`$. The decoupling is always preceded by spallation.

The upper bound on the neutron Lorentz factor due to decoupling is essentially determined by the rate of baryon outflow per unit solid angle $`\dot{M}_\mathrm{\Omega }`$ \[g/s\], and it is useful to rewrite $`\mathrm{\Gamma }_{\mathrm{dec}}`$ in terms of $`\dot{M}_\mathrm{\Omega }`$. We substitute $`n_i=n_b=(\dot{M}_\mathrm{\Omega }/R^2\mathrm{\Gamma }m_pc)`$ to equation (87), and then the decoupling condition $`\mathrm{\Gamma }\tau _{\mathrm{coll}}=R/c`$ gives the maximum Lorentz factor of neutrons

$$\mathrm{\Gamma }_{\mathrm{dec}}\left(\frac{\sigma _0\dot{M}_\mathrm{\Omega }}{r_0m_pc}\right)^{1/3}300\left(\frac{\dot{M}_\mathrm{\Omega }}{10^{26}\mathrm{g}\mathrm{s}^1}\right)^{1/3}\left(\frac{r_0}{3\times 10^6\mathrm{cm}}\right)^{1/3}.$$
(94)
### 4.2 Internal shocks

The $`\alpha `$-particles can also be destroyed later on, when internal shocks develop in the fireball. Lorentz factor $`\mathrm{\Gamma }_{\mathrm{max}}`$ fluctuates if the central engine is “noisy” during its operation $`0<t<t_b1`$ s. The fluctuations probably occur on timescales $`r_0/c10^4`$ s and longer, up to $`t_b`$. This leads to internal shocks (Rees & Mészáros 1994). Internal dissipation of the velocity fluctuations may give rise to the observed GRB, and this picture is plausible because it easily accounts for the observed short-timescale variability in GRB light curves. The amplitude of the fluctuations is described by the dimensionless rms of the Lorentz factor, $`A=\delta \mathrm{\Gamma }_{\mathrm{max}}/\mathrm{\Gamma }_{\mathrm{max}}`$. At modest $`A<1`$, the internal dissipation proceeds in a simple hierarchical manner (Beloborodov 2000). The internal collisions begin at

$$R_0=\frac{2\mathrm{\Gamma }_{\mathrm{max}}^2}{A}\lambda _0=5.4\times 10^{11}A^1\left(\frac{\mathrm{\Gamma }_{\mathrm{max}}}{300}\right)^2\left(\frac{\lambda _0}{3\times 10^6}\right)\mathrm{cm},$$
(95)
where $`\lambda _0`$ \[cm\] is a minimum length scale of the fluctuations, which is likely comparable to $`r_0`$. The shortest fluctuations are dissipated first, and at $`R>R_0`$ fluctuations with $`\lambda >\lambda _0`$ are dissipated in the hierarchical order.

#### 4.2.1 $`n`$-$`\alpha `$ collisions

Suppose $`\mathrm{\Gamma }_{\mathrm{max}}`$ is small enough, so that the neutrons decouple from the ions after the acceleration stage and have same $`\mathrm{\Gamma }_n=\mathrm{\Gamma }_{\mathrm{max}}`$. As soon as internal shocks develop in the ion component, the neutrons begin to drift with respect to the ions. Their relative velocity $`\stackrel{~}{\beta }A`$ is sufficient for spallation reactions if $`A>0.1`$ (and for $`A1`$ the pion production and multi-GeV neutrino emission takes place, see Mészáros & Rees 2000).

Let $`L_\mathrm{\Omega }=R^2n_bm_pc^3\mathrm{\Gamma }_{\mathrm{max}}^2`$ be the fireball kinetic luminosity per unit solid angle; then

$$n_b=\frac{L_\mathrm{\Omega }}{m_pc^3R^2\mathrm{\Gamma }_{\mathrm{max}}^2}.$$
(96)
The $`\alpha `$-particles are destructed by neutrons with rate $`\dot{Y}_\alpha Y_\alpha Y_nn_b\sigma _{\mathrm{sp}}c`$ and their lifetime in the fluid frame is $`\tau _{\mathrm{life}}=(Y_nn_b\sigma _{\mathrm{sp}}c)^1`$,

$$\tau _{\mathrm{life}}=\frac{m_pc^2R^2\mathrm{\Gamma }_{\mathrm{max}}^2}{\sigma _{\mathrm{sp}}Y_nL_\mathrm{\Omega }}.$$
(97)
It should be compared with the timescale of side expansion in the fluid frame, $`(R/c\mathrm{\Gamma }_{\mathrm{max}})`$. If $`\tau _{\mathrm{life}}<(R/c\mathrm{\Gamma }_{\mathrm{max}})`$ at $`R=R_0`$ then most of the $`\alpha `$-particles are destroyed by the internal shocks. This condition reads

$$\mathrm{\Gamma }_{\mathrm{max}}<\left[\frac{AY_n\sigma _{\mathrm{sp}}L_\mathrm{\Omega }}{2\lambda _0m_pc^3}\right]^{1/5}300A^{1/5}Y_n^{1/5}\left(\frac{L_\mathrm{\Omega }}{10^{52}}\right)^{1/5}\left(\frac{\lambda _0}{3\times 10^6}\right)^{1/5}.$$
(98)
#### 4.2.2 $`\alpha `$-$`\alpha `$ collisions

The $`\alpha `$-particles acquire random energy $`(A^2/2)4m_pc^2`$ in internal shocks, which easily exceeds the spallation threshold. The condition $`\tau _{\mathrm{life}}<(R/c\mathrm{\Gamma }_{\mathrm{max}})`$ for efficient $`\alpha `$-$`\alpha `$ spallation is similar to equation (98) (with $`Y_n`$ replaced by $`Y_\alpha `$). There is, however, one more condition. The shocked $`\alpha `$-particles are cooled by Coulomb interactions with $`e^{}`$ (or $`e^+`$) on a timescale $`\tau _{\mathrm{Coul}}`$, and an efficient spallation requires $`\tau _{\mathrm{life}}<\tau _{\mathrm{Coul}}`$ in addition to $`\tau _{\mathrm{life}}<(R/c\mathrm{\Gamma }_{\mathrm{max}})`$.

The $`e^\pm `$ can be considered as targets at rest for the hot ions because their radiative losses (synchrotron and/or inverse Compton) are rapid compared to the expansion rate. In the fluid frame, hot ions with mass $`m_i`$ and charge $`Ze`$ loose their random velocity $`\stackrel{~}{\beta }`$ on timescale (Ginzburg & Syrovatskii 1964)

$$\tau _{\mathrm{Coul}}=\stackrel{~}{\beta }\left(\frac{\mathrm{d}\stackrel{~}{\beta }}{\mathrm{d}\tau }\right)^1=\frac{2\stackrel{~}{\beta }^3m_i}{3Z^2\mathrm{ln}\mathrm{\Lambda }\sigma _\mathrm{T}m_ecn_e},$$
(99)
where $`n_e=n_{}+n_+`$ is the total density of $`e^{}`$ and $`e^+`$, and $`\mathrm{ln}\mathrm{\Lambda }20`$ is Coulomb logarithm. The internal shocks give $`\stackrel{~}{\beta }A`$ and, for $`\alpha `$-particles ($`m_i=4m_p`$, $`Z=2`$), we find $`\tau _{\mathrm{life}}=(n_\alpha \sigma _{\mathrm{sp}}c)^1<\tau _{\mathrm{Coul}}`$ if

$$\frac{n_e}{n_\alpha }<\frac{2A^3}{3\mathrm{ln}\mathrm{\Lambda }}\frac{m_p}{m_e}\frac{\sigma _{\mathrm{sp}}}{\sigma _\mathrm{T}}10A^3.$$
(100)
A crucial factor in this condition is the $`e^\pm `$ density. The postshock matter emits radiation, and $`e^\pm `$ pairs can be produced by $`\gamma `$-rays ($`\gamma +\gamma e^{}+e^+`$) that have energy $`h\nu >m_ec^2`$ in the fluid frame. Suppose a fraction $`\eta `$ of the internal energy density $`(A^2/2)n_bm_pc^2`$ is converted into radiation and a fraction $`f`$ of this radiation is above the threshold for pair production. A maximum pair density is evaluated assuming that all photons emitted above the threshold $`m_ec^2`$ are converted into pairs,

$$\frac{n_e}{n_b}f\eta \frac{A^2}{2}\frac{m_p}{m_e}.$$
(101)
From equations (100) and (101) we conclude that the Coulomb losses of $`\alpha `$-particles on $`e^\pm `$ can prevent the efficient $`\alpha `$-$`\alpha `$ spallation if $`f\eta >\kappa =3\times 10^3AX_\alpha `$. In this case, only a fraction $`\tau _{\mathrm{Coul}}/\tau _{\mathrm{life}}(\kappa /f\eta )`$ of the $`\alpha `$-particles are destroyed. The released neutrons can then continue the spallation process via $`n`$-$`\alpha `$ collisions (§ 4.2.1).

Equation (101) assumes an optical depth for $`\gamma `$-$`\gamma `$ interactions $`\tau _{\gamma \gamma }>1`$. We now check this assumption. The optical depth can be estimated as $`\tau _{\gamma \gamma }0.1(w_1/m_ec^2)\sigma _\mathrm{T}R/\mathrm{\Gamma }_{\mathrm{max}}`$ where $`w_1=f\eta (A^2/2)n_bm_pc^2`$. This yields

$$\tau _{\gamma \gamma }0.1f\eta \frac{A^2}{2}\frac{m_p}{m_e}n_b\sigma _\mathrm{T}\frac{R}{\mathrm{\Gamma }_{\mathrm{max}}}.$$
(102)
It is easy to see that $`\tau _{\gamma \gamma }>1`$ at all radii where efficient $`\alpha `$-$`\alpha `$ spallation can take place \[i.e. where $`\tau _{\mathrm{life}}<(R/c\mathrm{\Gamma }_{\mathrm{max}})`$\] if $`f\eta >5\times 10^4`$.

## 5 Conclusions

1. Reactions of $`e^\pm `$ capture on nucleons, $`e^{}+pn+\nu `$ and $`e^++np+\overline{\nu }`$, operate in GRB central engines and set an equilibrium proton fraction $`Y_e=n_p/(n_n+n_p)`$. If the engine is an accretion disk around a black hole of mass $`M`$, the equilibrium $`Y_e`$ is established at accretion rates $`\dot{M}>\dot{M}_{\mathrm{eq}}10^{31}(\alpha /0.1)^{9/5}(M/M_{})^{6/5}`$ g/s (eq. 48) where $`\alpha =0.010.1`$ is a viscosity parameter of the disk.

2. Of great importance for the explosion dynamics is whether $`Y_e<0.5`$, i.e., there is an excess of neutrons. A general analysis shows that a neutrino-transparent matter has equilibrium $`Y_e<0.5`$ if $`\mu >Q/2`$, and a similar condition for a $`\nu `$-opaque matter reads $`\mu >Q`$, where $`\mu `$ is the electron chemical potential in units of $`m_ec^2`$ and $`Q=(m_nm_p)/m_e=2.53`$. This condition is satisfied below a critical “neutronization” temperature $`T_n`$ (eqs.  and ). We find $`T<T_n`$ for plausible central engines of GRBs. In particular, accretion disks have a neutron excess at $`\dot{M}>\dot{M}_n10^{31}(M/M_{})^2(\alpha /0.1)`$ g/s (eq. 50).

3. Fireballs produced by neutron-rich engines should be also neutron rich. A major threat for neutrons in the escaping fireball is the neutrino flux from the central engine as absorption of a neutrino converts the neutron into a proton. This process is slower than the fireball expansion if the neutrino luminosity is below $`10^{53}`$ erg/s. A neutrino luminosity above $`10^{53}`$ erg/s would require a very powerful and $`\nu `$-opaque central engine, which is possible. In this case, however, absorption of $`\overline{\nu }`$ by the fireball protons takes place as well as absorption of $`\nu `$ by the neutrons. The balance between $`\nu `$ and $`\overline{\nu }`$ absorptions establishes a new equilibrium $`Y_e`$ in the fireball, which depends on the emitted spectra of $`\nu `$ and $`\overline{\nu }`$ and is likely below 0.5 (§ 2.3).

4. As the fireball expands and cools, the ejected free nucleons tend to recombine into $`\alpha `$-particles. This process competes, however, with rapid expansion and can freeze out. For this reason, nucleosynthesis is suppressed in fireballs with a high photon-to-baryon ratio $`\varphi =n_\gamma /n_b`$ (or, equivalently, high entropy per baryon $`s/k=3.6\varphi `$). We find that, in radial fireballs, more than half of nucleons can recombine only if $`\varphi <3\times 10^4(r_0/3\times 10^6)^{0.85}`$ where $`r_0`$ \[cm\] is the size of the central engine. In fireballs with parabolic collimation, the efficient recombination requires $`\varphi <6\times 10^5(r_0/3\times 10^6)^{0.9}`$. The typical GRB parameters $`\varphi 10^5`$ and $`r_03\times 10^6`$ cm are just marginal for nucleosynthesis.

5. Even in the case of efficient nucleon recombination, there are still leftover neutrons because of the neutron excess ($`Y_e<0.5`$). The minimum mass fraction of leftover neutrons is $`X_n=12Y_e`$.

6. The nucleosynthesis also produces deuterium, which is next abundant element after the free nucleons and $`\alpha `$-particles. Its typical mass fraction is a few per cent. The abundances of tritium, 3He, and all elements with mass number greater than 4 are negligible.

7. Synthesized $`\alpha `$-particles can be spalled later on, and then the populations of free neutrons and protons are increased. There are at least two possible mechanisms of spallation. (1) Energetic $`n`$-$`\alpha `$ collisions with a relative velocity $`\stackrel{~}{\beta }>\stackrel{~}{\beta }_{\mathrm{sp}}0.1`$ take place before the end of the fireball acceleration if $`\varphi >4\times 10^4(L_{\mathrm{th}}/L)^{3/4}`$. This mechanism works for fireballs with Lorentz factors $`\mathrm{\Gamma }_{\mathrm{max}}>160(kT_0/1\mathrm{MeV})`$ (eq. 92) where $`T_0`$ is an initial temperature (eq. 66). (2) Energetic $`n`$-$`\alpha `$ and $`\alpha `$-$`\alpha `$ collisions occur when the fireball is reheated by internal shocks. This mechanism can be efficient at modest $`\mathrm{\Gamma }_{\mathrm{max}}<300(L_\mathrm{\Omega }/10^{52})^{1/5}`$ (eq. 98) where $`L_\mathrm{\Omega }`$ \[erg/s\] is the fireball kinetic luminosity per unit solid angle.

The presence of a neutron component in GRB fireballs has quite spectacular implications. Beside making the fireball an interesting source of multi-GeV neutrinos (Derishev et al. 1999, Bahcall & Mészáros 2000, Mészáros & Rees 2000) the neutrons survive and play a dramatic role for the explosion development at large radii $`R10^{16}10^{17}`$ cm. When the fireball begins to decelerate as a result of the interaction with an external medium, neutrons continue to coast with a high Lorentz factor $`\mathrm{\Gamma }_n`$ and form a leading front. One can easily show that this front has kinetic energy much larger (by a factor of $`X_n\mathrm{\Gamma }_n\mathrm{\Gamma }_{\mathrm{max}}`$) than the rest-mass of the ambient medium. It leaves behind a relativistic trail loaded with the products of the neutron decay until the neutron front decays completely, which happens at $`R10^{17}`$ cm (about $`\mathrm{log}X_n\mathrm{\Gamma }_n\mathrm{\Gamma }_{\mathrm{max}}10`$ times the mean decay radius). An external shock wave — the customary source of GRB afterglows — has to form in the neutron trail rather than in a normal ambient medium. The mechanism of the fireball deceleration in this situation is elaborated in Paper 2.

The neutron features of GRB explosions appear inevitably in the standard fireball scenario. They would be absent, however, if the GRBs are produced by magnetized winds with extremely low baryon loading, where Poynting flux carries much more energy than matter (e.g. Usov 1994, Lyutikov & Blandford 2002). A detection or non-detection of neutron effects will constrain the level of baryon loading in GRBs.

At the final stages of the preparation of this manuscript, a paper by Pruet, Woosley, & Hoffman (2002) appeared. They evaluate $`Y_e`$ for numerical accretion models of Popham et al. (1999). The results are consistent with our analysis in § 2 (note a slightly different definition of the viscosity parameter $`\alpha _{\mathrm{our}}=[3/2]\alpha _{\mathrm{Popham}}`$).

I thank Raymond Sawyer for discussions of the neutronization process and Matthias Liebendoerfer for an introduction to the Lattimer-Swesty Equation of State. I acknowledge the hospitality of the ITP at Santa Barbara, where part of this work was done. This research was supported by NSERC, and in part by the NSF grant PHY99-07949 and RFBR grant 00-02-16135.


text:
# Structure of Multi-meron Knot Action

## 1 Introduction

Faddeev’s conjecture  that the energy of the ground state and properties of low-energy excitations in the non-Abelian field theory are determined by topological invariants of knots and links has been recently developed  using the $`𝐧`$-field model. One of these topological parameters is the Hopf invariant, which determines the knotting degree of filamental manifolds, where the unit vector $`𝐧`$ is defined.

The study of the behavior of a vortex filament tangle is an active area of research in the field theory and attracts attention due to several reasons . The topological order associated with linking at short distances exists against the background of disorder due to free motion of the tangle filaments. Therefore, such systems of entangled filaments contain data on their behavior in the ultraviolet, as well as infrared limit, while systems of point particles do not contain. In the present paper we discuss some properties of field configurations within the $`SU(2)`$ Yang-Mills theory and within the $`CP^1`$ GL model from this point of view.

## 2 The $`SU(2)`$ Yang-Mills theory

Let us consider the Langrangian density,

$$L=\frac{1}{4g^2}Tr\widehat{F}_{\mu \nu }^2,$$
(1)
$$\widehat{F}_{\mu \nu }=_\mu \widehat{A}_\nu _\nu \widehat{A}_\mu +[\widehat{A}_\mu ,\widehat{A}_\nu ],\widehat{A}_\mu =A_\mu ^a\frac{\sigma ^a}{2i}.$$
(2)
We will use special anzatz  for the potential $`A_\mu ^a`$:

$$A_\mu ^a=\epsilon ^{abc}_\mu n^bn^c(K+1).$$
(3)
Here $`g`$ is the bare coupling constant, $`K(𝐫,t)`$ is the uknown scalar function, $`𝐧`$ is the unit vector field, and $`\sigma ^a`$ denotes Pauli matrices. Using Eq.(3), one can get

$`L={\displaystyle \frac{1}{2g^2}}\left\{\delta _{\mu \nu }(_\lambda 𝐧)^2(_\mu 𝐧)(_\nu 𝐧)\right\}^\mu K^\nu K+`$
$`{\displaystyle \frac{1}{4g^2}}\left\{𝐧[_\mu 𝐧,_\nu 𝐧]\right\}^2\left(K^21\right)^2.`$
The Lagrangian density in the complete parametrization of the potential $`A_\mu ^a`$ was computed in Refs. It is seen from Eq.(4) that multipliers in curly brackets of the kinetic and potential parts of Eq.(4) play the role of the coupling constants for the $`K`$-field and *vice versa*. This takes place when $`𝐧𝐫/r`$ and the foregoing soft version is applicable. Let us fix the $`𝐧`$-field dynamics by the condition $`𝐧=𝐫/r`$. Then we will get a significantly simplified problem. Indeed, in this case $`A_k^a=\epsilon ^{akc}n^c(K+1)/r`$ and

$$\left\{\delta _{\mu \nu }(_\lambda 𝐧)^2(_\mu 𝐧)(_\nu 𝐧)\right\}=\frac{\delta _{\mu \nu }}{r^2},(𝐧[_\mu 𝐧,_\nu 𝐧])^2=\frac{1}{r^4}.$$

The Lagrangian density and the equation of motion for the scalar field $`K(r,t)`$ takes the form

$$L=\frac{1}{2g^2r^2}(_\mu K)^2\frac{1}{4g^2r^4}\left(K^21\right)^2,$$
(5)
$$r^2\mathrm{}K=K(K^21),\mathrm{}=\frac{^2}{t^2}\frac{^2}{r^2}.$$
(6)
The simplest solution  is the meron configuration $`K=t/\sqrt{t^2r^2}.`$ Multi-meron configurations of the $`K`$-field were studied in detail in Ref. Common property of the solutions of Eq.(6) is existence of singularities of the field $`K`$ on closed surfaces located at a finite distance from the point $`𝐫=0`$. Their existence makes it impossible to consider the infrared limit in the foregoing case, when degrees of freedom of the $`𝐧`$-field are frozen.

A fixed scalar degree of freedom, which is connected with the $`K`$\- field can be obtained by integrating over this variable in the Feynman integral. This problem was investigated in Refs. As a result, the main contribution to the model in Ref. is characterized by the Lagrangian density

$$L=c_1(_\mu 𝐧)^2+c_2\left(𝐧[_\mu 𝐧,_\nu 𝐧]\right)^2,$$
(7)
where $`c_1`$ and $`c_2`$ are effective coupling constants. The first term in Eq.(8) describes the infrared limit of the $`𝐧`$-field dynamics, while the second determines the behavior at a short distance, where topological effects of links are important.

Now let us proceed to the general case. This can be done by considering the self-dual equations using the anzatz in Eq.(3). Note that in this case of low symmetry  these equations differ from those obtained for spherically-symmetrical  or axially-symmetrical  cases. Thus, we have

$$[_\tau 𝐧,_k𝐧]=\frac{1}{2}\epsilon _{kms}[_m𝐧,_s𝐧].$$
(8)
$$_\tau K_k𝐧_kK_\tau 𝐧=\frac{1}{2}\epsilon _{kms}\left\{_mK_s𝐧_sK_m𝐧\right\},$$
(9)
where $`\tau =it`$ is the Euclidean time. Since the homotopic group $`\pi _4(SU(2))=_2`$ is non-trivial , we hope that there exists at least one non-trivial class of configurations of fields $`K`$ and $`𝐧`$.

## 3 The $`CP^1`$ Ginsburg-Landau model

We will use the GL functional,

$`F={\displaystyle }d^3x\left[{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{1}{2m}}\right|(\mathrm{}_k+i{\displaystyle \frac{2e}{c}}A_k)\mathrm{\Psi }_\alpha |^2+`$
$`{\displaystyle \underset{\alpha }{}}(b_\alpha |\mathrm{\Psi }_\alpha |^2+{\displaystyle \frac{c_\alpha }{2}}|\mathrm{\Psi }_\alpha |^4)+{\displaystyle \frac{𝐁^\mathrm{𝟐}}{8\pi }}],`$
with a two-component order parameter $`\mathrm{\Psi }_\alpha =\sqrt{2m}\rho \chi _\alpha ,\chi _\alpha =|\chi _\alpha |e^{i\phi _\alpha },`$ which satisfies the $`CP^1`$ condition, $`|\chi _1|^2+|\chi _2|^2=1`$. In the (3+0)-dimensional case, Eq.(10) has the meaning of free energy and $`\mathrm{\Psi }_\alpha `$ is a two-component order parameter, which is used either in the context of two-gap supercondutivity  or in the non-Abelian field theory as a Higgs doublet . In (2+1)D systems, Eq.(10) has the meaning of action, and the order parameter $`\mathrm{\Psi }_\alpha `$ gives a two-dimensional non-Abelian representation of the braid group, which is a specific group of symmetry for planar systems at permutations of particles inside them.

As was shown in Ref., there exists exact mapping of model (10) into the following version of the $`𝐧`$-field model:

$`F=`$
$`{\displaystyle d^3x\left[\frac{1}{4}\rho ^2\left(_k𝐧\right)^2+\left(_k\rho \right)^2+\frac{1}{16}\rho ^2𝐜^2+\left(F_{ik}H_{ik}\right)^2+V(\rho ,n_3)\right]}.`$
The equation (11) was obtained with the use of gauge invariant order parameter fields of the unit vector $`n^a=\overline{\chi }\sigma ^a\chi `$, where $`\overline{\chi }=(\chi _1^{},\chi _2^{})`$ and the velocity $`𝐜=𝐉/\rho ^2`$. The total current, $`𝐉=2\rho ^2(𝐣4𝐀)`$, has paramagnetic $`(𝐣=i[\chi _1\chi _1^{}c.c.+(12)])`$ and diamagnetic $`(4𝐀)`$ parts. Besides, we used in Eq.(11) the following notations: $`F_{ik}=_ic_k_kc_i`$, and $`H_{ik}=𝐧[_i𝐧\times _k𝐧]:=_ia_k_ka_i`$ and the dimensionless units: $`L=(\xi _1+\xi _2)/2`$ as the unit of length, $`\xi _\alpha =\mathrm{}/\sqrt{2mb_\alpha }`$ as the coherence length; $`\mathrm{}/L`$ as the unit of c; $`c^2/(512\pi e^2L^2)`$ as the unit of $`\rho ^2`$ and, finally, $`\gamma /L`$ with $`\gamma =\left(c\mathrm{}/e\right)^2/512\pi `$ as the unit of energy.

Let us now enumerate some non-trivial ($`𝐧const`$) situations that follow from Eq.(11) as a result of competition of the order parameters $`\rho `$, $`𝐧`$, and $`𝐜`$: (i) $`𝐜=0`$, $`\rho =const`$; (ii) $`𝐜=0`$, $`\rho const`$; (iii) $`𝐜0`$, $`\rho =const`$; (iv) $`𝐜0`$, $`\rho const`$. In the case when $`𝐜0`$, $`\rho const`$, and $`𝐧=const`$ we have the problem of a vortex state structure in the classical GL model.

In the limit (i) we get Eq.(7) . A numerical study of knotted configurations of the $`𝐧`$-field in this model was performed in Refs. . The lower energy bound in this case,

$$F32\pi ^2|Q|^{3/4},$$
(12)
is determined  by the Hopf invariant,

$$Q=\frac{1}{16\pi ^2}d^3x\epsilon _{ikl}a_i_ka_l.$$
(13)
If the compactification $`^3S^3`$ is used and $`𝐧S^2`$, then the integer $`Q\pi _3(S^2)=`$ shows the degree of linking or knotting of filamental manifolds, where the vector field $`𝐧(x,y,z)`$ is defined. In particular, for two linked rings (Hopf linking) $`Q=1`$, for the trefoil knot $`Q=6`$, etc. It is important to note that $`\pi _3(CP^M)=0`$ at $`M>1`$ and $`\pi _3(CP^1)=\pi _3(S^2)=`$ . In the latter case, the order parameter is two-component , and linked or knotted soliton configurations are labeled by the Hopf invariant in Eq.(13).

The characteristic size $`R_Q`$ of a knot can be found using Eq.(12) and the minimum value of the free energy $`F_{min}=2\sqrt{c_1c_2}`$ for the radius $`R_{min}=\sqrt{c_2/c_1}`$. Therefore,

$$R_Q=\frac{c_2}{16\pi ^2\left|Q\right|^{3/4}}.$$
(14)
We see that the characteristic size of such a knot is determined by a combination of dynamical ($`c_2`$) and topological ($`Q`$) features of the system. In case (ii), we have the problem of a soft version of the $`𝐧`$-field model . Case (iv) is most general.

We will focus now on limit (iii). Let us assume that $`\rho `$ can be found from the minimum value of the potential $`V(\rho )`$ and the velocity $`𝐜`$ is not equal to zero. Eq.(11) in this case takes the following form:

$`F=F_n+F_cF_{int}=`$
$`{\displaystyle d^3x\left[\left(\left(_k𝐧\right)^2+H_{ik}^2\right)+\left(\frac{1}{4}𝐜^2+F_{ik}^2\right)2F_{ik}H_{ik}\right]}.`$
It is seen from Eq.(15) that the superconducting state with $`𝐜0`$ has the energy, which is less than the minimum in the case (i) due to renormalization of the coefficient $`(=1)`$ in the second term of the functional $`F_n`$. In order to find the lower free energy bound in the superconducting state with $`𝐜0`$, we will use the following inequality :

$$F_n^{5/6}F_c^{1/2}(32\pi ^2)^{4/3}|L|,$$
(16)
where

$$L=\frac{1}{16\pi ^2}d^3x\epsilon _{ikl}c_i_ka_l$$
(17)
is the degree of mutual linking  of the velocity $`𝐜`$ lines and of the magnetic field $`𝐇=[\times 𝐚]`$ lines. It is also an integral of motion.

Using the Schwartz-Cauchy-Bunyakovsky inequality, one will find  that $`FF_{min}=(F_n^{1/2}F_c^{1/2})^2.`$ Using Eq.(16) and the condition $`Q0`$ the last inequality can be rewritten in the equivalent form

$$F32\pi ^2|Q|^{3/4}(1|L|/|Q|)^2.$$
(18)
The minimum value of free energy in Eqs.(18) corresponds to the field configurations with $`L=Q0`$ or $`𝐚=𝐜`$. They satisfy the self-dual condition $`F_n=F_c`$.

It is clear from Eq.(18) that for all numbers $`L<Q`$ the energy of the ground state is less than that described in the model , for which the inequality (12) is valid. The origin of the energy decrease can be easily understood. Even under the conditions of existence of the paramagnetic part $`𝐣`$ of the current $`𝐉`$, the diamagnetic interaction in the superconducting state cancels the own current energy and a part of the energy related to the $`𝐧`$-field dynamics for all state classes with $`L<Q`$. It is also important that the current (total momentum of a superconducting pair) is non-zero in the superconducting state. In this regard, the considered inhomogeneous state with the current is similar to the state  investigated in Refs. .

## 4 Discussion

Up to now the vector $`𝐀`$ has characterized the internal charge $`U(1)`$ gauge symmetry. If we apply the external electromagnetic field, then the vector potential $`𝐀`$ will be equal to the sum of internal and external gauge potentials. As a result, the velocity $`𝐜`$ decreases due to the diamagnetism of the superconducting state. This leads to suppression of the superconducting gap. The answer to the question on existence of full or partial Meisner screening in these states depends on the result of the competition between contributions from neutral $`𝐣`$ and charged $`4𝐀`$ parts to the total current $`𝐉`$.

In the $`(3+0)D`$ case of free energy (10), Hopf invariant (13) is analogous to the Chern-Simons action $`(k/4\pi )𝑑td^2x\epsilon _{\mu \nu \lambda }a_\mu _\nu a_\lambda `$, which determines strong correlations of $`(2+1)D`$ states  for $`k2`$. In planar systems, this coefficient has the sense of braiding degree of excitation world lines. In particular, for the semion $`k=2`$. Keeping in mind the relation of spatial dimensionality of the systems in their quantum and statistical descriptions, we note that the $`(2+1)D`$ dynamical case $`k=2`$ of open world line ends of excitations is equivalent to the compact $`(3+0)D`$ statistical example of the Hopf linking $`Q=1`$.

In conclusion, we have considered the structure of the multi-meron knot action and found the self-dual equation. Using the topological invariants of knots for the analysis of coherent phase states, we have also determined the conditions of applying the theory significantly based on the dynamics and topology of the $`𝐧`$-field model.

## Acknowledgements

We are very grateful to A. G. Abanov, S. A. Brazovskii, L. D. Faddeev, A. I. Larkin, Y.-S. Wu and G. E. Volovik for useful discussions. One of the authors (AP) thanks Prof. M. Boiti for kind hospitality. This work was partially supported by the RFBR under Grant # 01-02-17225.