text:
# 1 Introduction

## 1 Introduction

Spin dependent deep-inelastic lepton-hadron scattering has received much interest during the past years, see Refs. for reviews. Presently, next-to-leading order analyses of the experimental data for the structure function $`g_1`$ allow for the extraction of polarized parton distributions. However, these are considerably less constrained than their unpolarized counterparts. This is particularly true for the gluon spin distribution $`\mathrm{\Delta }\varphi _{g/P}`$, which is currently determined from scaling violations. Numerous polarized hard scattering processes are sensitive to the gluon spin distribution, and the first moment of $`\mathrm{\Delta }\varphi _{g/P}`$ contributes to the sum rule for the nucleon’s spin content .

Deep-inelastic leptoproduction of heavy quarks proceeds predominantly through the boson-gluon fusion process, and as such yields important information about the gluon distribution inside the proton. Polarized leptoproduction of charm quarks is therefore promising for an experimental determination of $`\mathrm{\Delta }\varphi _{g/P}`$. Presently, this measurement can be realized by the HERMES experiment . The upcoming COMPASS experiment  is dedicated to measure $`\mathrm{\Delta }\varphi _{g/P}`$ directly through the boson-gluon fusion process and subsequent polarized open charm production. COMPASS will probe both the region of photoproduction and deep-inelastic scattering and extend the kinematical range accessible to the HERMES experiment.

In general, fixed target experiments, such as HERMES or COMPASS with centre-of-mass energies of $`\sqrt{S}=7.4\mathrm{GeV}`$ and $`14.1\mathrm{GeV}`$ respectively, operate rather close to the threshold for charm pair-production. In this region of phase space, the perturbatively calculable hard scattering cross section has potentially large higher order double-logarithmic threshold corrections. Thus, a resummation of these Sudakov logarithms to all orders in perturbation theory is required in order to regain control over the perturbative expansion.

The purpose of this letter is to carry out this Sudakov resummation, which is done in section 2. We work to next-to-leading logarithmic (NLL) accuracy (see Ref. for a recent review), and our final result provides an approximate sum of the complete perturbative expansion if, at each order, the Sudakov corrections dominate. Previously, studies of threshold resummation for polarized scattering processes have been performed for the polarized Drell-Yan cross section in Ref.. By re-expanding our resummed cross section in $`\alpha _s`$ to next-to-leading (NLO) and next-to-next-to-leading order (NNLO), we provide in section 3 estimates of the size of the exact higher order corrections111For the exact NLO corrections to polarized heavy quark photoproduction cf. Refs.., the calculation of which has not yet been completed . Our estimates are valid in the region dominated by threshold logarithms and we study their impact on the charm structure function $`g_1^\mathrm{c}`$. In the kinematical region relevant to HERMES and COMPASS, we investigate the sensitivity to the gluon spin distribution $`\mathrm{\Delta }\varphi _{g/P}`$ and the dependence on the factorization scale. Section 4 contains our conclusions.

## 2 Resummed differential cross section

We begin with the definition of the exact single-particle inclusive (1PI) kinematics. We study the reaction of lepton ($`\mathrm{l}`$) – proton ($`\mathrm{P}`$) scattering,

$$\mathrm{l}(l)+\mathrm{P}(p)\mathrm{l}(lq)+\mathrm{Q}(p_1)+X[\overline{\mathrm{Q}}](p_2^{}),$$
(1)
where $`\mathrm{Q}`$ is a heavy quark and $`X[\overline{\mathrm{Q}}]`$ denotes any allowed hadronic final state containing at least the heavy antiquark. Neglecting electroweak radiative corrections222 QED radiative corrections have been studied in the leading-logarithmic approximation in Ref.., for longitudinally polarized protons, the difference of cross sections for reaction (1) may after summation over $`X`$ be written as

$`{\displaystyle \frac{d^5\sigma ^{}}{dxdyd\varphi dT_1dU_1}}{\displaystyle \frac{d^5\sigma ^{}}{dxdyd\varphi dT_1dU_1}}=`$
$`{\displaystyle \frac{4\alpha ^2}{Q^2}}\left[\left(2y{\displaystyle \frac{2M^2x^2y^2}{Q^2}}\right){\displaystyle \frac{d^2g_1(x,Q^2,m^2,T_1,U_1)}{dT_1dU_1}}{\displaystyle \frac{4M^2x^2y^2}{Q^2}}{\displaystyle \frac{d^2g_2(x,Q^2,m^2,T_1,U_1)}{dT_1dU_1}}\right],`$
where $`\alpha `$ and $`M`$ are the fine structure constant and the proton mass, and $`\varphi `$ denotes the angle between the proton spin and the momentum of the scattered lepton, $`(lq)`$ in Eq.(1). The lower arrow on $`d^5\sigma `$ indicates the polarization of the incoming lepton in the direction of its momentum $`l`$, whereas the upper arrow on $`d^5\sigma `$ represents the polarization of the proton, which is parallel or anti-parallel to the polarization of the incoming lepton.

The functions $`d^2g_k/dT_1dU_1,k=1,2`$ are the double-differential deep-inelastic polarized heavy quark structure functions. The kinematic variables $`Q^2,x,y`$ are defined by

$$Q^2=q^2>0,x=\frac{Q^2}{2pq},y=\frac{pq}{pl}.$$
(3)
Additionally, we also have the overall invariants

$$S=(p+q)^2S^{}Q^2,T_1=(pp_1)^2m^2,U_1=(qp_1)^2m^2,S_4=S+T_1+U_1+Q^2.$$
(4)
Only the structure function $`g_1`$ is relevant asymptotically for $`M/Q1`$. It contains leading twist two operators and has a parton model interpretation. We therefore restrict ourselves to $`g_1`$ only. Furthermore, we neglect the contributions from initial state quarks and antiquarks. This is justified as they vanish at leading order. We consider only the gluon-initiated partonic subprocess

$`\gamma ^{}(q)+g(k)`$ $``$ $`\mathrm{Q}(p_1)+X^{}[\overline{\mathrm{Q}}](p_2^{}),`$ (5)
where $`k=zp`$. The partonic invariants are

$$s=(k+q)^2s^{}+q^2,t_1=(kp_1)^2m^2,u_1=(qp_1)^2m^2,s_4=(p_2^{})^2m^2.$$
(6)
The invariant $`s_4`$ parametrizes the inelasticity of the partonic reaction (5) in 1PI kinematics. This distance from threshold is conveniently measured in terms of a dimensionless weight function $`w_S`$ ,

$$w_S=\frac{s_4}{m^2}\frac{2\overline{p}_2k_S}{m^2}\frac{2\zeta k_S}{m}.$$
(7)
The 1PI kinematics vector $`\zeta `$ is defined by splitting the recoil momentum $`p_2^{}=\overline{p}_2+k_S`$ into the momentum $`\overline{p}_2`$ at threshold, and the momentum $`k_S`$ of any additional (soft) radiation above threshold.

In double-differential form, $`d^2g_1/dT_1dU_1`$ factorizes as

$`S^2{\displaystyle \frac{d^2g_1(x,Q^2,m^2,T_1,U_1)}{dT_1dU_1}}`$ $`=`$ $`{\displaystyle \underset{z^{}}{\overset{1}{}}}{\displaystyle \frac{dz}{z}}\mathrm{\Delta }\varphi _{g/P}(z,\mu ^2)\mathrm{\Delta }\omega _g({\displaystyle \frac{x}{z}},s_4,t_1,u_1,Q^2,m^2,\mu ^2,\alpha _s(\mu )).`$ (8)
The polarized gluon distribution in the proton is denoted by $`\mathrm{\Delta }\varphi _{g/P}`$ and $`z`$ is the momentum fraction with $`z^{}=U_1/(S^{}+T_1)`$. The dimensionless function $`\mathrm{\Delta }\omega _g`$ describes the hard scattering process and depends on the partonic invariants $`t_1,u_1`$ and $`s_4`$, defined in Eq.(6). The factorization scale is denoted by $`\mu `$ and throughout this paper taken equal to the renormalization scale.

Sufficiently close to threshold, $`g_1`$ in Eq.(8) is dominated by those higher order contributions that contain plus-distributions of the type

$$\left[\frac{\mathrm{ln}^l(s_4/m^2)}{s_4}\right]_+=\underset{\mathrm{\Delta }0}{lim}\left\{\frac{\mathrm{ln}^l(s_4/m^2)}{s_4}\theta (s_4\mathrm{\Delta })+\frac{1}{l+1}\mathrm{ln}^{l+1}\left(\frac{\mathrm{\Delta }}{m^2}\right)\delta (s_4)\right\}.$$
(9)
They result from imperfect cancellations between soft and virtual contributions to the cross section. At order $`𝒪(\alpha _s^{i+2}),i=0,1,\mathrm{}`$ we refer to the corrections with $`l=2i+1`$ as leading logarithmic (LL), to the ones with $`l=2i`$ as next-to-leading logarithmic (NLL), etc.

The Sudakov resummation of these singular functions in Eq.(9) uses the methods and results of Refs.. It rests upon the factorization of $`\mathrm{\Delta }\omega _g`$ into separate functions $`\mathrm{\Delta }\psi _{g/g}`$, $`S`$, and $`H_g`$ for the jetlike-, soft, and off-shell quanta. This factorization, valid in the threshold region of phase space, implies a decomposition  of the total weight function $`w`$,

$`w=w_\psi +w_S=w_1\left({\displaystyle \frac{u_1}{m^2}}\right)+w_S,`$ (10)
which expresses the overall weight in terms of the individual weights $`w_\psi `$ and $`w_S`$. Each of the functions $`\mathrm{\Delta }\psi _{g/g}(w_1)`$, $`S(w_S)`$, and $`H_g`$, which are conveniently computed in $`\zeta A=0`$ gauge, organizes contributions corresponding to a particular set of quanta and thereby only depends on its own individual weight function.

It is convenient to consider these functions in moment space, defined by the Laplace transform with respect to the overall weight $`w`$,

$$\stackrel{~}{f}(N)=\underset{0}{\overset{\mathrm{}}{}}𝑑we^{Nw}f(w).$$
(11)
Under the assumption of factorization the partonic cross section $`\mathrm{\Delta }\omega _g`$ in moment space can then be written in a factorized form, up to corrections of order $`1/N`$, as

$`\mathrm{\Delta }\stackrel{~}{\omega }_g(N,t_1,u_1,Q^2,m^2,\mu ^2,\alpha _s(\mu ))=`$
$`H_g(\zeta ,Q^2,m^2,\alpha _s(\mu ))\left[{\displaystyle \frac{\mathrm{\Delta }\stackrel{~}{\psi }_{g/g}(N_u,k\zeta /\mu )}{\mathrm{\Delta }\stackrel{~}{\varphi }_{g/g}(N_u,\mu )}}\right]\stackrel{~}{S}({\displaystyle \frac{m}{N\mu }},\zeta ).`$
where $`N_u=N(u_1/m^2)`$. The $`N`$-dependence in each of the functions of Eq.(2) exponentiates.

The jet and soft functions, $`\mathrm{\Delta }\psi _{g/g}`$ and $`S`$ in Eq. (2), can each be represented as operator matrix elements . Dependence on the spin degrees of freedom can be kept explicit in these elements, such that the methods of Refs. for the resummation of the singular functions in Eq.(9) via appropriate evolution equations can be generalized in a straightforward way to account for initial state polarizations. The jet-function $`\mathrm{\Delta }\stackrel{~}{\psi }_{g/g}`$ obeys two evolution equations, namely the renormalization group equation and an equation describing the energy dependence via gauge-dependence . Solving both equations sums Sudakov double logarithms. Mass factorization introduces the density $`\mathrm{\Delta }\stackrel{~}{\varphi }_{g/g}`$ as a counterterm and requires the choice of a scheme. In the $`\overline{\mathrm{MS}}`$-scheme $`\mathrm{\Delta }\stackrel{~}{\varphi }_{g/g}`$ has no double logarithms. The soft function $`S`$ obeys a renormalization group equation, which resums all NLL logarithms originating from soft gluons not already accounted for in $`\mathrm{\Delta }\stackrel{~}{\psi }_{g/g}`$. In general, the soft and the hard function $`S`$ and $`H_g`$ depend also on the colour structure of the underlying scattering reaction, but for deep-inelastic production of heavy quarks this dependence is trivial .

The final result for the hard scattering function $`\mathrm{\Delta }\stackrel{~}{\omega }_g`$ in the $`\overline{\mathrm{MS}}`$-scheme in moment space resums all $`\mathrm{ln}N`$ in single-particle inclusive kinematics. We obtain

$`\mathrm{\Delta }\stackrel{~}{\omega }_g^{\mathrm{res}}(N,t_1,u_1,Q^2,m^2,\mu ^2,\alpha _s(\mu ))`$ $`=`$
$`H_g(\zeta ,Q^2,m^2,\alpha _s(m))\stackrel{~}{S}(1,\zeta )\mathrm{exp}\left\{2{\displaystyle \underset{\mu }{\overset{m}{}}}{\displaystyle \frac{d\mu ^{}}{\mu ^{}}}\gamma _g\left(\alpha _s(\mu ^{})\right)\right\}`$
$`\times \mathrm{exp}\left\{{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{dw}{w}}(1\mathrm{e}^{N_uw})\left[{\displaystyle \underset{w^2}{\overset{1}{}}}{\displaystyle \frac{d\lambda }{\lambda }}A^g(\alpha _s(\sqrt{\lambda }2k\zeta ))+{\displaystyle \frac{1}{2}}\nu _g(\alpha _s(w2k\zeta ))\right]\right\}`$
$`\times \mathrm{exp}\left\{{\displaystyle \underset{m}{\overset{m/N}{}}}{\displaystyle \frac{d\mu ^{}}{\mu ^{}}}2\mathrm{Re}\mathrm{\Gamma }_S\left(\alpha _s(\mu ^{})\right)\right\}\mathrm{exp}\left\{2{\displaystyle \underset{\mu }{\overset{2k\zeta }{}}}{\displaystyle \frac{d\mu ^{}}{\mu ^{}}}\left(\gamma _g\left(\alpha _s(\mu ^{})\right)\gamma _{g/g}(N_u,\alpha _s(\mu ^{}))\right)\right\},`$
where $`N_u`$ is defined below Eq.(2). To NLL accuracy, as defined below Eq.(9), the product $`H_gS`$ on the second line of Eq.(2) is determined from matching to the Born result at the scale $`\mu =m`$. To this accuracy the product $`H_gS`$ is also insensitive to the choice of treatment of $`\gamma _5`$.

The second exponent in Eq.(2) gives the leading $`N`$-dependence of the ratio $`\mathrm{\Delta }\stackrel{~}{\psi }_{g/g}/\mathrm{\Delta }\stackrel{~}{\varphi }_{g/g}`$ with

$`A_g(\alpha _s)=C_A{\displaystyle \frac{\alpha _s}{\pi }}+{\displaystyle \frac{1}{2}}C_AK\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2+\mathrm{},\nu ^g(\alpha _s)=\mathrm{\hspace{0.17em}2}C_A{\displaystyle \frac{\alpha _s}{\pi }}+\mathrm{},`$ (14)
where $`K=C_A(67/18\pi ^2/6)5/9n_f`$ can be inferred from Refs. using Ref. (cf. also Ref.). The scale evolution of the ratio $`\mathrm{\Delta }\stackrel{~}{\psi }_{g/g}/\mathrm{\Delta }\stackrel{~}{\varphi }_{g/g}`$ is governed by

$`\gamma _g(\alpha _s)`$ $`=`$ $`b_2{\displaystyle \frac{\alpha _s}{\pi }}+\mathrm{},`$ (15)
$`\gamma _{g/g}(N,\alpha _s)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{\pi }}\left(C_A\mathrm{ln}Nb_2\right)\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2\left({\displaystyle \frac{1}{2}}C_AK\mathrm{ln}N\right)+\mathrm{},`$ (16)
with $`\gamma _{g/g}`$ calculated in Refs., and the soft anomalous dimension to order $`\alpha _s`$ is

$`\mathrm{\Gamma }_S(\alpha _s)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{\pi }}\{({\displaystyle \frac{C_A}{2}}C_F)(L_\beta +1){\displaystyle \frac{C_A}{2}}(\mathrm{ln}\left({\displaystyle \frac{4(k\zeta )^2}{m^2}}\right)+\mathrm{ln}{\displaystyle \frac{m^4}{t_1u_1}})\}+\mathrm{},`$ (17)
with $`\beta =\sqrt{14m^2/s}`$ and $`L_\beta =(12m^2/s)\{\mathrm{ln}(1\beta )/(1+\beta )+\mathrm{i}\pi \}/\beta `$. We note that the various anomalous dimensions in Eqs.(14)–(17) are identical to NLL accuracy to the corresponding quantities of soft gluon resummation for unpolarized scattering .

Eq.(2) represents the central result of this paper. It provides the sum of Sudakov logarithms due to soft gluon emission to all orders in the perturbative expansion and accurate to the next-to-leading logarithm.

## 3 Finite order results

We expand our resummed result in $`\alpha _s`$ up to second order so as to provide NLO and NNLO approximations to NLL accurracy for the partonic single-particle inclusive double-differential cross section difference $`\mathrm{\Delta }\sigma _{\gamma ^{}g}`$.

For the process $`\gamma ^{}+g\mathrm{Q}+\overline{\mathrm{Q}}`$, the single-particle inclusive double-differential cross section difference is related to the hard function $`\mathrm{\Delta }\omega _g`$ of Eq.(8). To NLL accuracy it can be written in a factorized form as

$`s^{\mathrm{\hspace{0.17em}2}}{\displaystyle \frac{d^2\mathrm{\Delta }\sigma _{\gamma ^{}g}(s^{},t_1,u_1)}{dt_1du_1}}=\mathrm{\hspace{0.17em}8}\pi ^2\alpha {\displaystyle \frac{x}{Q^2}}\mathrm{\Delta }\omega _g(x,s_4,t_1,u_1,Q^2,m^2,\mu ^2,\alpha _s(\mu ))`$
$`\mathrm{\Delta }B_{\gamma ^{}g}^{\mathrm{Born}}(s^{},t_1,u_1)\left\{\delta (s^{}+t_1+u_1)+{\displaystyle \frac{\alpha _s(\mu )}{\pi }}K^{(1)}(s^{},t_1,u_1)+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)^2K^{(2)}(s^{},t_1,u_1)\right\},`$
with the Born level hard part $`\mathrm{\Delta }B_{\gamma ^{}g}^{\mathrm{Born}}`$ given by

$`\mathrm{\Delta }B_{\gamma ^{}g}^{\mathrm{Born}}(s^{},t_1,u_1)=`$
$`\alpha \alpha _se_\mathrm{q}^2\pi \left[{\displaystyle \frac{t_1}{u_1}}{\displaystyle \frac{u_1}{t_1}}2m^2\left({\displaystyle \frac{1}{t_1}}+{\displaystyle \frac{1}{u_1}}+{\displaystyle \frac{t_1}{u_1^2}}+{\displaystyle \frac{u_1}{t_1^2}}\right)+2q^2\left({\displaystyle \frac{1}{t_1}}+{\displaystyle \frac{1}{u_1}}+{\displaystyle \frac{2}{s^{}}}\right)\right].`$
The NLO and NNLO soft gluons corrections to NLL accuracy are

$`K^{(1)}(s^{},t_1,u_1)=`$
$`2C_A\left[{\displaystyle \frac{\mathrm{ln}(s_4/m^2)}{s_4}}\right]_++\left[{\displaystyle \frac{1}{s_4}}\right]_+\left\{C_A\left(\mathrm{ln}\left({\displaystyle \frac{t_1}{u_1}}\right)+\mathrm{Re}L_\beta \mathrm{ln}\left({\displaystyle \frac{\mu ^2}{m^2}}\right)\right)2C_F\left(\mathrm{Re}L_\beta +1\right)\right\}`$
$`+\delta (s_4)C_A\mathrm{ln}\left({\displaystyle \frac{u_1}{m^2}}\right)\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{m^2}}\right),`$
$`K^{(2)}(s^{},t_1,u_1)=`$
$`2C_A^2\left[{\displaystyle \frac{\mathrm{ln}^3(s_4/m^2)}{s_4}}\right]_++\left[{\displaystyle \frac{\mathrm{ln}^2(s_4/m^2)}{s_4}}\right]_+\{3C_A^2(\mathrm{ln}\left({\displaystyle \frac{t_1}{u_1}}\right)+\mathrm{Re}L_\beta \mathrm{ln}\left({\displaystyle \frac{\mu ^2}{m^2}}\right))`$
$`2C_A(b_2+3C_F(\mathrm{Re}L_\beta +1))\}+[{\displaystyle \frac{\mathrm{ln}(s_4/m^2)}{s_4}}\left]_+\mathrm{ln}\right({\displaystyle \frac{\mu ^2}{m^2}}\left)\right\{C_A^2(2\mathrm{ln}\left({\displaystyle \frac{t_1}{u_1}}\right)`$
$`2\mathrm{R}\mathrm{e}L_\beta +2\mathrm{ln}\left({\displaystyle \frac{u_1}{m^2}}\right)+\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{m^2}}\right))+2C_A(b_2+2C_F(\mathrm{Re}L_\beta +1))\}`$
$`+\left[{\displaystyle \frac{1}{s_4}}\right]_+\mathrm{ln}^2\left({\displaystyle \frac{\mu ^2}{m^2}}\right)\left\{C_A^2\mathrm{ln}\left({\displaystyle \frac{u_1}{m^2}}\right){\displaystyle \frac{1}{2}}C_Ab_2\right\},`$
with $`\mu `$ the $`\overline{\mathrm{MS}}`$-mass factorization scale, $`b_2=(11C_A2n_f)/12`$ and $`L_\beta `$ given below Eq.(17). We have checked that to NLL accuracy, the result in Eq.(3) agrees with the exact $`𝒪(\alpha _s)`$ corrections calculated in Refs..

Next, we perform a numerical investigation of the results (3) and (3). To that end, we relate the total partonic cross section difference to dimensionless coefficient functions $`\mathrm{\Delta }c_g^{(k,l)}`$,

$`\mathrm{\Delta }\sigma _{\gamma ^{}g}(s^{},q^2,m^2)`$ $`=`$ $`{\displaystyle \frac{\alpha \alpha _se_\mathrm{q}^2}{m^2}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(4\pi \alpha _s(\mu ))^k{\displaystyle \underset{l=0}{\overset{k}{}}}\mathrm{\Delta }c_g^{(k,l)}(\eta ,\xi )\mathrm{ln}^l{\displaystyle \frac{\mu ^2}{m^2}}.`$ (22)
As completely inclusive quantities the functions $`\mathrm{\Delta }c_g^{(k,l)}`$ depend only on the scaling variables

$`\eta `$ $`=`$ $`{\displaystyle \frac{s}{4m^2}}1,\xi ={\displaystyle \frac{Q^2}{m^2}},`$ (23)
where $`\eta `$ is a measure of the distance to the partonic threshold.

In Fig.1 we plot all coefficient functions $`\mathrm{\Delta }c_g^{(k,0)};k=0,1,2`$, i.e. those not accompanied by scale-dependent logarithms for $`Q^2=10\mathrm{GeV}^2`$ and $`m=1.5\mathrm{GeV}`$. Only the Born function $`\mathrm{\Delta }c_g^{(0,0)}`$ in Fig.1a is known exactly . For $`\mathrm{\Delta }c_g^{(1,0)}`$ in Fig.1a and $`\mathrm{\Delta }c_g^{(2,0)}`$ in Fig.1b we can give estimates to LL and NLL accuracy. In the LL case, we keep only the $`[\mathrm{ln}(s_4/m^2)/s_4]_+`$ and $`[\mathrm{ln}^3(s_4/m^2)/s_4]_+`$ terms in Eqs.(3) and (3) respectively. We find both $`\mathrm{\Delta }c_g^{(1,0)}`$ and $`\mathrm{\Delta }c_g^{(2,0)}`$ to be sizable near threshold where the large logarithms dominate, while they tend to be numerically rather small at larger values of $`\eta `$. Moreover, we also investigated some next-to-next-to leading (NNLL) logarithmic contributions333 For complete resummation to NNLL accuracy one would need to match the resummed result Eq.(2) at NLO, which requires knowledge of all one-loop soft and virtual corrections. These are not yet available ., such as the Coulomb corrections  for $`\mathrm{\Delta }c_g^{(1,0)}`$ and for $`\mathrm{\Delta }c_g^{(2,0)}`$ those NNLL terms, which we obtain from the expansion of the resummed result Eq.(2). Numerically, we found these NNLL terms to have an effect on $`\mathrm{\Delta }c_g^{(1,0)}`$ or $`\mathrm{\Delta }c_g^{(2,0)}`$ of the order of $`5\%`$ as compared to our NLL corrections.

Since no exact results for $`\mathrm{\Delta }c_g^{(1,0)}`$ and $`\mathrm{\Delta }c_g^{(2,0)}`$ are available we are unable to judge the quality of our approximation. However, similar investigations in the case of unpolarized heavy quark production , where exact one-loop results are known , suggest that generally NLL accuracy provides a very good approximation of the true result for $`Q^2`$ not too large444In the regime of large $`Q^2/m^2`$ on the other hand, the coefficient functions $`\mathrm{\Delta }c_g^{(k,0)}`$ can be approximated with different methods based on the operator product expansion ., because $`K^{(1)}`$, $`K^{(2)}`$ are the same as for the unpolarized case, only $`B_{\gamma ^{}g}^{\mathrm{Born}}`$ differs from $`\mathrm{\Delta }B_{\gamma ^{}g}^{\mathrm{Born}}`$.

Let us now turn to those coefficient functions multiplying scale-dependent logarithms. These are $`\mathrm{\Delta }c_g^{(1,1)}`$ in Fig.2a and $`\mathrm{\Delta }c_g^{(2,1)}`$, $`\mathrm{\Delta }c_g^{(2,2)}`$ in Fig.2b, which we plot for the same values of parameters as in Fig.1. We find in particular the estimates for $`\mathrm{\Delta }c_g^{(1,1)}`$ and $`\mathrm{\Delta }c_g^{(2,1)}`$ to NLL accuracy to be large near threshold. Additionally, we are able to obtain the exact result for the function $`\mathrm{\Delta }c_g^{(1,1)}`$ by means of renormalization group methods,

$`\mathrm{\Delta }c_g^{(1,1)}(\eta (x),\xi )`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle \underset{ax}{\overset{1}{}}}𝑑z\left(b_2\delta (1z){\displaystyle \frac{1}{2}}\mathrm{\Delta }P_{gg}^{(0)}(z)\right)\mathrm{\Delta }c_g^{(0,0)}(\eta \left({\displaystyle \frac{x}{z}}\right),\xi ),`$ (24)
with $`a=1+4m^2/Q^2`$ and the one-loop splitting function $`\mathrm{\Delta }P_{gg}^{(0)}`$ given in Refs. . Fig.2a clearly shows that the approximation based on NLL accuracy traces the exact result (24) extremely well even at larger $`\eta `$. For $`\mathrm{\Delta }c_g^{(2,1)}`$ we have checked that NNLL terms obtained from the expansion of the resummed result Eq.(2) have only a small effect.

Let us now investigate the impact of the approximate higher order perturbative corrections on the inclusive hadronic structure function $`g_1`$. In terms of coefficient functions, the charm structure function $`g_1^\mathrm{c}`$ can be expanded as follows:

$`g_1^\mathrm{c}(x,Q^2,m^2)`$
$`={\displaystyle \frac{\alpha _s\left(\mu \right)e_\mathrm{c}^2Q^2}{8\pi ^2m^2x}}{\displaystyle \underset{ax}{\overset{1}{}}}𝑑z\mathrm{\Delta }\varphi _{g/P}(z,\mu ^2){\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(4\pi \alpha _s\left(\mu \right))^k{\displaystyle \underset{l=0}{\overset{k}{}}}\mathrm{\Delta }c_g^{(k,l)}(x/z,\xi )\mathrm{ln}^l{\displaystyle \frac{\mu ^2}{m^2}}`$
$`={\displaystyle \frac{\alpha _s\left(\mu \right)e_\mathrm{c}^2}{2\pi ^2}}{\displaystyle \underset{\mathrm{}}{\overset{A}{}}}d\left(\mathrm{log}_{10}\eta \right)\mathrm{ln}10\eta \mathrm{\Delta }\varphi _{g/P}(\eta ,x,\mu ^2){\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(4\pi \alpha _s\left(\mu \right))^k{\displaystyle \underset{l=0}{\overset{k}{}}}\mathrm{\Delta }c_g^{(k,l)}(\eta ,\xi )\mathrm{ln}^l{\displaystyle \frac{\mu ^2}{m^2}},`$
where $`a=1+4m^2/Q^2`$, $`A=\mathrm{log}_{10}(\xi \{1/x1\}/41)`$, and $`\mathrm{\Delta }\varphi _{g/P}`$ represents the polarized gluon distribution. We recall that the $`\overline{\mathrm{MS}}`$-scheme has been chosen and that contributions from light initial state quarks are neglected. For $`\mathrm{\Delta }\varphi _{g/P}`$ we compare the parametrizations of Refs.. For our analysis at NLO (LO) we use a 2-loop (1-loop) running coupling with $`n_f=3`$ light flavours, a charm (pole) mass of $`m=1.5\mathrm{GeV}`$ , and $`\mathrm{\Lambda }_{\mathrm{QCD}}=0.232\mathrm{GeV}`$.

In Eq.(3), we have chosen to rewrite the standard expression for the convolution of $`g_1^\mathrm{c}`$ in a form, that facilliates the investigation of partonic threshold effects on $`g_1^\mathrm{c}`$. To do so, we plot $`\eta \mathrm{\Delta }\varphi _{g/P}`$ as a function of $`\eta `$ in Fig.3, and examine the support it gives to the coefficient functions in the second expression of Eq.(3). By comparing Fig.3 with Figs.1 and 2 one can directly judge over which ranges of $`\eta `$ the function $`\eta \mathrm{\Delta }\varphi _{g/P}`$ becomes large and hence samples the partonic coefficient functions.

Fig.3 reveals that for all parametrizations this happens indeed in the threshold region for $`x\text{ }>\mathrm{\hspace{0.17em}0.01}`$ and at scales around $`Q^210\mathrm{GeV}^2`$. We therefore expect our estimates for $`\mathrm{\Delta }c_g^{(k,l)}`$ to provide a good description of the true higher order corrections for $`g_1^\mathrm{c}`$ for $`x\text{ }>\mathrm{\hspace{0.17em}0.01}`$. The gluon densities from set A of Ref. and the valence scenario of Ref. are both positive and similar in shape, while the gluon density from set C of Ref. relaxes the positivity constraint on $`\mathrm{\Delta }\varphi _{g/P}`$ and oscillates.555Recent results from the HERMES collaboration  indicate a positive ratio of polarized over unpolarized gluon distribution $`\mathrm{\Delta }\varphi _{g/P}/\varphi _{g/P}`$ at intermediate $`x`$.

Next, we investigate the $`x`$-dependence of $`g_1^\mathrm{c}`$. In particular we are interested in a comparision of our NLO and NNLO estimates (3) and (3) with the leading order result . In Fig.4a we plot the leading order result for $`xg_1^\mathrm{c}`$ as well as the NLO and the NNLO approximations over a range $`0.007x1`$. We choose $`Q^2=10\mathrm{GeV}^2`$, a fixed value of the factorization scale $`\mu =m=1.5\mathrm{GeV}`$ and the gluon parametrization GS A of Ref.. We obtain similar results with the gluon distribution in the valence scenario of Ref.. We find that the perturbative corrections are sizable, both in the region of intermediate $`x`$, around $`0.05`$ and at smaller $`x`$, where however, the approximation is less certain to work well. To assess the quality of our approximation, we compare at each order our results to LL and to NLL accuracy. At intermediate $`x`$, the small differences between LL and NLL accuracy show a very good stability of the threshold approximation for the description of $`g_1^\mathrm{c}`$. Towards smaller $`x`$ however, these deviations increase and indicate that our approximate methods fail for $`x0.01`$.

In Fig.4b we repeat the plot of $`xg_1^\mathrm{c}`$ for the same parameters as in Fig.4a but with the gluon density GS C of Ref.. Again, the perturbative corrections are sizable over the whole range in $`x`$, but the shape of the curves is completely different. Both the coefficient functions and the gluon density oscillate, which leads in particular at intermediate $`x`$ around $`0.1`$ to destructive interference, with $`g_1^\mathrm{c}`$ being only marginally different from zero.

The final issue we study is the dependence of $`g_1^\mathrm{c}`$ on the mass factorization scale. In general, leading order calculations exhibit a strong sensitivity on the factorization scale, which is usually reduced once higher order corrections are taken into account. In addition, there are general arguments supporting a reduction in scale dependence from including soft gluon effects . Therefore, we are interested in the effect of our NLL-approximate NLO and NNLO results (3) and (3) on the scale dependence of $`g_1^\mathrm{c}`$ in comparision to the leading order calculation . Note however, that the arguments of Ref. leading to genuine NLO soft gluon resummations rely on NNLL accuracy, which is not yet available for $`g_1^\mathrm{c}`$, see footnote 3.

We plot the ratio $`g_1^\mathrm{c}\left(\mu \right)/g_1^\mathrm{c}\left(\mu =m\right)`$ over a range $`m\mu 10\mathrm{G}\mathrm{e}\mathrm{V}`$ and fix $`m=1.5\mathrm{GeV}`$, $`Q^2=10\mathrm{GeV}^2`$ and $`x=0.05`$. In Fig.5a we use the parametrization GS A of Ref., but the following conclusions hold also for the gluon density in the valence scenario of Ref.. We find at NLO that the corrections based on LL accuracy only are not sufficient to reduce the scale dependence. Clearly, soft gluon effects need to be approximated at least to NLL accuracy to achieve the desired result, a feature that has also been noticed in studies of unpolarized heavy quark production . Our best estimate for $`g_1^\mathrm{c}`$ makes use of the exact expression for the coefficient function $`\mathrm{\Delta }c^{(1,1)}`$ of Eq.(24). However, since it differs only slightly from the NLO result to NLL accuracy for the chosen parameters in Fig.5, we do not display it here.

In Fig.5b we repeat the analysis for the parametrization GS C of Ref.. In this case, the NLO and NNLO approximations do not reduce the scale dependence. This is due to the oscillating gluon density, which leads to $`g_1^\mathrm{c}`$ being close to zero at the $`x`$-value chosen, and even causes $`g_1^\mathrm{c}`$ to change sign, depending on the scale.

To summarize, the study of the charm structure function as in Figs.4 and 5 shows that the soft gluon estimates of higher order corrections to the coefficient functions are well under control and give stable predictions for $`g_1^\mathrm{c}`$ at scales $`Q^2`$ not too large and $`x\text{ }>\mathrm{\hspace{0.17em}0.01}`$. On the other hand, in the chosen kinematical range the effects of higher orders do not upset the sensitivity to the polarized gluon distribution function, with different gluon parametrizations leading to qualitatively different behaviour for $`g_1^\mathrm{c}`$. Therefore, measurements of the charm structure function allow to further constrain the polarized gluon density. Finally, we found that the sizable variation of $`g_1^\mathrm{c}`$ at leading order due to different values of the charm mass are not reduced by our approximate higher order corrections.

## 4 Conclusions

In this letter we have investigated the effects of soft gluons in the case of polarized leptoproduction of heavy quarks. Our study helps to keep theoretical control to next-to-leading logarithmic accuracy over large higher order corrections due to threshold logarithms and we have presented analytical results for the all-order resummed cross section in single-particle inclusive kinematics.

Moreover, we have provided both analytical and numerical results for the approximate NLO and NNLO perturbative corrections to the deep-inelastic charm structure function $`g_1^\mathrm{c}`$, to NLL accuracy. As a consequence of the scattering reaction being largely driven by initial state gluons with not much more energy than required to produce the final state, we have found our analysis of $`g_1^\mathrm{c}`$ to be well applicable in the kinematical range accessible to the HERMES and COMPASS experiments. In this region, our NLO and NNLO estimates can help to reduce theoretical uncertainties and may assist in the theoretical interpretation of future $`g_1^\mathrm{c}`$ measurements.

### Acknowledgments

We would like to thank E. Laenen, P. Mulders and J. Smith for fruitful discussions. We are also grateful to E. Laenen for comments on the manuscript and to J. Smith for providing us with Ref..

This work is part of the research programme of the Foundation for Fundamental Research of Matter (FOM) and the National Organization for Scientific Research (NWO).


format: text

text:
# A cannonball model of gamma-ray bursts: superluminal signatures

## 1 Introduction

The rapid directional localization of gamma-ray bursts (GRBs) by the satellites BeppoSAX (e.g. Costa et al. 1997), Rossi (e.g. Levine et al. 1996) and by the IPN (Inter-Planetary Network) of spacecrafts (e.g. Cline et al. 1999) led to the discovery of long-duration GRB afterglows, to determinations of their redshifts that verified their cosmological origin (Meegan et al. 1992), to the identification of their birthplaces —mainly star formation regions in normal galaxies— and to evidence of a possible association of GRBs with supernova explosions (Galama et al. 1998; Kulkarni et al. 1998a). The enormous isotropic energies inferred from the redshifts and fluences and the short-time variability have demonstrated that GRBs are produced by highly collimated emissions from gravitational stellar collapse (Shaviv and Dar 1995; Dar 1998a and references therein), but the specific progenitors and production mechanisms are unknown. In Table I we summarize some properties of the localized GRBs, for a currently favoured $`\mathrm{\Omega }=1`$ Universe with cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, matter density $`\mathrm{\Omega }_\mathrm{M}=0.3`$ (all in critical units), and Hubble constant $`\mathrm{H}_0=65`$ km/(s Mpc).

The prevalent widespread belief is that GRBs are generated by synchrotron emission from fireballs or firecones produced by mergers of compact stars (Paczynski 1986, Goodman et al. 1987) or by hypernova explosions (Paczynski 1998). However, various observations suggest that most GRBs are produced in supernova events by highly collimated superluminal jets (for more details, see Dar and De Rújula 2000a).

In this paper, we outline a relativistic-cannonball model of GRBs, which, we contend, explains the main observed features of GRBs and their afterglows. We concentrate here on a major prediction of the model: the superluminal motion of the GRBs afterglow. This specific signature, as we shall see, can still be searched for —with a measure of urgency— in the very particular case of SN1998bw/GRB 980425.

## 2 The supernova–GRB association

There is mounting evidence for an association of supernova (SN) explosions and GRBs. The first example was GRB 980425 (Soffita et al. 1998; Kippen 1998), within whose error circle SN1998bw was soon detected optically (Galama et al. 1998) and at radio frequencies (Kulkarni et al. 1998a). The chance probability for a spatial and temporal coincidence is less than $`10^4`$ (e.g. Galama et al. 1998), or much smaller if the revised BeppoSAX position (e.g. Pian, 1999) is used in the estimate. The unusual radio (Kulkarni et al. 1998a; Wieringa et al. 1999)) and optical (Galama et al. 1998; Iwamoto et al. 1998) properties of SN1998bw, which may have been blended with the afterglow of GRB 980425, support this association. The exceptionally small fluence and redshift of GRB 980425 make this event very peculiar, a fact that we discuss in detail below.

Evidence for a SN1998bw-like contribution to a GRB afterglow (Dar 1999a) was first found by Bloom et al. (1999) for GRB 980326, but the unknown redshift prevented a quantitative analysis. The afterglow of GRB 970228 (located at redshift $`\mathrm{z}=0.695`$) appears to be overtaken by a light curve akin to that of SN1998bw (located at $`\mathrm{z}_{\mathrm{bw}}=0.0085`$), when properly scaled by their differing redshifts (Dar 1999b). Let the energy flux density of SN1998bw be $`\mathrm{F}_{\mathrm{bw}}[\nu ,\mathrm{t}]`$. For a similar supernova located at z:

$`\mathrm{F}[\nu ,\mathrm{t}]=`$ $`{\displaystyle \frac{1+\mathrm{z}}{1+\mathrm{z}_{\mathrm{bw}}}}{\displaystyle \frac{\mathrm{D}_\mathrm{L}^2(\mathrm{z}_{\mathrm{bw}})}{\mathrm{D}_\mathrm{L}^2(\mathrm{z})}}\times `$ (1)
$`\mathrm{F}_{\mathrm{bw}}[\nu {\displaystyle \frac{1+\mathrm{z}}{1+\mathrm{z}_{\mathrm{bw}}}},\mathrm{t}{\displaystyle \frac{1+\mathrm{z}_{\mathrm{bw}}}{1+\mathrm{z}}}]\mathrm{A}(\nu ,\mathrm{z}),`$
where $`\mathrm{A}(\nu ,\mathrm{z})`$ is the extinction along the line of sight. The SN–GRB association in the case of GRB 970228 was reconfirmed by Reichart (1999) and by Galama et al. (2000). Evidence of similar associations is found for GRB 990712 (Hjorth et al. 1999; Sahu et al. 2000), GRB 980703 (Holland 2000; Dar and De Rújula 2000a) and GRB 000418: an example that we show in Fig. 1. In the case of GRB 990510, the evidence is marginal (Dar and De Rújula 2000a). For the remaining cases in Table I the data preclude a conclusion, for one or more reasons: the late afterglow is not measured; $`\mathrm{F}_{\mathrm{bw}}[\nu ^{}]`$ is not known for large $`\nu ^{}\nu (1+\mathrm{z})`$; the GRB’s afterglow or the host galaxy are much brighter than the SN. The case of GRB 970508, for which the afterglow is brighter than a SN contribution given by Eq. (1), is shown in Fig. 2.

All in all, it is quite possible that a good fraction of GRBs be associated with core-collapse SNe, perhaps even all of the most frequent, long-duration GRBs. The converse statement —that most SNe of certain types are associated with GRBs— appears at first sight to be untenable. The rate of Type Ib/Ic/II SNe has been estimated at $`\mathrm{R}_{\mathrm{SN}}1/50\mathrm{y}^1`$ per $`\mathrm{L}_{}`$ galaxy (e.g. Van den Bergh & Tammann 1991) while the rate of GRBs (for our adopted cosmological parameters) is a mere $`\mathrm{R}_{\mathrm{GRB}}5\times 10^9\mathrm{y}^1`$ per $`\mathrm{L}_{}`$ galaxy. Thus, very few of these SNe produce visible GRBs. But, if the SN-associated GRBs were beamed within an angle $`\theta 10^3`$, only a fraction $`\pi \theta ^2/4\pi 2.5\times 10^7`$ would be visible, making the observed rates compatible and making possible a rough one-to-one SN–GRB association (or a ten-to-one association for $`\theta 3\times 10^3`$). Beaming within such narrow cones also solves the GRB “energy crisis” (Dar and Plaga 1999).

To proceed, we need some understanding of the generation of superluminal jets in SNe. Lacking definite theoretical guidance, we turn to observations in search for enlightment.

## 3 Superluminal jets

Relativistic jets seem to be emitted by all astrophysical systems wherein mass is accreted at a high rate from a disk onto a central compact object (for a recent review, see Mirabel and Rodriguez 1999a). Highly relativistic jets have been observed with high resolution in galactic superluminal sources, such as the microquasars GRS 1915+105 (Mirabel and Rodriguez 1994, 1999a,b; Rodriguez and Mirabel 1999) and GRO J165-40 (Tingay et al. 1995) where mass is accreted onto a stellar black hole, and in many active galactic nuclei (AGNs) hosting a massive black hole. These jets are not continuous streams: they consist of individual blobs of plasma (plasmoids or “cannonballs”). The firing of these cannonballs (CBs) coincides with a sudden removal of the accretion-disk material (Mirabel and Rodriguez 1999b).

As they travel, microquasar CBs expand at the speed of a relativistic plasma (c/$`\sqrt{3}`$ in their rest system) presumably due to the pressure of their enclosed radiation. As they become transparent and cool down, their lateral size stabilizes to a roughly constant value, presumably constrained by magnetic self-containment and/or by the ram pressure of the ambient material. Quasar CBs show no measurable expansion as they travel, sometimes for as long as Mpcs. Galactic and quasar CBs expand explosively when finally stopped by the material they traverse.

Microquasars CBs are emitted in pairs, at similar velocities, $`\mathrm{v}=\beta \mathrm{c}`$, and angles $`\theta `$ and $`\pi \theta `$ relative to the line of sight. Their apparent transverse velocities are $`\mathrm{v}_\mathrm{T}^\pm =\beta \mathrm{c}\mathrm{sin}\theta /(1\beta \mathrm{cos}\theta )`$. For $`\beta (\mathrm{sin}\theta +\mathrm{cos}\theta )>1`$ the approaching CB is superluminal: $`\mathrm{v}_\mathrm{T}^+>\mathrm{c}`$, a case of “relativistic aberration” (Rees 1966). For $`\theta 1`$ and $`\gamma 1`$:

$$\mathrm{v}_\mathrm{T}^+\frac{2\gamma ^2\theta }{(1+\gamma ^2\theta ^2)}\mathrm{c},$$
(2)
yielding superluminal velocities $`\mathrm{v}_\mathrm{T}2\mathrm{c}/\theta `$ for $`\gamma ^2\theta ^21`$, and $`\mathrm{v}_\mathrm{T}2\mathrm{c}\gamma ^2\theta `$, for $`\gamma ^2\theta ^21`$, provided $`2\gamma ^2\theta >1`$. Superluminal velocities of a few c have been observed in microquasars, and as high as 30 c in quasars (Gabudza et al. 1993).

## 4 The cannonball model of GRBs

The ejection of matter in a SN explosion is not fully understood. The known mechanisms for imparting the required kinetic energy to the ejecta are inefficient. Simulations of stellar collapse for parent stars rotating at a typical observed rate give rise to proto-neutron stars spinning much faster than young pulsars (Heger et al. 2000). Magnetic fields no doubt play an important role in the outwards transport of angular momentum, but current simulations (e.g. Spruit and Phinney, 1998) cannot fully cope with their effects. In the related problem of accretion disks, a freely parametrized bulk viscosity must be introduced to describe the efficient angular momentum transport, presumably due to turbulent magnetohydrodynamic effects (e.g. Longair, 1994). For all these reasons, our theoretical understanding of core-collapse SN events is still unsatisfying.

We posit that the result of a SN event is not just a compact object plus the ejecta: a fraction of the parent star falls back onto the newly born compact object (De Rújula 1987). The free-fall time from a parent stellar radius $`\mathrm{R}_{}`$ into an object of mass $`\mathrm{M}_\mathrm{c}`$ is:

$`\mathrm{t}_{\mathrm{fall}}`$ $`=\pi \left[{\displaystyle \frac{\mathrm{R}_{}^3}{8\mathrm{G}\mathrm{M}_\mathrm{c}}}\right]^{1/2}`$ (3)
$`1\mathrm{day}\left[{\displaystyle \frac{\mathrm{R}_{}}{10^{12}\mathrm{cm}}}\right]^{3/2}\left[{\displaystyle \frac{1.4\mathrm{M}_{}}{\mathrm{M}_\mathrm{c}}}\right]^{1/2}.`$
For a parent star that rotates, this time would be longer, but for the material that falls from the polar regions. It is quite natural to suppose that the infalling material settles into an orbiting disk, or a thick torus if its mass is comparable to $`\mathrm{M}_\mathrm{c}`$. We assume that, as in other cases of significant accretion into a compact object (microquasars and AGNs), jets of relativistic CBs of plasma are ejected. We presume their composition to be “baryonic”, as it is in the jets of SS 433, from which Lyα and Fe Kα lines have been detected (Margon 1984), although the violence of the relativistic jetting-process may break many nuclei into their constituents.

In our admittedly sketchy model of the GRB engine, the time structure of the $`\gamma `$-ray signal will be essentially that of the ejections. The energies of the GRB subpulses and the “waiting times” between them are nearly proportional (e.g. van Paradijs et al. 1988). This is indicative of a process, such as discontinuous accretion, in which the “reservoir” must be replenished.

Let “jet” stand for the ensemble of CBs emitted in one direction in a SN event. If a momentum imbalance between the opposite-direction jets is responsible for the large peculiar velocities $`\mathrm{v}_{\mathrm{NS}}450\pm 90\mathrm{km}\mathrm{s}^1`$ (Lyne and Lorimer 1994) of neutron stars born in SNe, the jet kinetic energy must be, as we shall assume for our GRB engine, $`10^{52}`$ erg (e.g. Dar and Plaga 1999). For this choice the jet mass is relatively small: $`\mathrm{M}_{\mathrm{jet}}1.5\times 10^6\mathrm{M}_{\mathrm{NS}}(10^3/\gamma )`$, or $`1.86\mathrm{M}_{}`$ for the Lorentz factor $`\gamma =𝒪(10^3)`$ that we independently estimate from the GRBs beaming angle, from the average fluence per pulse observed by BATSE in GRBs, and from the typical duration of these pulses.

The large peculiar velocities of pulsars may have been acquired well after their birth, in a “second” transition of the cooled and spun-down neutron star to a strange star or a quark star (Dar and De Rújula 2000b). In that case $`\mathrm{v}_{\mathrm{NS}}100`$ km s-1 ought to be used: the observed velocity of young pulsars, millisecond pulsars and stellar black holes in binaries and isolation (e.g. Toscano et al. 1999 and references therein). The typical jet energy would then be $`10^{51}`$ erg and our results would have to be accordingly scaled, but their dependence on $`\mathrm{E}_{\mathrm{jet}}`$ is weak and the conclusions are unaffected.

There are other events in which a variety of GRBs could be produced by mechanisms similar to the ones we have discussed: large mass accretion episodes in binaries including a compact object, mergers of neutron stars with neutron stars or black holes, the “second transitions” just discussed, etc. In each case, the ejected cannonballs would make GRBs by hitting stellar winds or envelopes, circumstellar mass or light. We discuss only core-collapse SN explosions, as the GRBs they would produce by our mechanism, although relatively “standard”, satisfactorily reproduce the general properties of the heterogeneous ensemble of GRBs.

## 5 The making of a GRB

The making of a GRB is illustrated in Fig. 3. In the time after core collapse $`\mathrm{t}_{\mathrm{fall}}=𝒪(1)`$ day, as in Eq. (3), that it takes to “fuel and start” the GRB engine, the SN shell will have acquired a radius $`\mathrm{R}_\mathrm{s}2.6\times 10^{14}`$ cm, for a typical velocity $`\mathrm{c}/10`$ of the ejecta (we do not scale all results with all the ratios of actual to typical values, as we are only sketching the model). If the CBs expand with a velocity $`\beta _{\mathrm{exp}}\mathrm{c}`$ in their rest system ($`\beta _{\mathrm{exp}}=1/\sqrt{3}`$ for a relativistic plasma), by the time they hit the SN shell their transverse radius has grown to $`\beta _{\mathrm{exp}}\mathrm{R}_\mathrm{s}/\gamma `$, or $`\mathrm{R}_{\mathrm{CB}}1.5\times 10^{11}\mathrm{cm}`$ for $`\beta _{exp}=1/\sqrt{3}`$ and $`\gamma 10^3`$.

The typical number of CBs per SN event we certainly cannot predict. But the observed number of bursts in intense GRBs averages to $`\mathrm{n}_{\mathrm{CB}}5`$: the number —in our model— of well aimed and relatively energetic (i.e. observable) GRBs per SN. Let a CB have a fraction f of the total jet energy and let $`\mathrm{M}_\mathrm{s}`$ be the mass of the SN shell. The ratio of the swept-up shell mass to the CB’s energy is $`\mathrm{M}_\mathrm{s}\mathrm{c}^2(\mathrm{R}_{\mathrm{CB}}/\mathrm{R}_\mathrm{s})^2/(4\mathrm{f}\mathrm{E}_{\mathrm{jet}})`$, that is $`6\times 10^4/\mathrm{f}`$ for $`\mathrm{M}_\mathrm{s}=10\mathrm{M}_{}`$ and $`\mathrm{E}_{\mathrm{jet}}=10^{52}`$ erg. The ensemble of CBs (for which $`\mathrm{f}=1`$) certainly blows a passage through the SN shell, but, for f sufficiently small, the leading CBs are significantly decelerated by the shell, opening up the possibility of CBs colliding and merging after they exit the shell, which complicates the nature and time structure of the GRB signal. The fraction of the SN shell surface that is blown up by a CB is only $`5\times 10^7(\beta _{\mathrm{exp}}/0.1)^2(10^3/\gamma )^2`$, so that it appears unlikely that the CBs are aimed precisely enough to hit the very same place: CB pile-up may be the exception, not the rule.

The collisions of the CB nuclei —whose energy is $`𝒪`$(1) TeV— with those in the SN shell produce hadronic cascades. The SN shell and the CB are transparent to the neutrinos, semitransparent to the muons and opaque to the photons constituting these cascades. The CB inpact on the shell generates a very intense, narrowly collimated and very short burst of sub-TeV neutrinos. The fraction of c.m. energy that does not escape in neutrinos or muons ($`1/3`$) is converted into internal black-body radiation of the CB. The CBs become visible in $`\gamma `$-rays when the obscuring column density of the shell approaches one absorption length, $`\mathrm{X}_\mathrm{s}(\mathrm{r})\mathrm{X}_\gamma 10\mathrm{g}\mathrm{cm}^2`$, resulting in a short burst of sub-TeV $`\gamma `$-rays from $`\pi ^0`$ decays. As a reheated CB expands, it becomes transparent and emits its internal radiation, collimated and boosted by the CB’s motion and observed in the form of MeV $`\gamma `$-rays: one of the pulses that constitute a GRB. Of all these issues (Dar and De Rújula 2000a), we only discuss here the total energy in GRBs and the time duration of their pulses.

Approximately a third of the collision c.m. energy is converted in the hadronic cascades into electromagnetic radiation (via $`\pi ^0`$ production). Energy and momentum conservation in the collision between a CB and the SN shell imply that the internal heat acquired by the CB is $`\mathrm{E}_\mathrm{i}(\beta _{\mathrm{exp}}/6\gamma )\sqrt{2\mathrm{f}\mathrm{E}_{\mathrm{jet}}\mathrm{M}_\mathrm{s}\mathrm{c}^2}`$. The corresponding radiation pressure makes the CB expand and cool at a rate $`\mathrm{T}_{\mathrm{CB}}1/\mathrm{R}_{\mathrm{CB}}1/\mathrm{t}`$, until it becomes optically thin and releases its enclosed radiation. For an ionized CB, this takes place when $`\mathrm{R}_{\mathrm{CB}}[3\mathrm{M}_{\mathrm{CB}}\sigma _\mathrm{T}/(4\pi \mathrm{m}_\mathrm{p})]^{1/2}`$, where $`\sigma _\mathrm{T}0.65\times 10^{24}`$ cm2 is the Thomson cross section111Nuclei are disassembled by the CB’s collision with the SN shell, and neutrons decay in $`1`$ s of observer’s time, thus the number of electrons in the CB is $`\mathrm{M}_{\mathrm{CB}}/\mathrm{m}_\mathrm{p}`$..

The internal heat, reduced by expansion, is the energy the CB radiates in its rest system, which is boosted and collimated by the CB’s motion. Local observers at angles $`\theta `$ to the CB’s direction would see this energy distributed as:

$$\frac{\mathrm{dE}}{\mathrm{d}\mathrm{\Omega }}\frac{\beta _{\mathrm{exp}}^2}{\gamma ^{3/2}}\left[\frac{2\gamma }{(1+\gamma ^2\theta ^2)}\right]^3\left[\frac{4\pi \mathrm{R}_\mathrm{s}^2}{3\sigma __\mathrm{T}}\right]^{1/2}\frac{\sqrt{2\mathrm{m}_\mathrm{p}\mathrm{M}_\mathrm{s}}\mathrm{c}^2}{24\pi }.$$
(4)
Remarkably, the CB’s mass and energy have dropped from the above expression, except for the fact that, for the result to be correct, they must be large enough for the CB to pierce the SN shell and remain relativistic. For a most likely viewing angle $`\theta 1/\gamma `$, and $`\gamma =10^3`$, one infers from Eq. (4) an isotropic energy $`\mathrm{E}_{\mathrm{isot}}(\beta _{\mathrm{exp}}/0.1)^2\mathrm{\hspace{0.17em}6.4}\times 10^{53}\mathrm{erg}`$. This energy —or $`\mathrm{n}_{\mathrm{CB}}`$ times this energy in a GRB with that many pulses— are in agreement with the range of observed isotropic energies listed in Table I. The light from the SN shell is also Compton up-scattered to MeV energies, but its contribution to a GRB is subdominant relative to that of Eq. (4) (Dar and De Rújula 2000a).

The duration of a single CB-induced GRB pulse is the light crossing time of the CB when it becomes transparent to its radiation, i.e. $`\mathrm{R}_{\mathrm{cb}}/\mathrm{c}`$ in its rest frame, or

$$\mathrm{t}\left[\frac{3\mathrm{M}_{\mathrm{CB}}\sigma _\mathrm{T}}{4\pi \mathrm{m}_\mathrm{p}}\right]^{1/2}\frac{1+\theta ^2\gamma ^2}{2\gamma \mathrm{c}}(1+\mathrm{z}).$$
(5)
for the observer. For $`z=1`$, a most probable $`\theta 1/\gamma `$, and our adopted values, $`\mathrm{E}_{\mathrm{jet}}=10^{52}`$ erg and $`\gamma =10^3`$, $`\mathrm{t}2.2\mathrm{f}^{1/2}`$ s, comparable to the observed duration of GRB pulses. The total duration of a GRB is that of the (unpredictable) firing-sequence time of CBs into the observer’s viewing cone; the duration of GRBs with only one observable pulse is given by Eq. (5).

The evidence for superluminal motion in GRBs, as we shall argue, can only be measured in their afterglow phase. The only case in which this can be done (SN1998bw/GRB 980425) is already some 850 days old: it may be dimming out of observability. Hence, we skip a more detailed discussion of the spectral properties of GRBs (Dar and De Rújula, 2000a) to concentrate on the issue of superluminal motion of the afterglow’s source.

## 6 GRB afterglows

Far from their parent SNe, the CBs are slowed down by the interstellar medium (ISM) they sweep, which has been previously ionized by the forward-beamed CB radiation (travelling essentially at $`\mathrm{v}=\mathrm{c}`$, the CB is “catching up” with this radiation, so that the ISM has no time to recombine). As in the jets and lobes of quasars, a fraction of the swept-up ionized particles are “Fermi accelerated” to cosmic-ray energies and confined to the CB by its turbulent magnetic field, maintained by the same confined cosmic rays (Dar 1998b; Dar and Plaga 1999). The synchrotron emission from the accelerated electrons, boosted by the relativistic bulk motion of the CB, produces afterglows in all bands between radio and X-rays, collimated within an angle $`1/\gamma (t)`$, that widens as the ISM decelerates the CB. When the CBs finally stop, they enter a Sedov–Taylor phase, causing a temporal “break” in their afterglow.

A CB of roughly constant cross section, moving in a previously ionized ISM of roughly constant density, would lose momentum at a roughly constant rate, independent of whether the ISM constituents are re-scattered isotropically in the CB’s rest frame, or their mass is added to that of the CB. The pace of CB slowdown is $`\mathrm{d}\gamma /\mathrm{dx}=\gamma ^2/\mathrm{x}_0`$, with $`\mathrm{x}_0=\mathrm{M}_{\mathrm{CB}}/(\pi \mathrm{R}_{\mathrm{CB}}^2\mathrm{n}\mathrm{m}_\mathrm{p})`$ and n the number density along the CB trajectory. For $`\gamma ^21`$, the relation between the length of travel dx and the (red-shifted, relativistically aberrant) time of an observer at a small angle $`\theta `$ is $`\mathrm{dx}=[2\mathrm{c}\gamma ^2/(1+\theta ^2\gamma ^2)][\mathrm{dt}/(1+\mathrm{z})]`$. Inserting this into $`\mathrm{d}\gamma /\mathrm{dx}`$ and integrating, we obtain:

$$\frac{1+3\theta ^2\gamma ^2}{3\gamma ^3}=\frac{1+3\theta ^2\gamma _0^2}{3\gamma _0^3}+\frac{2\mathrm{c}\mathrm{t}}{(1+\mathrm{z})\mathrm{x}_0},$$
(6)
where $`\gamma _0`$ is the Lorentz factor of the CB as it exits the SN shell. The real root $`\gamma =\gamma (\mathrm{t})`$ of the cubic Eq. (6) describes the CB slowdown with observer’s time.

Let the radiation emitted by a CB in its rest system (bremsstrahlung, synchrotron, Compton-boosted self-synchrotron) be isotropic and have an energy flux density $`\mathrm{F}_0(\nu _0)`$. This radiation is boosted and collimated by the CB’s motion, and its time-dependence is modified by the observer’s time flowing $`(1+z)(1+\gamma ^2\theta ^2)/(2\gamma )`$ times faster than in the CB’s rest system. For $`\gamma 1`$, an observer at small $`\theta `$ sees an energy flux density:

$`\mathrm{F}[\nu ]\left({\displaystyle \frac{\mathrm{R}_{\mathrm{CB}}}{\mathrm{D}_\mathrm{L}}}\right)^2`$ $`\left[{\displaystyle \frac{2\gamma }{1+\gamma ^2\theta ^2}}\right]^3\times `$ (7)
$`\mathrm{F}_0\left(\nu [1+\mathrm{z}]{\displaystyle \frac{1+\gamma ^2\theta ^2}{2\gamma }}\right)\mathrm{A}(\nu ,\mathrm{z}),`$
with $`\gamma =\gamma (\mathrm{t})`$ as in Eq. (6) and $`\mathrm{A}(\nu ,\mathrm{z})`$ an eventual absorption dimming.

During the afterglow regime, the CBs emission is dominated by electron synchrotron radiation from the magnetic field in the CB. At frequencies above the plasma frequency, a CB becomes transparent to this radiation when the attenuation length due to electron free–free transitions exceeds its radius. The spectral shape is $`\mathrm{F}_0\nu _0^\alpha `$, with $`\alpha =(\mathrm{p}1)/2`$ and p the spectral index of the electrons. For equilibrium between Fermi acceleration and synchrotron and Compton cooling, $`\mathrm{p}3.2`$ and $`\alpha 1.1`$, while for small cooling rates, $`\mathrm{p}2.2`$ and $`\alpha 0.6`$ (Dar and De Rújula 2000c), or $`p1.2`$ and $`\alpha 0.1`$ if Coulomb losses dominate. At very low radio frequencies self-absorption becomes important and $`\alpha 1/3(2.1)`$ for optically thin (thick) CBs. For a detailed modelling of synchrotron radiation from quasar lobes, see, for instance, Meisenheimer et al. (1989).

Since $`\gamma (\mathrm{t})`$, as in Eq. (6), is a decreasing function of time, the afterglow described by Eq. (7) may have a very interesting behaviour. An observer may initially be outside the beaming cone: $`\theta ^2\gamma ^2>1`$, as we shall argue to be the case for GRB 980425, for which we estimate $`\gamma _0^2\theta ^2100`$ (other relatively dim GRBs in Table I, such as 970228 and 970508, may also be of this type). The observed afterglow would then initially rise with time. As $`\gamma `$ decreases, the cone broadens, and at $`\gamma \theta 1`$ the time dependence changes towards $`\mathrm{F}_\nu \mathrm{t}^{(3+\alpha )/3}\mathrm{t}^{1.4}`$, while the $`\nu `$ dependence stays put at $`\nu ^\alpha `$. Such a behaviour may explain the puzzling initial rise of the optical afterglow of GRBs 970228 and 970508, as well as the second peak around $`\mathrm{t}_\mathrm{p}3`$$`4\times 10^6\mathrm{s}`$ in the unresolved radio emission from SN1998bw/GRB 980425 (Kulkarni et al. 1998a; Frail et al. 1999), if it corresponds to the GRBs afterglow.

The afterglow of GRB 970508 is shown in Fig. 2, and compared with our prediction Eqs. (6,7), for the measured index $`\alpha =1.1`$. The adjusted parameters are the height and time of the afterglow’s peak and $`\theta =10/\gamma _0`$. The figure is for a single CB; with a few of them at chosen times and relative fluxes, it would be easy to explain the early “warning shots” at $`\mathrm{t}<1`$ day and the abrupt rise at $`\mathrm{t}=1`$ to 2 days. At $`\mathrm{t}\mathrm{t}_\mathrm{p}`$, however, they would add up to a single curve like the one shown in the figure.

As a CB is finally stopped by the ISM, it enters a Sedov–Taylor phase. The CB’s radius increases as $`\mathrm{t}^{2/5}.`$ The Lorentz factor of electrons decreases like $`\mathrm{t}^{6/5}`$. In equipartition, the magnetic field decreases like $`\mathrm{t}^{3/5}`$. Wijers et al. (1997) have shown that these facts lead to an afterglow decline $`\mathrm{F}_\nu \mathrm{t}^{(15\alpha 3)/5}\mathrm{t}^{2.7}`$ for $`\alpha 1.1`$, as observed for the late-time afterglows of microquasars and of some GRBs.

## 7 GRB 980425: a special case

In the list of Table I, GRB 980425 stands out in two —apparently contradictory— ways: it is, by far, the closest ($`\mathrm{z}=0.0085`$, $`\mathrm{D}_\mathrm{L}=39`$ Mpc) and it has, by far, the smallest implied spherical energy: $`8.1\times 10^{47}`$ erg, 4 to 6 orders of magnitude smaller than that of other GRBs. If GRB 980425 was not abnormal, Eq. (4) tells us that it must have been beamed at an angle $`\theta 10/\gamma _0`$ relative to our line of sight. The probability of a jetted GRB to be pointing to us within an angle $`\theta 10/\gamma _010^2`$ is $`\mathrm{P}\pi \theta ^2/(4\pi )2.5\times 10^5`$. For this GRB, BATSE observed a single wide peak of $`5`$ s duration and a maximum flux $`(3\pm 0.3)\times 10^7\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ (Soffita et al. 1998). With the BATSE sensitivity of $`6\times 10^2\mathrm{cm}^2\mathrm{s}^1`$ in the range $`50<\mathrm{E}_\gamma <300\mathrm{keV}`$, such a burst can be detected up to $`150`$ Mpc. For our choice of $`\mathrm{H}_0`$, the local luminosity density is $`1.3\times 10^2\mathrm{L}_{}\mathrm{Mpc}^3`$ (e.g. Ellis 1997), yielding $`3.7\times 10^3`$ SNe a year within 150 Mpc. Thus, the mean time to observe a SN–GRB association with $`\theta >10/\gamma _0`$ is $``$ 10 y, consistent with a single detection, since the launch of BATSE, of a GRB with the odd properties of 980425.

## 8 Is the superluminal motion of GRB 980425 observable?

As we have discussed, the peculiarities of of GRB 980425 can be understood if its source was “fired” from SN1998bw with a bulk-motion Lorentz factor $`\gamma 1000`$, at an angle $`\theta 10/\gamma `$ relative to our line of sight. Its transverse superluminal displacement, $`\mathrm{D}_\mathrm{T}`$, from the SN position can be obtained by time-integrating $`\mathrm{v}_\mathrm{T}`$, as in Eq. (2), using $`\gamma (\mathrm{t})`$ as in Eq. (6). The result can be reproduced, to better than 10% accuracy, by using the approximation $`\mathrm{v}_\mathrm{T}2\gamma ^2\theta \mathrm{c}`$, valid for $`\gamma <1/\theta `$ (or $`\mathrm{t}>\mathrm{t}_\mathrm{p}`$, the afterglow’s peaktime):

$$\mathrm{D}_\mathrm{T}\frac{2\mathrm{c}\mathrm{t}_\mathrm{p}}{\theta }\left[\frac{\mathrm{t}}{\mathrm{t}_\mathrm{p}}\right]^{1/3}.$$
(8)
At present ($`\mathrm{t}850\mathrm{d}`$), the displacement of the GRB from the initial SN/GRB position is $`\mathrm{D}_\mathrm{T}20(\gamma _0/10^3)`$ pc, corresponding to an angular displacement $`\mathrm{\Delta }\alpha 100(\gamma _0/10^3)`$ mas, for $`\mathrm{z}=0.0085`$, $`\mathrm{t}_\mathrm{p}3.5\times 10^6`$ s.

The two sources, a few tens of mas away, may still be resolved by HST. In Fig. 4 we show our prediction for the late-time V-band light curve of SN1998bw/GRB 980425. The SN curve is a fit by Sollerman et al. (2000) for energy deposition by $`{}_{}{}^{60}\mathrm{Co}`$ decay in an optically thin SN shell. The GRB light curve is our predicted afterglow for GRB 9980425, as given by Eqs. (6,7), and constrained to peak at the position of the second peak in the radio observations (Kulkarni et al. 1998a; Frail et al. 1999). The fitted normalization is approximately that of the mean GRB afterglow of the GRBs in Table I suppressed by the same factor as its $`\gamma `$-ray fluence relative to the mean $`\gamma `$-ray fluence in Table I. The joint system has at present (day $`850`$) an extrapolated magnitude $`\mathrm{V}26`$ (Fynbo et al. 2000). It can still be resolved from its host galaxy ESO-184-G82 by HST and, perhaps, by VLT in good seeing conditions. An extrapolation of the V-band late time curve of SNR1998bw (Sollerman et al. 2000) suggests that the present magnitude of the SN is $`\mathrm{V}28`$ , which is near the detection limit of HST and is dimming much faster than $`{}_{}{}^{60}\mathrm{Co}`$-decay would imply.

In the GHz radio band, the system has been last observed in February 1999 by ATCA (Australia Telescope Compact Array) to have $``$ mJy flux density (Frail et al. 2000). It approached a power-law time decline with a power-law index -1.47. If a Sedov-Taylor break in the radio afterglow of GRB 980425 has not occurred yet, its spectral density may still be strong enough to determine its position with VLA and VLBI to better than mas precision. If the second peak is the GRB’s afterglow, the two radio centroids should now be separated by $`\gamma _0/10`$ mas (or by $`\gamma _0/20`$ mas, if the first peak is the afterglow). A refined location from a reanalysis of the early ATCA observations (Kulkarni et al 1998a; Frail et al. 1999) of the initial SN and of the late afterglow may also reveal a superluminal displacement.

If the GRB afterglow has entered the late fast-decline phase seen in some GRB afterglows and in quasar and microquasar “afterglows” from jetted ejections (observed power-law index $`2.7\pm 0.3`$), a further delay in follow-up observations can make it very difficult or impossible to detect and resolve the GRB/SNR radio image into its two predicted images.

A GRB as close as GRB 980425 (z = 0.0085) should occur only once every $`10`$ years. For typical GRBs ($`\mathrm{z}1`$) there is no hope of resolving them with HST into two separate SN and GRB images. Resolving them with VLBI would also be arduous. For these reasons, we exhort interested observers to consider immediate high-resolution optical (STIS) and radio (VLA and VLBI) follow-up observations of SN1998bw and the afterglow of GRB980425.

## 9 Conclusions

We have argued that GRBs and their afterglows may be produced by the firing of extremely relativistic cannonballs in SN explosions. For the GRB to be observable, the CBs must be close to the line of sight, implying that their afterglows would appear to move superluminally. Only one of the located GRBs (980425) is close enough to us for this superluminal displacement to be observable with the currently available resolution. Its afterglow may by now have reached a “break” and be too dim to be seen. Or it may not. If observed, the superluminal displacement of this GRB’s afterglow would be a decisive card in favour of cannonballs, as opposed to stationary fireballs.

Table I - Gamma ray bursts of known redshift z

Comments: $`a`$: Luminosity distance in Gpc (for $`\mathrm{\Omega }_\mathrm{m}=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\mathrm{H}_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. $`b`$: BATSE $`\gamma `$–ray fluences in units of $`10^5`$ erg cm-2. $`c`$: (Spherical) energy in units of $`10^{53}`$ ergs. $`d`$: R-magnitude of the host galaxy, except for GRB 990510, for which the V-magnitude is given. References: 1: Djorgovski et al. 1999; 2: Metzger et al. 1997; 3: Djorgovski et al. 2000; 4: Kulkarni et al. 1998b; 5: Tinney et al. 1998; 6: Djorgovski et al. 1998a; 7: Djorgovski et al. 1998b; 8: Kelson et al. 1999; 9: Vreeswijk et al. 1999a; 10: Galama et al. 1999; 11: Dodonov et al. 1999; 12: Vreeswijk et al. 1999b; 13: Feng et al. 2000; 14: Bloom et al. 2000; 15: Fruchter et al. 1999b; 16: Bloom et al. 1998a; 17: Djorgovski et al. 2000a; 18: Odewahn et al. 1999; 19: Galama et al. 1998; 20: Djorgovski et al. 2000b; 21: Bloom et al. 1998b; 22: Fruchter et al. 1999c; 23: Fruchter et al. 1999d; 24: Hjorth et al. 1999; 25: Diercks et al. 2000; 26: Djorgovski et al. 1999; 27: Fruchter et al. 2000; 28: Metzger et al. 2000.


text:
# Rippling Instability of a Collapsing Bubble

## Abstract

When a bubble of air rises to the top of a highly viscous liquid, it forms a dome-shaped protuberance on the free surface. Unlike a soap bubble, it bursts so slowly as to collapse under its own weight simultaneously, and folds into a striking wavy structure. This rippling effect occurs in fact for both elastic and viscous sheets, and a theory for its onset is formulated. The growth of the corrugation is governed by the competition between gravitational and bending (shearing) forces and is exhibited for a range of densities, stiffnesses (viscosities), and sizes – a result which arises less from dynamics than from geometry, suggesting a wide validity. A quantitative expression for the number of ripples is presented, together with experimental results which are in agreement with the theoretical predictions.

Every day, nature surprises us with structures and patterns of such beauty as to fill the scientist with wonder and the artist with envy. In the present paper, we address an instability which turns a hemispherical, smooth, liquid bubble into a striking wrinkled structure, first observed by Debrégeas, de Gennes, and Brochard-Wyart . In their experiment, 0.1 to 10 cm3 of air injected into a highly viscous liquid ($`\eta 10^3`$ Pa$``$s) rises to the free surface, imprisoned in a hemispherical bubble of thickness $`t110\mu `$m. If the latter is punctured at its apex by a needle, surface tension drives the rapid expansion of a circular opening. The retraction velocity soon (after about $`1030`$ ms) saturates to a constant, owing to the high viscous resistance. In the meantime, the air flow through the hole equilibrates the pressure difference, allowing the bubble to collapse under its own weight. As it deflates, an instability appears: the fluid sheet folds into a wavy structure, with radial ripples that break the original axisymmetry. In the absence of a detailed theory, Ref.  proposes a scaling estimate for the number of ripples $`n^{}\left(\mu gR^3/K\right)^{1/2}`$, where $`\mu `$ is the mass of the film per unit area, $`g`$ the gravitational acceleration, $`R`$ the radius of the hole, and $`K`$ an effective bending rigidity of the sheet which was assumed to be elastic during the early stages of the rippling.

The rippling results from the competition between compression, bending, and gravity. Each fluid element tends to fall under its own weight, but experiences a viscous resistance from its neighborhood. If the bubble were to collapse in a uniform, symmetric way, it would occupy a progressively reduced area, leading to an in-plane compression which would require forces that far exceed the scale set by gravity. Instead, the film deforms in a nearly inextensional fashion by undergoing pure bending. Equivalently, for a given (gravitational) force, the relative time scale associated with stretching is much larger than that for bending, and the surface therefore corrugates over short times, before eventually relaxing into a uniform, thicker membrane.

This instability is reminiscent of buckling phenomena , originally studied in the context of elastic rods but occurring also in the creeping flows of viscous liquid filaments (a striking everyday example being the coiling of a stream of honey when it reaches a toast ). For an elastic rod, buckling occurs at the longest possible wavelength in order to minimize the bending energy. In the bubble problem, however, gravity plays a distinctive role in determining the configuration. For a given amplitude, bending still favors large scale deformations, while gravitational energy is minimized for an almost flat sheet with as many tiny ripples as possible; the optimal wavelength results from a compromise between the two. Such an argument, however, does not fully characterize the effect. Unlike the above examples, here the system under consideration is a curved two-dimensional sheet, and the associated geometry constrains the rippling both qualitatively and quantitatively, as we shall see.

The instability occurs both an elastic (solid) film and a viscous (liquid) one. The elastic case corresponds to a convex curved shell with an apical hole of radius $`R`$ that is allowed to collapse under its own weight. In the viscous case, an additional complication arises because the radius of the hole changes during the rippling. Following a short initial transient, the hole grows steadily at a rate $`v\sigma /\eta `$ resulting from the balance of surface tension ($`\sigma 20`$mN/m) and viscous stress. It thus takes a time $`\tau \eta t/\sigma `$ for the opening radius to increase by $`t`$. During this time, the liquid acquires a velocity $`Vg\tau `$ due to gravity, larger than $`v`$ by a factor $`V/v10^7`$. Even if the liquid is viscoelastic, so that the retraction velocity is enhanced by a factor $`R/t`$ ($`1010^4`$) , the hole radius remains essentially constant while the instability occurs (cf. Fig. 1). We may therefore treat the hole radius $`R`$ as a given parameter in the theory .

Although the bubble has the geometry of a sphere before collapsing, it is quite flattened by the time the ripples appear (Fig. 1(B)). For simplicity, we consider the unperturbed configuration to be a shallow cone, of slope $`\alpha 1`$, described by its height above the surface,

$$h=\alpha \left(r_0r\right),$$
(1)
where $`r`$ is the cylindrical radial coordinate and $`r_0`$ the radius of the base. Any deformation of $`h`$ introduced by the rippling may be written, without loss of generality, as

$$h+\delta h=\alpha \left(r_0r\right)\delta \alpha \left(r\right)+\underset{n1}{}\left\{\delta \beta _n^{\left(1\right)}\left(r\right)\mathrm{cos}\left(n\theta \right)+\delta \beta _n^{\left(2\right)}\left(r\right)\mathrm{sin}\left(n\theta \right)\right\},$$
(2)
where $`\theta `$ is the azimuthal angle. The perturbation $`\delta \alpha `$ represents a uniform ($`n`$ independent) flattening accompanying the growth of ripples of amplitude $`\delta \beta _n^{\left(i\right)}`$, and a crucial step consists in understanding their form and interdependence. In the case of a thin elastic (viscous) sheet, the two primary modes of deformation are in-plane stretching (shearing) and out-of-plane bending. A generic deformation of an elastic cone (made of a material with Young modulus $`Y`$), of amplitude $`\zeta `$ on a scale $`\mathrm{}`$, requires stretching forces (per unit surface) of order $`Yt\zeta ^2/\mathrm{}^2`$, but significantly smaller stretching forces (per unit surface) of order $`Yt^3\zeta ^2/\mathrm{}^4`$ , so that for a given external force, here gravity, inextensional deformations are greatly preferred . In the case of a highly viscous sheet, forces arise from velocity gradients, thus introducing a dynamical element in the problem. However, their dependence on $`t`$ and $`\mathrm{}`$ (essentially due to the variation of the strain across the film) is similar, so that inextensional deformations are again largely favored if $`tlr_0/n^{}`$. (This condition is satisfied if the selected number of ripples $`n^{}`$ is small compared to $`10^3`$, which is the case as we shall see below.) Equivalently, for a given loading, the time scale corresponding to bending is smaller than that for stretching by a factor $`\left(t/l\right)^2`$ . Thus, at the onset of the instability, perturbations of the cone must preserve its metric. This requirement translates into the constraints

$$\delta \beta _n^{\left(i\right)}\left(r\right)=\delta \beta _n^{\left(i\right)}\times r+\delta \beta _n^{\left(i\right)},$$
(3)
(where $`\delta \beta _n^{\left(i\right)}`$ and $`\delta \beta _n^{\left(i\right)}`$ are constants) and

$$4\alpha \delta \alpha \left(r\right)=\underset{n1,i}{}\left\{\left(n^21\right)\delta \beta _n^{\left(i\right)2}\left(rr_0\right)+n^2\delta \beta _n^{\left(i\right)2}\left(\frac{1}{r}\frac{1}{r_0}\right)\right\}.$$
(4)
In the following, we elucidate the elastic (solid) case before extending our treatment to the viscous (liquid) case. The energy functional of a perturbed elastic cone reads

$`E\left[h+\delta h\right]`$ $`=`$ $`{\displaystyle _{\text{cone}}}d\left(\text{surface}\right)\times \left(\text{gravitational potential energy }+\text{ bending potential energy}\right)`$ (5)
$`=`$ $`{\displaystyle _R^{r_0}}r𝑑r{\displaystyle _0^{2\pi }}𝑑\theta \sqrt{1+\left[\left(h+\delta h\right)\right]^2}\left\{\mu g\left(h+\delta h\right)+{\displaystyle \frac{K}{2}}\left(^2\delta h\right)^2\right\}.`$ (6)
where $`K=Yt^3/12\left(1\nu ^2\right)`$ is the rigidity and $`\nu `$ the Poisson ratio. Only bending elastic energy appears in $`E`$, since we have confined ourselves to the class of inextensible deformations.

If the elastic cone is attached to the plane on which it rests so that $`\delta h\left(r=r_0\right)=0`$, Eq. (3) yields $`\delta \beta _n^{\left(i\right)}=\delta \beta _n^{\left(i\right)}r_0`$. On substituting Eqs. (24) into Eq. (6) we then obtain, to lowest order in the perturbation,

$`\delta E`$ $``$ $`E\left[h+\delta h\right]E\left[h\right]`$ (7)
$`=`$ $`{\displaystyle \frac{\pi }{2}}Kf\left({\displaystyle \frac{r_0}{R}}\right){\displaystyle \underset{n1,i}{}}\left(\delta \beta _n^{\left(i\right)}\right)^2\left\{g({\displaystyle \frac{r_0}{R}},\gamma R^3)+2\left[\phi \left({\displaystyle \frac{r_0}{R}}\right)\gamma R^3\psi \left({\displaystyle \frac{r_0}{R}}\right)\right]n^2+n^4\right\}`$ (8)
$``$ $`{\displaystyle \underset{n1}{}}\delta E_n.`$ (9)
$`\gamma ^{1/3}=\left(\alpha K/\mu g\right)^{1/3}`$ is an intrinsic length scale arising from the competition between gravity and bending elasticity .

Each mode contributes an amount $`\delta E_n`$ to the change in energy, and rippling occurs if $`\delta E_n<0`$ for some integer. In general, $`\delta E_n<0`$ for a range of different $`n`$’s; the most negative variation corresponds to the maximally growing perturbation, and thus sets the wavelength of the instability. The formulation also yields a “threshold condition” $`\phi \left(r_0/R\right)<\gamma R^3\psi \left(r_0/R\right)`$ for the occurrence of rippling. This condition involves the three independent quantities $`\gamma `$, $`r_0`$, and $`R`$, and may be translated into three corresponding statements. (i) Rippling is suppressed if $`\gamma <\gamma _\text{c}(r_0,R)=R^3\phi /\psi `$, i.e. if the cone is too light or too rigid. (ii) Similarly, no rippling occurs if the hole, or equivalently the cone, is too small, $`r_0<r_{0\text{c}}(\gamma ,R)`$. Azimuthal continuity requires the wavelength of the deformation to be at most of order $`r_0`$, resulting in a forbidding bending cost if $`r_0`$ becomes small compared to the intrinsic (energetically determined) scale $`\gamma ^{1/3}`$. (iii) Finally, the threshold depends, quite unexpectedly, on the ratio $`r_0/R`$. The dependence of the symmetric ($`n=0`$) mode on the radial coordinate $`r`$ is different from that of the rippling ($`n0`$) modes, so that the high elastic cost cannot be justified by gravitational gain anymore if the hole is reduced beyond a critical size. Minimizing $`\delta E`$ in Eq. (9) yields the selected number of ripples as

$$n^{}=\text{Int}\sqrt{\frac{\mu gR^3}{K}\frac{1}{\alpha }\psi \left(\frac{r_0}{R}\right)\phi \left(\frac{r_0}{R}\right)},$$
(10)
where Int$`\left\{x\right\}`$ is the integer closest to $`x`$. This relation improves on the estimate of Ref. , where the authors consider the short time elastic behavior, and establishes its domain of validity.

For an elastic (solid) sheet, the rippling phenomenon is of an essentially static nature; upon increasing, say, the mass of the sheet, the equilibrium configuration is shifted from symmetric to rippled. Approaching the problem from a dynamical perspective by considering the elastic forces and torques rather than the corresponding energies results in an evolution equation $`\pi \mu P(r_0,R)d^2\left(\delta \beta _n^{\left(i\right)}\right)/dt^2=\delta E/\delta \beta _n^{\left(i\right)}`$ for each mode. Here $`P`$ is a polynomial function independent of $`n`$, so that the energetically optimal mode, with number $`n^{}`$, is indeed the fastest growing one. In the case of a viscous liquid, the effect is intrinsically dynamical: bending occurs only on short times, while the equilibrium configuration is ultimately reached by a slow thickening. Nevertheless, the motion of a viscous film satisfies a formulation close to that of an elastic sheet, as can be shown by integrating the Stokes equation through the thickness . Indeed, it is easy to see that bending results from a torque $`\frac{\eta t^3}{4\left(1\nu ^2\right)}\times d(`$curvature$`)/dt`$ analogous to an elastic torque $`K\times (`$curvature$`)`$, so that a highly viscous film may be described by an effective bending modulus $`K_l=\eta t^3/3\tau `$ ($`\nu =1/2`$ for an incompressible medium), where $`\tau `$ is a time scale associated with the falling velocity. Thus, all the conclusions of the stability analysis for the elastic cone, and in particular the expression for the number of ripples (Eq. (10)), may be transposed to the case of the bubble modulo a certain time scale related to the gravity-induced velocity of the fluid. Comparing the nascent ripples’ amplitude to the film thickness yields an estimate of this time scale as $`(t/g)^{1/2}`$ .

In order to check our results against experiment, we visualized the bursting of silicone oil bubbles. Once the bubble is punctured with a sharp needle, its evolution is followed using a high-speed camera capable of recording up to 1000 frames per second. The resulting video is then analyzed to determine the radius $`r_0`$ of the bubble, the hole size $`R`$ at which the ripples are first observed, along with the number of ripples $`n^{}`$. Since the hole expands very fast at first, $`R`$ is much larger than $`R_\text{c}`$ by the time the bubble begins to collapse. To compare the experiments with the theory, in which $`R`$ enters as a parameter, the latter is measured at the onset of the instability for each given size of the bubble. The quantitative measurements are compared to the theoretical predictions for the dependence of $`n^{}`$ on the bubble size on Fig. 2. On a more qualitative level, the experiments show a suppression of the instability for small bubbles, in agreement with the threshold conditions above.

We conclude with a discussion of possible refinements of the theory and their relation to the geometric nature of the problem. A more complete theory would incorporate a (flattened) hemisphere as the initial condition, rather than a cone. Also, due to the progressive drainage of the liquid, the thickness $`t`$ acquires a dependence on $`r`$ (and time). This in turn implies non-uniform rigidity $`K\left(r\right)`$ and mass $`\mu \left(r\right)`$, leading to functions $`f`$, $`g`$, $`\phi `$, $`\psi `$, and $`P`$ of a more complicated form. On a more fundamental level, all these aspects should be addressed in terms of the coupled hydrodynamics of the slow viscous (liquid) flow and the rapid air flow; the corresponding calculations will be reported elsewhere . Yet, the strong geometrical constraints involved in the problem are suggestive of the robustness of the results. The question we have answered is akin to that of applying a curved surface unto a flat one in the most economical way, a problem which has taxed cartographers for many centuries and lies at the birth of differential geometry. It is also somewhat of an inverse counterpart to the problem of fitting a flat sheet to a three dimensional landscape, which has been studied in various contexts , and is an issue that still vexes fashion designers. The relevance of the geometrical constraints is manifest, for example, in the strong dependence of the rippling on the size of the opening, which is closely related to a well-known theorem due to Gauss , Jellett , and others, according to which (loosely put) a closed surface cannot be bent without being stretched, while an open surface can be bent inextensionally. Similarly, we find that a smaller hole implies a relatively stiffer bubble, and hampers the rippling. While the detailed form of the functions $`\phi `$ and $`\psi `$ arise from the physical constraints and dynamics imposed by the forces and various boundary conditions, the essence is in the geometry.

FIG. 1. Stroboscopic images of a collapsing liquid bubble of size is $`r_0=1\text{cm}`$ and thickness $`t100\mu \text{m}`$. The silicone oil has viscosity $`\eta =10^3\text{Pa}\text{s}`$, surface tension $`\sigma =21\text{mN/m}`$, and mass density $`0.98\text{g/cm}\text{3}`$. (A) $`30ms`$ after the film is punctured by a sharp needle, the bubble shows a retracting hole of radius $`R=1.4\text{mm}`$ at its apex, but no ripples yet. (B) $`30ms`$ later, the bubble loses its axisymmetric shape. The radius of the hole remains essentially constant, at $`R=1.6\text{mm}`$, while the ripples grow. The inset displays a schematic side view of the essentially conical deflating bubble at the onset of the instability, with the important quantities involved in the phenomenon. The extreme shallowness allows for a perturbative treatment in the slope $`\alpha `$ of the cone.

FIG. 2. Plot of the number of ripples $`n^{}`$ as a function of the bubble radius $`r_0`$, comparing the experimental measures (points) with the theoretical predictions (solid lines). These data were gathered using silicone oil of viscosity $`\eta =600\text{Pa}\text{s}`$ and bubbles of thickness $`t30\mu \text{m}`$. The errors in the measurement of $`r_0`$ arise from meniscus effects which are more important in smaller bubbles. The bursting time elapsed up to rippling is measured to be of order of 1 to 5 times $`\left(t/g\right)^{1/2}`$, consistent with our proposed mechanism for the formation of the corrugation. For each experimental realization, the ratio $`r_0/R`$ was measured at the onset of the instability, and the corresponding dependence of $`R`$ on $`r_0`$ was used to get a theoretical curve $`n^{}=n^{}\left(r_0\right)`$. The thinner line displays the prediction for an elastic sheet attached to the plane on which it rests. The liquid, however, can in no way be clamped, and one has to relax the boundary conditions at the base. This leads to a vanishing of the unprimed modes (cf. Eq. (3)), which are unfavorable in terms of both graviational and viscous forces; the fastest growing primed modes lead to the behavior represented by the thicker line. The latter is plotted here for a slope $`\alpha 3^{}(0.05\text{rad})`$ of the cone, consistent with our perturbative treatment, and in agreement with direct observation. The dashed line represents the best fit of the scaling form $`n^{}\left(\mu gR^3/K\right)^{1/2}`$ , where $`R`$ is chosen as the relevant length scale. If $`R`$ is replaced by $`r_0`$, the above expression for $`n^{}`$ may be closely fitted (up to an overall multiplicative factor) to our predicted curve, showing that the size of the bubble is the dominant length scale within the present experimental range and conditions. This is consistent with the relaxed boundary conditions, which allow the ripples to be significant close to the outer edge of the bubble (cf. also Fig. 1.(B)). In this way, the ripples trade a bulk gain in gravitational and bending stresses against a cost in stretching in a thin rim close to the outer edge. The increased thickness of the liquid film close to the base further emphasizes this effect, as it reduces the difference in magnitude between a typical stretching and a typical bending stress. The inset shows a top view of the fully developed ripples, from which $`n^{}`$ is measured.