© ESO 2020

1. Introduction

By the time the connection between long gamma-ray bursts (GRBs) and core-collapse supernovae had finally been established (with the association of SN1998bw tothe association of SN1998bw to GRB 980425 – Galama et al. 1998; Patat et al. 2001 – and that of SN2003dh to GRB 030329 – Hjorth et al. 2003; Stanek et al. 2003; Matheson et al. 2003), the community started to realise that long GRB (LGRB) jets must confront themselves with the stellar envelope before being able to expand in the interstellar space up to transparency. This early phase of interaction with the ambient material can leave imprints on the jet properties, whose modelling can be used to constrain the properties of the progenitor (Matzner 2003). Numerical simulations (e.g. Aloy et al. 2000; Zhang et al. 2003; Morsony et al. 2007; Xie & MacFadyen 2019) and (semi-)analytical models (e.g. Ramirez-Ruiz et al. 2002; Matzner 2003; Morsony et al. 2007; Bromberg et al. 2011) have been used to investigate the dynamics of this interaction. In the following years, evidence accumulated in favour of the idea that short GRBs (SGRBs) had a different progenitor (see e.g. Fox et al. 2005, whose deep upper limits could exclude the presence of a supernova associated to GRB 050709); the most favoured option deals with a compact object merger involving at least one neutron star (Narayan et al. 1992). The first investigation of the interaction of a SGRB jet with a post-neutron star merger environment was presented by Aloy et al. (2005), whose numerical simulations found the effect of the interaction with the environment to be significant on the jet propagation and on the determination of its structure. Subsequent, increasingly refined simulations (e.g. Nagakura et al. 2014; Just et al. 2016; Gottlieb et al. 2018b; Xie et al. 2018; Duffell et al. 2018; Geng et al. 2019) confirmed these results, showing that in some cases the merger ejecta could even be dense enough so as to choke the jet.

All of these simulations, regardless of the progenitor type, also show that the distribution of energy per unit solid angle and the average Lorentz factor of the jet after breakout are strong functions of the angle from the jet axis, that is, realistic jets are structured (as opposed to the uniform, or “top-hat”, jet widely assumed in GRB modelling) and the structure is determined by both the properties of the jet at launch and its interaction with the progenitor ambient medium. Lipunov et al. (2001) and Rossi et al. (2002) were the first to realise that the presence of jet structure could lead to a radically different appearance of GRB jets when seen under different viewing angles and that this could be one of the main effects in shaping the GRB luminosity function. The LGRB luminosity function was later shown to be consistent with a single quasi-universal jet structure with a uniform core and a steep decrease of energy density away from the jet axis (Pescalli et al. 2015; Salafia et al. 2015). Kumar & Granot (2003), Granot & Kumar (2003), and Rossi et al. (2004) were the first to provide a detailed modelling of the afterglow emission from a structured jet, showing that the most prominent differences, with respect to the uniform jet model, are for off-axis observers. Structured jets have been invoked to explain some particular LGRBs (such as GRB 030329, Berger et al. 2003) and more recently, in the context of SGRBs, electromagnetic counterparts of gravitational wave events (e.g. Lazzati et al. 2017a; Lamb & Kobayashi 2018; Kathirgamaraju et al. 2018; Salafia et al. 2018).

Despite these advancements, the uniform jet model prevailed until very recently for its simplicity and for its success in describing at least a fraction of GRB afterglows. The first compelling evidence of a structured jet, indeed, came only from observations of GRB 170817A (e.g. Ghirlanda et al. 2019; D’Avanzo et al. 2018; Lazzati et al. 2018; Lyman et al. 2018; Mooley et al. 2018a; Margutti et al. 2018; Resmi et al. 2018; Troja et al. 2018), the GRB associated to GW170817, the first neutron star merger detected in gravitational waves by LIGO and Virgo (Abbott et al. 2017a,b). This has been the first GRB jet to be conclusively observed off-axis, despite an associated prompt emission was detected. Its afterglow multi-wavelength light curves, the apparently superluminal motion (Mooley et al. 2018a) and small projected size (Ghirlanda et al. 2019) of its radio image, and the late time steep decay (Lamb et al. 2019; Mooley et al. 2018b) could only be explained invoking a structured jet seen under a 15–25° viewing angle. Intriguingly, both the prompt and afterglow emission of a jet with the same structure, located at the median distance and within the median interstellar medium density of the previously known SGRB population, if seen on-axis, would have fallen right in the middle of the known SGRB population (Salafia et al. 2019), hinting at a quasi-universal SGRB jet structure.

In this work we present a semi-analytical model of the interaction of the jet with the ambient medium, which represents an attempt at identifying and modelling the main aspects of this complex physical process. We used this model to construct synthetic populations of LGRB and SGRB jets, under the assumption that the properties of these jets at launch and those of their progenitors vary little within the population. This leads naturally to a quasi-universal structure. The comparison of the luminosity distributions of these synthetic populations with the observation-based reconstructed luminosity functions of GRBs shows a general agreement. This provides support to the idea that GRB jets share a quasi-universal structure, even though several uncertainties and caveats have yet to be addressed.

2. Model of the jet propagation through the ambient medium

2.1. Initial remarks

In this and the following section, we describe our model for the jet propagation through, and breakout from, the ambient medium. The model is purely hydrodynamic, that is, we do not include a description of the effects of the magnetic field. Since the main candidate mechanism for jet launching (Blandford & Znajek 1977) requires a strong magnetic field at the base, this may seem a major oversight. The propagation of a Poynting-flux-dominated jet into an ambient medium can be in principle much different from that of a hydrodynamic jet (e.g. Komissarov 1999; Lyubarsky 2009; Bromberg et al. 2014). On the other hand, as initially argued by Levinson & Eichler (2005) based on analytical considerations, and later demonstrated by Bromberg & Tchekhovskoy (2016) by means of relativistic MHD simulations, when such a Poynting-flux-dominated jet propagates in a dense ambient medium and starts being collimated by the cocoon, it becomes unstable to internal kink modes, which lead to magnetic reconnection at the collimation nozzle. This dissipates the magnetic field down to equipartition¹, leading to a propagation that resembles that of a purely hydrodynamic, hot jet (Levinson & Eichler 2005; Bromberg & Tchekhovskoy 2016). External kink modes can still cause the jet head to oscillate in the transverse direction, effectively increasing the jet working surface and therefore slowing down the propagation (Bromberg & Tchekhovskoy 2016), but this effect is progressively reduced during the propagation as long as the ambient material density profile decreases faster than²z⁻² (which should always be the case in GRB progenitors). We thus expect only a minor impact of this simplification on our conclusions, and we defer a more detailed consideration of the inclusion of the magnetic field to a future work.

Let us therefore consider a relativistic jet launched by some central engine located at the origin of our coordinate system. The jet is represented by an outflow that initially moves radially within a cone of half-opening angle θ_j, 0, directed towards the z axis (we employ cylindrical coordinates), with a constant luminosity L_j. The base of the jet is located at a height z_base. If the central engine is surrounded by some ambient material, the jet will collide with it and form a forward shock (that propagates into the ambient medium) and a reverse shock (where the jet material enters the shocked region). We call “head” the region comprised between the forward and reverse shock. Figure 1 shows a sketch of the described setting, with key quantities annotated, to be used as reference throughout this section.

Fig. 1. Sketch of the various components of the jet propagation model, with some key quantities annotated. All quantities are defined and described in the text. Typical radii (order of magnitude in cm) are shown in red on the right, corresponding to the progenitors and average jet properties we employ when computing our synthetic jet populations (note that the collimation-related radii zc, zcoll and $\widehat{z}$ change during the jet propagation).

2.2. Jet head advancement

As shown for example by Marti et al. (1995) and Matzner (2003), the jet head velocity can be estimated by balancing the ram pressure of the jet and that of the ambient material in the head rest frame. In the case of a moving ambient medium, Murguia-Berthier et al. (2017) showed that the resulting head velocity is given by

$$\begin{array}{c}\hfill {\beta }_{\mathrm{h}}=\frac{{\beta }_{\mathrm{a}}+{\beta }_{\mathrm{j}}\sqrt{\stackrel{\sim }{L}}}{1+\sqrt{\stackrel{\sim }{L}}}\end{array}$$(1)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{L}=\frac{{\rho }_{\mathrm{j}}\prime h_{\mathrm{j}}\prime {\mathrm{\Gamma }}_{\mathrm{j}}^2}{{\rho }_{\mathrm{a}}\prime {\mathrm{\Gamma }}_{\mathrm{a}}^2}=\frac{L_{\mathrm{j}}}{\pi {\theta }_{\mathrm{j}}^2{z}^2{\rho }_{\mathrm{a}}{\mathrm{\Gamma }}_{\mathrm{a}}{c}^3}·\end{array}$$(2)

Here the subscripts a and j indicate respectively the ambient and jet material, primed quantities are measured in the comoving frame of the respective fluid, β indicates the velocity in units of the speed of light c, Γ = (1 − β²)^−1/2 is the Lorentz factor, ρ is the rest mass density, h is the specific dimensionless enthalpy, and θ_j is the jet head half-opening angle.

Numerical relativistic hydrodynamics simulations generally show that the jet head proceeds slower than predicted by Eq. (1) in presence of jet collimation by the ambient medium. As detailed in Harrison et al. (2018), the agreement with simulations in the case of a static ambient medium can be restored by replacing $\stackrel{\sim }{L}$ in Eq. (1) with the effective value ${\stackrel{\sim }{L}}_{\mathrm{eff}}=A^2(\stackrel{\sim }{L},{\theta }_{\mathrm{j},0},\mathrm{\Omega })\stackrel{\sim }{L}$, where

$$\begin{array}{c}\hfill A(\stackrel{\sim }{L},{\theta }_{\mathrm{j},0},\mathrm{\Omega })=\{\begin{array}{cc}0.35\hfill & \stackrel{\sim }{L}\le 1\hfill \\ 0.35{\stackrel{\sim }{L}}^{\frac{0.46}{{log}_{10}{\stackrel{\sim }{L}}_{\mathrm{col}}}}\hfill & 1<\stackrel{\sim }{L}\le {\stackrel{\sim }{L}}_{\mathrm{col}}\hfill \\ 1\hfill & \stackrel{\sim }{L}>{\stackrel{\sim }{L}}_{\mathrm{col}}\hfill \end{array}\end{array}$$(3)

where $\stackrel{\sim }{L}$ is computed using Eq. (2) evaluated at the jet head, and ${\stackrel{\sim }{L}}_{\mathrm{col}}={\theta }_{\mathrm{j},0}^{-4/3}{(16\mathrm{\Omega }/3)}^{2/3}$. Here Ω is a functional of the ambient medium density profile whose exact form is described in Harrison et al. (2018). In what follows, we employ this correction, assuming it to be valid also in the case of a moving ambient medium.

2.2.1. The cocoon

As the jet head advances, the swept ambient material is cast aside, forming an over-pressured region surrounding the jet, which is usually called “cocoon”. We assumed the energy flow from the head to the cocoon Ė_c to be a fraction η of the jet energy flow through the reverse shock, that is ηL_j(β_j − β_h). We estimated η(t) as the fraction of the head volume that is in causal contact with the cocoon at a given time t. Let us assume the sound speed in the head to be $c_{\mathrm{s}}=c/\sqrt{3}$. The time-scale associated to the head advancement, as measured in the head comoving frame, is t_h ∼ z_h/β_hΓ_hc, so the region of the head in causal contact with the cocoon extends to a distance $l\sim z_{\mathrm{h}}/\sqrt{3}{\beta }_{\text{h}}{\mathrm{\Gamma }}_{\text{h}}$ inside the head. Assuming the head to be a cylinder of radius r_j = θ_jz_h, the fraction of its volume being in causal contact with the cocoon at a given time is then given by

$$\begin{array}{c}\hfill \eta =\{\begin{array}{cc}\frac{2}{\mu }-\frac{1}{{\mu }^2}\hfill & \hfill \mu >1\\ 1\hfill & \hfill \mu \le 1\end{array}\end{array}$$(4)

where $\mu =\sqrt{3}{\theta }_{\text{j}}{\mathrm{\Gamma }}_{\text{h}}{\beta }_{\text{h}}$. This is similar to the prescription given in Bromberg et al. (2011), with the advantage that it is a continuous function of the jet head velocity.

We assumed the cocoon to be at rest in a frame that moves upwards at a speed ${\overline{\beta }}_{\mathrm{a}}$ which is the rest-mass averaged speed of the ambient medium from z₀ (the height of its base) to z_h, namely

$$\begin{array}{c}\hfill \overline{{\beta }_{\mathrm{a}}}=\frac{{\int }_{z_0}^{z_{\mathrm{h}}}{\rho }_{\mathrm{a}}(z){\beta }_{\mathrm{a}}(z)\mathrm{d}z}{{\int }_{z_0}^{z_{\mathrm{h}}}{\rho }_{\mathrm{a}}(z)\mathrm{d}z}\end{array}$$(5)

and we define ${\overline{\mathrm{\Gamma }}}_{\mathrm{a}}={(1-{\overline{{\beta }_{\text{a}}}}^2)}^{-1/2}$. The energy of the cocoon is given by

$$\begin{array}{c}\hfill E_{\mathrm{c}}(t)={\displaystyle {\int }_0^t\eta (t)L_{\mathrm{j}}({\beta }_{\mathrm{j}}-{\beta }_{\mathrm{h}}(t))\mathrm{d}t}\end{array}$$(6)

where t is measured in the central engine rest frame and the jet luminosity is evaluated at the retarded time t − z_h/β_jc. The term in parentheses accounts for the relative velocity between the jet and the head.

2.2.2. Cocoon pressure

We assumed the cocoon to be cylindrical, extending from the base of the ambient material located at a height z₀ (which could be above the jet base if a cavity is present around the central engine, and moves with the ambient material if the wind stops being injected from the central region) and extending up to the jet head at z_h (where the head bow shock is present). Let us indicate the cocoon cylindrical radius as r_c. We assumed the cocoon to be radiation-dominated and neglected pressure gradients, so that its pressure is $P_{\mathrm{c}}=E_{\mathrm{c}}/3\pi r_{\mathrm{c}}^2(z_{\mathrm{h}}-z_0){\overline{\mathrm{\Gamma }}}_{\mathrm{a}}$. The cocoon expands sideways at a speed β_c which is obtained by balancing the cocoon pressure and the ram pressure from the external medium. To keep with our assumption of cylindrical cocoon shape, we averaged out this speed vertically, namely ${\beta }_{\mathrm{c}}={(1+{\overline{\rho }}_{\mathrm{a}}{c}^2/P_{\mathrm{c}})}^{-1/2}$, where ${\overline{\rho }}_{\mathrm{a}}={(z_{\mathrm{h}}-z_0)}^{-1}{\int }_{z_0}^{z_{\mathrm{h}}}{\rho }_{\text{a}}(z,t)\mathrm{d}z$ is the altitude-averaged ambient density.

2.2.3. Reconfinement shock

The jet is surrounded by the cocoon, which exerts an approximately uniform lateral pressure P_c inwards on the jet. If this pressure is higher than thermal pressure in the jet $P_{\text{j}}~L{z}_{\text{base}}^2/4\pi {\theta }_{\text{j},\text{0}}^2{z}^{4}c$ (Matzner 2003, assuming an adiabatic jet), a reconfinement shock forms (Komissarov & Falle 1997; Morsony et al. 2007). The shape of the shock is given by the balance between the component of the jet ram pressure normal to the shock and the cocoon pressure. In the frame of the ambient material, within small angles from the jet axis, the balance equation is given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{j}}\prime h_{\mathrm{j}}\prime {\mathrm{\Gamma }}_{\mathrm{j}}^2{\mathrm{\Gamma }}_{\mathrm{a}}^2{({\beta }_{\mathrm{j}}-{\beta }_{\mathrm{a}})}^2{\mathrm{\Gamma }}_{\mathrm{a}}^{-2}{(\frac{r_{\mathrm{s}}}{z}-\frac{\mathrm{d}{r}_{\mathrm{s}}}{\mathrm{d}z})}^2+P_{\mathrm{j}}=P_{\mathrm{c}}\end{array}$$(7)

where r_s is the cylindrical radius of the reconfinement shock, and z is measured in the central engine frame (the ${\mathrm{\Gamma }}_{\text{a}}^{-2}$ factor above accounts for its transformation to the ambient medium frame). Assuming P_c ≫ P_j we neglected the latter and wrote the solution of the above ordinary differential equation as

$$\begin{array}{c}\hfill r_{\mathrm{s}}(z)={\theta }_{\mathrm{j},0}(1+A{z}_{\mathrm{c}})z-{\theta }_{\mathrm{j},0}A{z}^2\end{array}$$(8)

where $A=\sqrt{\pi c{\beta }_{\text{j}}{P}_{\mathrm{c}}/L_{\mathrm{j}}{({\beta }_{\text{j}}-{\beta }_{\text{a}})}^2}$, and z_c is the height at which the jet thermal pressure equals the cocoon pressure (Bromberg et al. 2011). The reconfinement shock thus converges to the jet axis at a height $\widehat{z}=A^{-1}+z_{\mathrm{c}}$. Following Bromberg et al. (2011), owing to the parabolic shape of the reconfinement shock, we assumed the jet cross section radius r_j to be constant above $(\widehat{z}+z_{\mathrm{c}})/2$, being given by $r_{\mathrm{j}}={\theta }_{\mathrm{j},0}(\widehat{z}+z_{\mathrm{c}})/2$. For simplicity, any variation in the conditions (especially the cocoon pressure) near z_c is assumed to affect instantaneously the reconfinement shock and thus the jet head opening angle.

3. Jet breakout and development of jet structure

3.1. Breakout condition

Most realistic progenitor density profiles feature a smooth decrease of the mass density beyond some radius, therefore we need an unambiguous way of defining the breakout radius. Since the jet head forward shock is radiation-dominated, it remains confined to the progenitor interior as long as the Thomson optical depth $\tau (z_{\text{h}})={\int }_{z_{\text{h}}}^{\infty }({\sigma }_{\text{T}}/m_{\text{p}}){\rho }_{\text{a}}(z)\text{d}z$ from the head height z_h outwards is larger than unity³ (here σ_T and m_p indicate the Thomson cross section and the proton mass respectively). We therefore took τ = 1 as our breakout condition.

3.2. Post-breakout evolution

As soon as the jet head breaks out of the ambient medium (i.e. the progenitor star in the long GRB scenario, or the outer edge of the compact binary merger ejecta in the short GRB scenario), at a lab-frame time t_bo and a height z_bo, both the jet and the cocoon are free to expand in the interstellar medium (ISM). As the cocoon material starts to flow out of the new open channel, the shock at its sides stalls and its pressure drops, making the jet collimation at its base soon ineffective (Lazzati & Begelman 2005). Since the jet in our model is effectively collimated at a height $z_{\mathrm{coll}}=(\widehat{z}+z_c)/2$, the information on the cocoon pressure loss starts affecting the jet opening angle only after a delayed time t_delay = (z_bo − z_coll)/c_s, where $c_{\mathrm{s}}=c/\sqrt{3}$ is the sound speed⁴. After this time, we assumed the cocoon pressure at z_coll to drop exponentially as

$$\begin{array}{c}\hfill P_{\mathrm{c}}(t)=P_{\mathrm{c},\mathrm{bo}}exp(-\frac{c_{\mathrm{s}}(t-t_{\mathrm{bo}}-t_{\mathrm{delay}})}{z_{\mathrm{bo}}})\end{array}$$(9)

due to the cocoon material flowing at sound speed out of the open channel. As detailed in Lazzati & Begelman (2005), the pressure drop causes the jet half-opening angle to increase exponentially

$$\begin{array}{c}\hfill {\theta }_{\mathrm{j}}(t)={\theta }_{\mathrm{j},\mathrm{bo}}exp\left(\alpha \frac{c_{\mathrm{s}}(t-t_{\mathrm{bo}}-t_{\mathrm{delay}})}{z_{\mathrm{bo}}}\right)\end{array}$$(10)

until it eventually reaches the base jet half-opening angle θ_j, 0, or until the jet injection stops. Here θ_j, bo is the half-opening angle of the head at breakout, and α is a parameter that depends on the details of the transient acceleration phase of the jet after breakout (Lazzati & Begelman 2005). We set this value to 1/8 (it was 1/4 in Lazzati & Begelman 2005) based on a comparison with numerical simulations (Appendix A.2). Material that breaks out within an angle ${\theta }_{\mathrm{j}}(t)<{\mathrm{\Gamma }}_{\mathrm{j}}^{-1}$ subsequently expands laterally to ${\mathrm{\Gamma }}_{\mathrm{j}}^{-1}$ (Lazzati & Begelman 2005), so we actually took θ_j(t) as the maximum between the value computed from Eq. (10) and ${\mathrm{\Gamma }}_{\mathrm{j}}^{-1}$.

This simple modelling neglects the entrainment (and the possible development of instabilities) between the jet and the surrounding cocoon material, which can transfer some of the jet energy outwards. As a simple effective description of this transfer, we assumed that the jet energy that flows out of the progenitor (after breakout) during each time interval dt is spread over the latitudinal angle in a Gaussian fashion, with a sigma equal to θ_j(t). The final jet structure is then obtained by integrating over t, namely

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{jet}}}{\mathrm{d}\mathrm{\Omega }}(\theta )={\displaystyle {\int }_{t_{\mathrm{bo}}}^{T_{\mathrm{jet}}+z_{\mathrm{bo}}/{\beta }_{\mathrm{j}}c}{L}_{\mathrm{j}}\mathcal{G}({\theta }_{\mathrm{j}}(t),\theta )\mathrm{d}t}\end{array}$$(11)

where

$$\begin{array}{c}\hfill \mathcal{G}({\theta }_{\mathrm{j}},\theta )\propto exp[-\frac{1}{2}{\left(\frac{\theta }{{\theta }_{\mathrm{j}}}\right)}^2]{cos}^2\theta \end{array}$$(12)

is the assumed angular energy spread function, and $2\pi {\int }_0^1\mathcal{G}({\theta }_{\text{j}},\theta )\text{d}\cos \theta =1$. The cos²θ term is inserted to ensure that the function goes to zero as θ → π/2. The quantity T_jet in Eq. (11) represents the jet duration in the central engine frame, so that T_jet + z_bo/βjc is the time when the jet ceases to flow out of the ambient medium. With this description, a jet that propagates uncollimated by the ambient material develops a Gaussian structure $\text{d}E/\text{d}\mathrm{\Omega }\underset{\sim }{\propto }\exp [-{(\theta /{\theta }_{\text{j},\text{0}})}^2/2]$. Similarly, a jet that is collimated, but whose duration is short compared to the timescale t_θ = t_delay + 8z_bo/c_s over which the jet cross section increases, has $\text{d}E/\text{d}\mathrm{\Omega }\underset{\sim }{\propto }\exp [-{(\theta /{\theta }_{\text{j},\text{bo}})}^2/2]$. Cases in between will feature a shallower decrease followed by a Gaussian cut-off. For both short and long GRBs, typical values of t_θ are longer than the jet duration, thus the typical GRB jet structure is most likely narrow with a steep fall off outside the jet core.

4. Cocoon structure

After breakout, the shocked ambient material in the cocoon is free to expand under the effect of its own internal pressure. As discussed in the preceding section, its material initially blows out of the open channel formed after the jet breakout, and then expands (anisotropically) in the surrounding space. Based on the standard theory of radiation-dominated relativistic fireballs (Cavallo & Rees 1978; Piran et al. 1993), we can estimate its average terminal Lorentz factor as ${\overline{\mathrm{\Gamma }}}_{\mathrm{c}}\sim 1+E_{\mathrm{c}}/M{c}^2$, where

$$\begin{array}{c}\hfill M={\displaystyle {\int }_0^{t_{\mathrm{bo}}}\pi r_{\mathrm{j}}^2(t){\rho }_{\mathrm{a}}(z_{\mathrm{h}}(t))({\beta }_{\mathrm{h}}(t)-{\beta }_{\mathrm{a}}(t))c\mathrm{d}t}\end{array}$$(13)

is the mass in the inner cocoon (i.e. the ambient mass swept and cast aside as the jet head propagates). Initially, it expands within an angle ∼r_c, bo/z_bo, where r_c, bo is the cocoon cylindrical radius at breakout, but its internal pressure causes it to expand laterally to an angle ${\theta }_{\mathrm{c}}\sim max(r_{\mathrm{c},\mathrm{bo}}/z_{\mathrm{bo}},\arcsin ({\overline{\mathrm{\Gamma }}}_{\mathrm{c}}^{-1}))$. These simple arguments give us the Lorentz factor and angular scale of the cocoon structure.

4.1. Energy angular structure

Numerical relativistic hydrodynamical simulations (Lazzati et al. 2017b; Lazzati & Perna 2019) suggest that the cocoon typically features an approximately exponential distribution of energy per unit solid angle. Inspired by that, we assume the following ansatz energy distribution

$$\begin{array}{c}\hfill \frac{\mathrm{d}{E}_{\mathrm{c}}}{\mathrm{d}\mathrm{\Omega }}(\theta )=K(E_{\mathrm{c}})exp(-\frac{\theta }{{\theta }_{\mathrm{c}}}){cos}^2\theta \end{array}$$(14)

where θ_c is defined in the preceding section, and K(E_c) is chosen so that the total kinetic energy in the cocoon is E_c, namely

$$\begin{array}{c}\hfill K(E_{\mathrm{c}})=\frac{E_{\mathrm{c}}}{2\pi {\int }_0^{1}exp(-\frac{\theta }{{\theta }_{\mathrm{c}}}){cos}^2\theta d(cos\theta )}\end{array}$$(15)

which implies the assumption that all internal energy is converted to kinetic energy.

4.2. Lorentz factor angular structure

We still need to define how the average Lorentz factor varies with the angular distance from the jet axis. We assume the jet material, where not mixed with the cocoon, to reach a terminal Lorentz factor Γ_j (which is a free parameter of the model). Based on the natural expectation that the cocoon material closer to the jet axis moves faster than that at larger angles, we assign to the cocoon material the ansatz Lorentz factor profile (similar to Lazzati & Perna 2019)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{c}}(\theta )-1=({\overline{\mathrm{\Gamma }}}_{\mathrm{c}}-1)exp(-\omega \frac{\theta }{{\theta }_{\mathrm{c}}})\end{array}$$(16)

where ${\overline{\mathrm{\Gamma }}}_{\mathrm{c}}$ and θ_c have been defined in the preceding section, and ω is a free parameter which we set to ω = 1/3 based on comparison with numerical simulations (Appendix A.3)

5. Final structure

After breakout, the jet and cocoon form a single structure (often referred to as a structured jet), which soon reaches homologous (i.e. ballistic) expansion. While some rearrangement of the internal structure is still possible in this phase, we assume it to be negligible and simply compute the final angular energy distribution as the sum of the jet and cocoon energies, namely

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}\mathrm{\Omega }}(\theta )=\frac{\mathrm{d}{E}_{\mathrm{jet}}}{\mathrm{d}\mathrm{\Omega }}(\theta )+\frac{\mathrm{d}{E}_{\mathrm{c}}}{\mathrm{d}\mathrm{\Omega }}(\theta ).\end{array}$$(17)

For what concerns the Lorentz factor of the structured jet, we compute it as the mass-weighted average of those of the jet and cocoon material, namely

$$\begin{array}{c}\hfill \mathrm{\Gamma }(\theta )=\frac{\frac{\mathrm{d}{E}_{\mathrm{jet}}}{\mathrm{d}\mathrm{\Omega }}(\theta )+\frac{\mathrm{d}{E}_{\mathrm{c}}}{\mathrm{d}\mathrm{\Omega }}(\theta )}{{\mathrm{\Gamma }}_{\mathrm{j}}^{-1}\frac{\mathrm{d}{E}_{\mathrm{jet}}}{\mathrm{d}\mathrm{\Omega }}(\theta )+{\mathrm{\Gamma }}_{\mathrm{c}}^{-1}(\theta )\frac{\mathrm{d}{E}_{\mathrm{c}}}{\mathrm{d}\mathrm{\Omega }}(\theta )}\end{array}$$(18)

6. The luminosity distribution of GRBs from structured jets

Equipped with the model presented in the preceding sections, we proceeded to constructing a synthetic population of GRB jets, each with its own structure set by the interaction with the ambient medium. We treated long GRBs (LGRBs) and short GRBs (SGRBs) as separate, independent populations. For each, we chose a single, representative progenitor model, into which we injected jets whose properties are distributed within a narrow range. This results in a population of structured jets (those which successfully punch through the ambient medium) and choked jets (those which do not). Based on a simple modelling of the prompt emission, we computed the luminosity distribution of each population, assuming isotropic viewing angles. Finally, we compared the results with the luminosity functions derived from the observations.

6.1. Prompt emission model

In order to construct our luminosity distributions, we needed to model the jet prompt emission isotropic-equivalent peak luminosity L_iso that results from a given jet structure, as a function of the viewing angle. Given the jet injection duration T_jet, we computed the GRB duration as T_GRB = T_jet − (t_bo − z_bo/β_jc) and we assumed it to be independent of the viewing angle. This assumption is based on the idea that the gamma-ray emission is powered by some form of energy dissipation within the jet (e.g. internal shocks or magnetic reconnection), and that it consists of a series of short episodes taking place at some typical radius. Under these assumptions, the duration of the GRB is set by the post-breakout jet injection duration, and does not depend strongly on the viewing angle (see e.g. Salafia et al. 2016). As a caveat, we note that this requires the angular timescale R/Γ²c at the photospheric radius R (at the angles that contribute significantly to the emission for a given viewing angle) to be shorter than the duration of the burst itself, otherwise the angular time spread would dominate.

We computed the prompt emission isotropic equivalent energy E_iso(θ_v) following Salafia et al. (2015, 2019), that is, we assumed 10 percent of the kinetic energy at each angle to be converted into gamma-ray radiation, and we integrated the emission over the jet accounting for relativistic beaming. We assumed the prompt emission light curve to have a triangular shape, so that L_iso(θ_v) = 2E_iso(θ_v)/T_GRB. This is a crude approximation, given the diversity of GRB light curves (especially in LGRBs), but we adopted it for its simplicity.

For the present study, we only applied this model to cases where the jet successfully breaks out of the progenitor. We defer the implementation of more detailed emission models and the inclusion of the possible prompt emission in the case of a choked jet to future work.

6.2. Long GRBs

6.2.1. Ambient medium density profile

For our representative LGRB progenitor ambient medium, we took the density profile from the stellar model 16TI of Woosley & Heger (2006). This model has been used in several previous studies of LGRB jets (e.g. Morsony et al. 2007; López-Cámara et al. 2013). It represents a star with an initial mass of M = 16 M_⊙, low metallicity Z = 10⁻² Z_⊙, and a large initial angular momentum J = 3 × 10⁵² erg s, which is evolved to pre-core collapse using the KEPLER code (Weaver et al. 1978). The original simulation data⁵ features a low-mass (ΔM ∼ 0.1 M_⊙) extended (ΔR ∼ 7 × 10¹¹ cm) envelope – represented by a single point in the adaptive grid – which is considered (Woosley & Heger 2006) an artefact due to the treatment of rotation in a one-dimensional simulation. We therefore remove this point from the density profile; we add, for consistency, a stellar wind with typical Wolf-Rayet mass loss rate Ṁ_w = 10⁻⁵ M_⊙ yr⁻¹ and velocity v_w = 1000 km s⁻¹, even though this has a negligible impact on the modelled jet dynamics. The final density profile is shown in Fig. 2. We note that our results, as detailed below, do not depend strongly on the detailed shape of the progenitor density profile: we obtain similar results employing the progenitor model of Duffell & MacFadyen (2015). A more systematic investigation of the dependence on the progenitor properties will be presented in a future work.

Fig. 2. Radial density profile from our representative LGRB progenitor (model 16TI of Woosley & Heger 2006), which represents a low-metallicity, high-angular momentum, massive star prior to collapse. The low-density power-law tail, visible in the lower-right corner of the figure, has been added by ourselves to the original density profile of the model. It represents a stellar wind with typical Wolf-Rayet mass loss rate Ṁw = 10−5 M⊙ yr−1 and velocity vw = 1000 km s−1.

6.2.2. Jet properties at launch

Given the large uncertainties on the jet launching mechanism and on its initial collimation, we simply assumed all jets to have a fixed half opening angle θ_j, 0 = 0.25 rad ≈ 14° at injection. We also assumed their terminal Lorentz factor to be Γ_j = 100 in all cases (while higher values are actually observed in some cases, this seems to be a typical value at least for LGRBs, see Ghirlanda et al. 2018). We extracted the jet luminosity L_j from a log-normal distribution centred at μ = 2 × 10⁴⁹ erg s⁻¹ ≡ ⟨L⟩ with a dispersion σ = 0.65 dex. Similarly, we extracted the duration from a log-normal distribution with μ = 30 s ≡ ⟨T⟩ and σ = 0.45 dex. These values are educated guesses, based on the fact that the typical collimation-corrected kinetic energy of LGRBs is around 10⁵¹ erg (Goldstein et al. 2016), and their average rest-frame duration is just below 30 s (Salafia et al. 2015). We chose z_base = 5 × 10⁶ cm as our injection height. We simulated a total of 1000 jets, 94% of which successfully broke out of the progenitor star.

6.2.3. Resulting LGRB structures

Figure 3 shows the distribution of jet structures for the successful jets in our synthetic LGRB population. The red solid line represents the median value of the kinetic energy per unit solid angle (upper panel) and Γ − 1 (lower panel, respectively) at a given angle θ. The orange and yellow shaded regions contain 50% and 90% of the values at each fixed θ, respectively. The figure shows that our LGRB jets are on average very narrow, featuring a core of ≲1.5 deg, with a steep fall off of the kinetic energy density outside. The cocoon is in general quite energetic, but wide, so that its typical energy density per unit solid angle is lower by around four orders of magnitude with respect to the jet core. The very narrow jet opening angle might seem at odds with those derived from jet breaks in LGRB afterglow light curves (Berger 2014), but we caution that (i) the jet structure (especially if narrow, see e.g. Granot & Piran 2012) is expected to evolve after breakout due to jet lateral spread, and it may thus be different at the afterglow stage (see also Gill et al. 2019a), and (ii) in the context of structured jet light curves, jet-break-like features are related to the viewing angle rather than to the jet opening angle (Kumar & Granot 2003; Rossi et al. 2004). Interestingly, the general features we found for the LGRB jet structure – a narrow opening angle and a faster angular decrease in energy than in Lorentz factor – are in agreement with the observational constraints derived by Beniamini & Nakar (2019).

Fig. 3. Distribution of jet structures in our LGRB population. In each panel, at each fixed angle θ, the yellow (resp. orange) filled area encloses 90% (resp. 50%) of the structures in the population, while the red line shows the median of the distribution at that angle. Top panel: energy structures (kinetic energy per unit solid angle as a function of the angle θ from the jet axis). Bottom panel: velocity structures (average Lorentz factor diminished by one, as a function of θ).

6.2.4. Luminosity distribution

Figure 4 shows the luminosity distribution of our LGRB population, computed as described in the preceding sections. The thick solid red histogram shows the local (z = 0) rate density per logarithmic luminosity bin, dR₀/d log₁₀(L_iso/(erg s⁻¹)), constructed from our model population, assuming a total local rate R₀ = 600 yr⁻¹ Gpc⁻³. The thin, coloured histograms show the luminosity distributions of jets seen within 2 deg from the jet axis (yellow), between 2 and 10 deg (dark blue) and at viewing angles larger than 10 deg (light blue). The black points show the actual luminosity distribution of LGRBs from Pescalli et al. (2015), who collected and updated binned rate estimates from previous works: the points to the right of L_iso = 10⁵⁰ erg s⁻¹ come from the reconstructed luminosity function from Wanderman & Piran (2009); the leftmost point is the low-lumnosity GRB rate from Soderberg et al. (2004), updated as described in Pescalli et al. (2015); the three remaining points are lower limits derived in Pescalli et al. (2015). The general features of the model population are in good agreement with the observed points, which is remarkable given the fact that we based our parameter distributions on educated guesses This shows that the luminosity distribution of LGRBs can be successfully interpreted in a scenario where both jets and progenitors have very similar properties, the main parameter behind the diversity being the viewing angle.

Fig. 4. LGRB luminosity distribution of our model population (thick solid red histogram) compared to the observed distribution (black points, from Pescalli et al. 2015). Thinner histograms show the luminosity distributions of three sub-classes corresponding to three viewing angle bins, reported in the legend.

6.2.5. Dependence on the parameters

Figure 5 gives a grasp of how the luminosity function depends on the jet properties at launch. The grey thick line is the fiducial luminosity function (the same as that shown in Fig. 4). Thinner coloured lines show the variation induced by changing one of the parameters of the distributions of jet properties at launch at a time. To produce each line we construct a synthetic population of 1000 jets whose properties are extracted from distributions with the same parameters as those described in Sect. 6.2.2, except one, which is displaced as reported in the legend. The resulting luminosity functions (always constructed only from jets that successfully break out) are normalised to the same total rate as the fiducial one. Here we describe the general trends, and some tentative qualitative explanations. Increasing the average jet luminosity ⟨L⟩ by a factor of 10 (pink line in the left-hand panel of Fig. 5) essentially shifts the whole luminosity function to higher luminosities, keeping the shape similar. On the contrary, lowering ⟨L⟩ by the same factor (purple line in the left-hand panel of Fig. 5), the high-luminosity end is shifted towards lower luminosities, but the low-luminosity end remains essentially unchanged: this is the combined effect of a lower average luminosity after breakout, and of an increased ratio of cocoon to jet energy, which makes the structure shallower. Lowering the average jet injection duration ⟨T⟩ (blue line in the central panel of Fig. 5) or increasing the jet opening angle at launch (yellow line in the right-hand panel) both make the luminosity function slope below L_iso ∼ 10^51.5 erg s⁻¹ steeper: this is due to the larger energy in the cocoon, which makes the structure, broadly speaking, shallower (Pescalli et al. 2015). Also the opposite is true, as shown by the remaining lines.

Fig. 5. Long GRB luminosity function dependence on the parameters of the jet at launch. Each panel shows how the fiducial luminosity function (thick grey solid line) is affected by a change in one of the parameters (coloured thinner lines). See text for an interpretation of the dependencies.

A more detailed account of the dependence of the luminosity function on the parameters and assumptions will be presented in a future work.

6.3. Short GRBs

6.3.1. Ambient medium density profile

Assuming that all (or most) SGRBs are produced by the merger of two neutron stars, and that there is little variation (among different events) in the properties of the various outflows produced during and after the merger, we assumed the homologously expanding ejecta cloud described in Xie et al. (2018) to be representative of the SGRB ambient medium. Current SGRB observations seem to actually indicate some degree of diversity in the associated kilonova emission (Gompertz et al. 2018), which would in turn suggest a variety outflow properties, but it is difficult at present to quantify it. We therefore opted for a single progenitor ambient medium model for simplicity. The density profile along the jet axis (solid red line) and velocity profile (dashed blue line) are shown in Fig. 6.

Fig. 6. Density and velocity profile of the SGRB ambient medium model. The red line shows a snapshot of the density profile (along the polar axis) of the neutron star merger ejecta model described in Xie et al. (2018) (also used in Kasliwal et al. 2017 and Gottlieb et al. 2018b). The blue dashed line shows the corresponding velocity profile (velocity axis on the right). The ejecta are assumed to expand homologously.

6.3.2. Jet properties at injection

As for the LGRB population, we assumed all jets to have a fixed half opening angle θ_j, 0 = 0.25 rad at injection and a terminal Lorentz factor Γ_j = 100. We extracted the jet luminosity L_j from a log-normal distribution centred at μ = 3 × 10⁴⁹ erg s⁻¹ with a dispersion σ = 0.85 dex; we extracted the duration from a log-normal distribution with μ = 0.3 s and σ = 0.45 dex. For SGRBs, since the ambient medium is not static, we needed to define an additional parameter, that is the delay between the neutron star merger (which corresponds to the time when all the ambient material is concentrated at r = 0) and the start of the jet injection. This represents the time it takes for the post-merger system to develop the necessary conditions to launch a jet (see Gill et al. 2019b for an interesting discussion of the delay between GW170817 and GRB170817A and its implications for that system). These conditions are very uncertain, as the actual jet launching mechanism is still debated, but the most likely option seems to be energy extraction from a spinning black hole surrounded by an accreting torus of highly magnetised material, through the mechanism first described in Blandford & Znajek 1977 (see also Tchekhovskoy et al. 2010). For this to happen, (i) the merger remnant must have collapsed to a black hole, (ii) the magnetic field must have been amplified in the torus (mainly by magneto-rotational instabilities) to reach a significant magnetisation, and (iii) an ordered large scale magnetic field structure with a significant poloidal component must have developed. The time for (ii) is likely very short (few dynamical times), as shown by high-resolution GRMHD simulations (e.g. Kiuchi et al. 2018; Kawamura et al. 2016). Process (iii) should take a few Alfvén times (tens of milliseconds, see e.g. Christie et al. 2019, whose simulations resolve the magnetic field amplification, its large scale organisation and the consequent jet launching). In absence of a prompt collapse to a BH (which would anyway leave little matter outside), therefore, the delay is likely dominated by (i). If the merger remnant is a hypermassive proto-neutron star (which is supported by differential rotation), the collapse to a BH takes tens to hundreds of milliseconds; on the other hand, if the merger remnant is a supra-massive neutron star (which is supported by solid-body rotation), the collapse to a BH takes place only after electromagnetic spin down slows down the remnant enough for it to become unable to support itself against self-gravity: in this case, the collapse could take place after several seconds. It is likely, though, that no accretion disc would be left after such a late collapse (Margalit et al. 2015). The most likely scenario for a successful SGRB jet to be launched therefore seems that of a short-lived hypermassive proto-neutron star (see also Shibata et al. 2006). Based on these arguments, we fixed our time delay at Δt_inj = 0.1 s. Of the total 1000 jets we simulated, 88% successfully broke out of the ambient medium.

6.3.3. Distribution of SGRB structures

The jet structure distribution for our synthetic SGRB population features a somewhat larger dispersion with respect to the LGRB population, as shown in Fig. 7. It is interesting to compare this distribution with the jet structure of GRB 170817A derived by Ghirlanda et al. (2019) based on multi-wavelength fitting of the afterglow light curves and of the centroid motion observed in VLBI images. The black lines in Fig. 7 represent the best-fit (solid), one sigma (dashed) and two sigma (dotted) contours of the derived jet structure. These show that a good fraction of jet structures in our synthetic population are compatible⁶ with GRB 170817A.

Fig. 7. Distribution of structures in our SGRB population. The filled areas and the red lines have the same meaning as in Fig. 3. The black lines show the confidence regions (dotted: two sigma, dashed: one sigma, solid: best fit) for the jet structure of GRB 170817A from Ghirlanda et al. (2019).

6.3.4. Luminosity distribution

Figure 8 shows the luminosity distribution of our SGRB population, computed as described in the preceding sections. The thick solid red histogram shows the rate density per logarithmic luminosity bin at redshift zero, dR₀/dlog₁₀(L_iso/(erg s⁻¹)), constructed from our model population, assuming a total local rate R₀ = 300 yr⁻¹ Gpc⁻³. This local rate has been chosen to make the high luminosity end comparable to estimates by Wanderman & Piran (2015) (W15) and Ghirlanda et al. (2016) (G16 – their case a). It falls on the low end of the binary neutron star merger rates derived by The LIGO Scientific Collaboration (2019). The shape of the simulated luminosity distribution is in good agreement with that of W15, which is again remarkable, given that we did not attempt to perform a fitting. In terms of the rate, the consistency with W15 would need a lower total rate (around 80 yr⁻¹ Gpc⁻³). The high-luminosity end (L >  10⁵²erg s⁻¹) is consistent with G16 for this choice of total rate, but the downward turn of the distribution at lower luminosities is not reproduced. However, as shown in the next section, a similar feature appears if the jet opening angle at launch is assumed to be narrower.

Fig. 8. SGRB luminosity distribution of our model population (thick solid red histogram) compared to two luminosity function models based on observations (Wanderman & Piran 2015 – green solid line and shaded area – and Ghirlanda et al. 2016 – purple solid line and shaded area). Thinner histograms show the luminosity distributions of sub-classes of jets belonging to three viewing angle bins, reported in the legend. The grey line and filled horizontal band show the binary neutron star merger local rate based on gravitational wave observations (The LIGO Scientific Collaboration 2019).

Figure 8 reveals another particular feature of our synthetic SGRB population: even the high-end of the luminosity function is dominated by jets seen (slightly) off-axis. Given the uncertainties in the modelling and the degeneration in the parameters (see also next subsection), though, we defer a deeper investigation of this feature to a future work.

6.3.5. Dependence on the parameters

Figure 9 illustrates the dependence of the luminosity function on the jet properties at launch. The meaning of the various curves is the same as in Fig. 5. As in the case of long GRBs, most of the trends can be understood intuitively. Similarly to LGRBs, a change in the average jet luminosity ⟨L⟩ at launch shifts the high-end of the luminosity function correspondingly (see Fig. 9, left-hand panel). On the other hand, the low-end, which is dominated by jets seen far off-axis, does not change appreciably. This is due to the homologous expansion of the ambient material (Duffell et al. 2018): jets with a lower luminosity produce slower-moving heads, but as the ambient medium expands these confront themselves with lower densities with respect to heads that move faster. As a consequence, it turns out that the breakout success (and therefore the energy stored in the cocoon after breakout) does not depend explicitly on the jet luminosity (at a fixed initial jet opening angle – Duffell et al. 2018). An increase (resp. decrease) in the average jet injection duration ⟨T⟩ shifts the whole luminosity function to lower (resp. higher) luminosities (central panel of Fig. 9). This happens because in this case the jet duration is comparable to the jet breakout time, which has two consequences: (i) the cocoon energy is comparable to that in the jet and (ii) the post-breakout duration is short, which boosts the cocoon luminosity (due to the particular way we model it, see Sect. 6.1). Since slightly off-axis jets contribute even at the high-end of the luminosity function (with this particular choice of parameters), this impacts the whole luminosity function. Increasing the jet opening angle at launch (yellow line in the right-hand panel of Fig. 9) increases the average energy in the cocoon, but it also increases its mass and therefore its opening angle. This shifts the low-luminosity peak of the luminosity function to the right, but leaves the high-end essentially unchanged. On the other hand, a decrease of θ_j, 0 (red line) significantly changes the typical energy structure after breakout, making the jet narrower and more powerful, while the cocoon becomes less energetic. The result is a luminosity function which resembles much more that of LGRBs, where the high-end is dominated by on-axis jets and the low-end by off-axis cocoons. With this choice of parameters, the shape resembles more that of the SGRB luminosity function by Ghirlanda et al. (2016).

Fig. 9. Short GRB luminosity function dependence on the parameters of the jet at launch. In each panel, the thick grey line shows the synthetic short GRB luminosity function with fiducial parameters. The thinner, coloured lines show the impact of changing a single parameter. See text for a description of the various dependencies.

Finally, Fig. 10 shows the dependence of the luminosity function on the delay Δt_inj between the onset of the expansion of the ambient medium and the start of the jet injection. Reducing the delay by a factor of 10 (red line) does not change the luminosity functionn significantly, since the fiducial value Δt_inj = 0.1 s is already short compared to the breakout time. Increasing it by the same factor (yellow line) again increases the energy in the cocoon, due to the reduced jet collimation by the ambient medium, which causes the head working surface to widen.

Fig. 10. Short GRB luminosity function dependence on the delay Δtinj between merger and jet injection.

7. Discussion

7.1. Caveats and open questions

While the results of our analysis are encouraging about the feasibility of a GRB unification under the quasi-universal jet scenario, we stress that several caveats remain to be addressed. Here is a non-exhaustive list of issues we think need to be investigated in the future:

– Our model of the jet propagation through the ambient medium, of its collimation by the cocoon, and of the development of structure builds on previous works (Matzner 2003; Bromberg et al. 2011; Lazzati & Perna 2019) which have shown a good agreement with selected numerical simulations. Nevertheless, the model is based on some simplifying assumptions whose impact, to our knowledge, has not been studied in detail yet. To name a few: (i) the jet is assumed to be uniform at launch, but no systematic investigation has been performed, to our knowledge, of how the dynamics and the resulting jet structure would change if the jet were assumed to be already structured at its base (which is a natural expectation in a realistic setting, see e.g. Kathirgamaraju et al. 2019; Christie et al. 2019); (ii) as discussed in Sect. 2, the inclusion of magnetic field might affect the jet dynamics and collimation (but probably not drastically – Levinson & Eichler 2005; Bromberg & Tchekhovskoy 2016). The impact of different degrees of magnetisation and magnetic field configurations on the final jet structure has yet to be addressed;

– There are some indications that the conversion efficiency of jet energy to radiation during the prompt emission might decrease away from the jet axis (e.g. Beniamini et al. 2019; Beniamini & Nakar 2019; Salafia et al. 2019). This could impact the luminosity function;

– In the present work, we assumed all progenitors to be identical (separately for LGRBs and SGRBs). Further studies are necessary to determine how the population would change if we included different progenitors, whether or not this would narrow down even more the required dispersion of jet properties at launch. Another important topic that has not been addressed here is how the jet properties at launch depend on the progenitor.

– An additional fundamental aspect that needs to be addressed is whether it is possible to explain the afterglows of LGRBs and SGRBs with a quasi-universal structure seen at different viewing angles.

Further theoretical work, and hopefully new precious observations from future detections of off-axis jets, will eventually allow us to address these questions and caveats.

7.2. The fraction of choked jets and the duration distribution

During the last few years, some authors (Bromberg et al. 2012; Moharana & Piran 2017; Sobacchi et al. 2017; Petropoulou et al. 2017) explored the idea that the duration distribution of GRBs could contain information about the minimum time needed for the jet to punch through the typical progenitor. The signature of this should be a flattening in the probability distribution of duration around such minimum time. These works assume an intrinsic power-law probability distribution of jet duration

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}T}\propto T^{-\alpha }\end{array}$$(19)

above some minimum duration T_min. Since the jets must spend an average time T_b to break out of the progenitors, the post-breakout duration (T_γ) distribution becomes

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}{T}_{\gamma }}\propto {(T_{\mathrm{b}}+T_{\gamma })}^{-\alpha }\end{array}$$(20)

which reflects the intrinsic duration distribution $\propto T_{\gamma }^{-\alpha }$ when T_γ ≫ T_b, but it is flat when T_γ ∼ T_b. The model of Petropoulou et al. (2017) includes a power-law distribution of luminosity, which makes then T_b luminosity-dependent, but the results are similar. These models reproduce well the observed GRB duration distribution. Since the values of α that reproduce the observed population are large (α ∼ 4), they predict a large fraction of choked jets (the fraction also depends on the minimum duration T_min and, in the model of Petropoulou et al. 2017, also on the minimum luminosity, but a small fraction of choked jets would require a fine tuning of T_min close to T_b). In contrast with that, our synthetic populations are dominated by jets that successfully break out from their progenitors (Sects. 6.2.2 and 6.3.2). It is natural, therefore, to ask whether the duration distribution of our populations does reflect the actual observed one. Figure 11 shows a comparison between the observed duration probability distribution and that of our synthetic sample. The blue histogram shows the estimated probability distribution of rest-frame duration (defined as T₉₀/(1 + z)) of the sample of Swift GRBs with measured redshift. For each bin, we compute the probability estimate as

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}T}(T_i)\approx \frac{N_i}{N_{\mathrm{tot}}\mathrm{\Delta }{T}_i},\end{array}$$(21)

Fig. 11. GRB duration probability distribution. The blue histogram shows the estimated probability distribution of intrinsic duration (defined as T90/(1 + z)) of Swift-detected GRBs with a measured redshift. Red points represent the corresponding distribution for our simulated populations (assuming an observed rate ratio of 13.5 between LGRB and SGRB, i.e. the same rate ratio as in the Swift sample), while the grey points show the distribution for the LGRB synthetic population alone (downward arrows are upper limits that correspond to empty bins). Vertical error bars show the Poisson uncertainty associated to each bin.

where T_i is the central duration of the bin, ΔT_i is the bin width, N_i is the number of events in that bin, and N_tot is the total number of events. The associated (Poisson) uncertainty is

$$\begin{array}{c}\hfill \delta \frac{\mathrm{d}p}{\mathrm{d}T}(T_i)\approx \frac{\sqrt{N_i}}{N_{\mathrm{tot}}\mathrm{\Delta }{T}_i}·\end{array}$$(22)

For empty bins, we take as upper limit the estimate obtained by assuming a single event in that bin. Red points show the corresponding probability distribution (for the duration defined in Sect. 6.1) of our synthetic samples, where we assumed a 13.5 observed rate ratio of LGRB over SGRB (i.e. the same as for the Swift sample used). Grey points show the LGRB synthetic sample alone. The comparison shows a good agreement, indicating that the observed GRB duration distribution is still consistent with the possibility that most GRB jets successfully break out of their progenitor. The observation of a successful GRB jet associated to GW170817 (Ghirlanda et al. 2019; Mooley et al. 2018a) suggests that this is the case for SGRBs (see Beniamini et al. 2019, for an excellent discussion of this latter point).

8. Conclusions

During the last two decades, evidence accumulated in favour of the presence of jet structure in both LGRB and SGRB, unveiling its important role in shaping the appearance of these sources. The limited range of properties of their putative progenitors seems to suggest the possibility that such jet structure is quasi-universal: this would have a great unifying impact. In this work, we set up a model of the interaction of a GRB jet with its environment, specifically aimed at predicting the jet structure after breakout. With this model, we computed self-consistently the quasi-universal structure that would result if the properties of GRB jets at launch and those of their progenitors were distributed in very narrow ranges. Strikingly, these quasi-universal structures can reproduce both the LGRB and SGRB luminosity functions. While we acknowledge that several issues remain to be addressed, we consider this results is encouraging towards the unification of GRBs within the quasi-universal hypothesis.

Acknowledgments

We thank the anonymous referee for raising important points that helped to improve the quality of this work. We thank R. Ciolfi for the discussion that sparked the research behind this work. We thank G. Ghirlanda and G. Ghisellini for stimulating discussions and for their continued support and motivation. We thank D. Lazzati and B. Giacomazzo for insightful comments on the first draft. O.S. acknowledges the MIUR (Italian Ministry of University and Research) grant “FIGARO” (1.05.06.13) for support. S. A. acknowledges the GRAvitational Wave Inaf TeAm – GRAWITA (P.I. E. Brocato) for support.

References

Appendix A: Comparison with simulations

In this section we compare our semi-analytical model with results of a number of numerical relativistic hydrodynamical simulations described in the literature.

A.1. Jet head propagation

Our modelling of the jet head propagation follows quite closely previous works (Matzner 2003; Bromberg et al. 2011), with some adjustment to account for a moving ambient medium. Similar models have been presented recently by Matsumoto & Kimura (2018), Gill et al. (2019b) and Lazzati & Perna (2019), with some minor differences. We investigate here the validity of this approach by comparing the breakout time, breakout radius and cocoon energy predicted by our model to those obtained by Gottlieb et al. (2018a, G18 hereafter) and Nagakura et al. (2014, N14 hereafter) in numerical simulations. Table A.1 shows the results of our comparison.

Configurations A, B and C in G18 all involve a relativistic jet of luminosity L_j = 2 × 10⁵⁰ erg s⁻¹ and initial opening angle θ_j, 0 = 10° launched in an ambient medium representing binary neutron star ejecta in homologous expansion, with a total mass M_ej = 10⁻² M_⊙ and a power-law density profile decreasing as r^−3.5. The configurations differ by the maximum velocity β_max of the ejecta (being 0.4 in configuration C and 0.2 in the other two) and by the delay t_inj between the start of the ejecta expansion and the jet injection (t_inj = 80, 240 and 40 ms for A, B and C respectively). We reproduced the ambient medium configuration and followed the jet evolution semi-analytically using our model, up to the jet breakout. The breakout times and cocoon energies predicted by our model agree with the simulations within 20%.

The initial configurations of N14 are very similar to those of G18. Their reference configuration (M-ref) is the same as case C of G18, but with t_inj = 50 ms and θ_j, 0 = 15°. The other configurations, listed in Table A.1, are variants over the reference one, exploring alternative values of some parameters, namely L_j = 4 × 10⁵⁰ erg s⁻¹ (M-L4), θ_j, 0 = 30° (M-th30) and 45° (M-th45), t_inj = 500 ms (M-ti500), and M_ej = 10⁻³ M_⊙ (M-M3), 2 × 10⁻² M_⊙ (M-M2-2) and 10⁻¹ M_⊙ (M-M1). Also in this case, the results of our model agree with the simulations within ∼20%, but they show a systematic trend towards faster breakouts, which may indicate that the $\stackrel{\sim }{L}$ correction adopted from Harrison et al. (2018) needs to be modified in the case of a moving ambient medium.

A.2. Collimation and jet opening angle

We performed an additional comparison with the time evolution of the average jet opening angles as reported in N14. Figure A.1 shows the time evolution of the height-averaged opening angle of the jet from N14 (bottom panel) and as predicted by our model (top panel). The angle was computed as

$$\begin{array}{c}\hfill {\theta }_{\mathrm{av}}(t)=\frac{{\int }_{z_0}^{z_{\mathrm{h}}(t)}{\theta }_{\mathrm{j}}(z,t)\mathrm{d}z}{z_{\mathrm{h}}(t)-z_0}·\end{array}$$(A.1)

Fig. A.1. Comparison between the height-averaged opening angle evolution from the numerical simulations in Nagakura et al. 2014 (right-hand panel, adapted from their Fig. 3) and those predicted by our semi-analytical model (left-hand panel). Each line refers to a different initial configuration (see text). Stars and upward triangles in the left-hand panel mark tbo and tbo + tdelay (see Sect. 3), respectively. The evolution for configuration M-th45 is not shown in the original figure.

In the simulations, θ_j(z, t) is the opening angle containing the “relativistic” (hΓ >  10) material located at a height z, at time t. In our model, according to our assumptions detailed above, the jet is conical up to $Z=(\widehat{z}+z_c)/2$ and cylindrical above that point, so that in our case

$$\begin{array}{c}\hfill {\theta }_{\mathrm{av}}(t)=\frac{{\theta }_{\mathrm{j},0}Z(t)}{z_{\mathrm{h}}(t)-z_0}[1-\frac{z_0}{Z(t)}+ln\left(\frac{z_{\mathrm{h}}(t)}{Z(t)}\right)]\end{array}$$(A.2)

before breakout (t ≤ t_bo), and

$$\begin{array}{c}\hfill {\theta }_{\mathrm{av}}(t)=\frac{{\theta }_{\mathrm{j}}(t)z_{\mathrm{bo}}}{z_{\mathrm{bo}}-z_0}[1-\left(\frac{{\theta }_{\mathrm{j},0}}{{\theta }_{\mathrm{j}}(t)}\right)\frac{z_0}{z_{\mathrm{bo}}}+ln\left(\frac{{\theta }_{\mathrm{j},0}}{{\theta }_{\mathrm{j}}(t)}\right)]\end{array}$$(A.3)

after breakout (t >  t_bo). Figure A.1 shows a remarkably good agreement (given the simplifications adopted) between our model and the simulations in most cases, after setting α = 1/8 in Eq. (10). With the initial configuration M-th45 (larger jet) Nagakura et al. (2014) find that their jet fails to penetrate the ejecta within their simulation time 1 s, so they do not report the breakout time and radius. The imperfect agreement in case M-th30 suggests that the small angle approximations adopted in the treatment of collimation may break down at jet opening angles as large as ≳30°.

A.3. Jet and cocoon structure

In order to compare our recipe for the jet structure to numerical simulations, we reproduced the charachteristics of the ambient media and the properties of the jet at its base from the simulations described in Lazzati et al. (2017b) and Xie et al. (2018). We then applied our semi-analytical model of the jet head propagation (Sect. 2) to find the breakout quantities. Finally, we used the latter as input to compute the jet and cocoon structures as detailed in Sect. 3, Sect. 4 and Sect. 4.2. Let us describe the three configurations and the results of the comparison.

Lazzati et al. (2017a,b). Lazzati et al. (2017b) used a modified version of the relativistic hydrodynamics code FLASH (Fryxell et al. 2000) to simulate a jet propagating and breaking out of an ambient medium with a static, spherically symmetric density distribution with a profile ρ(r) ∝ r²exp(−r). The jet is injected within an angle θ_j, 0 = 16° with a constant luminosity L_j = 10⁵⁰ erg s − 1 and lasts for T_jet = 1 s. Even though the authors state that the terminal Lorentz factor of the injected jet is Γ₃₀₀, the actual maximum Lorentz factor attained in the simulation is lower, presumably because its evolution was not followed up to the coasting phase. We thus choose to set Γ_j = 100. The resulting jet structure is shown in Fig. A.2. The general agreement is good, especially for what concerns the kinetic energy density distribution. The Lorentz factor profile is in quite good agreement in the inner jet and in the outer cocoon, but it falls off too fast in the transition from the jet to the cocoon.

Fig. A.2. Comparison between the jet structure predicted by our model and the results of numerical simulations by Lazzati et al. (2017b). Upper panel: distribution of kinetic energy per unit solid angle, while Lower panel: average Lorentz factor diminished by one. Grey dashed and dotted lines represent respectively the jet and the cocoon as computed by our model (Sects. 3, 4 and 4.2). Red lines show the combined structure. Blue thick lines show the results of the numerical simulations (the Lorentz factor profile is given in Lazzati et al. 2018).

Xie & MacFadyen (2019). The simulations in Xie et al. (2018) model two jets propagating in a moving ambient medium which represents the result of the merger of two neutron stars. The medium is in homologous expansion and features an inner denser cloud surrounded by a fast tail with a steep density fall-off. The two cases modelled, dubbed “narrow engine” and “wide engine”, differ by the opening angle of the jet at its base, and by the jet enthalpy at injection. One subtlety is that the jets in the two models are injected following different techniques: in the wide engine case, the jet is injected through a cylindrical nozzle and thus has a quite well-defined opening angle θ_j, 0 = 0.35 rad; in the narrow engine case, the jet is injected through a circular nozzle as in Duffell & MacFadyen (2015), and features a non-uniform angular distribution of luminosity, with an approximately Gaussian form $\underset{\sim }{\propto }\exp [-{(\theta /{\theta }_{\text{0}})}^2/2]$, with θ₀ = 0.1 rad. We thus tried different values of the base jet half-opening angle θ_j, 0 in our model between θ₀ and 2θ₀, looking for the best agreement with their results. We found that θ_j, 0 = 1.38θ₀ produces a very good agreement in both the kinetic energy and Lorentz factor profile (setting Γ_j = 100, that is, equal to their terminal Lorentz factor), as shown in Fig. A.3. In the wide engine case this tuning is not necessary, and we simply adopt the value θ_j, 0 = 0.35 rad as in the simulation, and we set Γ_j = 20, again equal to their terminal Lorentz factor. The result is shown in Fig. A.4. In this latter case, the agreement in the Lorentz factor profile is not as good.

Fig. A.3. Same as Fig. A.2, but for the “narrow engine” case in Xie et al. (2018) (we choose the structure extracted from the 104 s snapshot as shown in their Fig. 2, right-hand panels).

Fig. A.4. Same as Fig. A.2, but for the “wide engine” case in Xie et al. (2018) (we choose the structure extracted from the 104 s snapshot as shown in their Fig. 7, right-hand panels).

These comparisons show that, while s omewhat crude, our recipe seems to capture the main trends in the development of the jet and cocoon structure, at least in the conditions covered by these simulations. We plan to run a series of dedicated numerical hydrodynamical simulations to investigate further these aspects.

All Tables

All Figures

Fig. 1. Sketch of the various components of the jet propagation model, with some key quantities annotated. All quantities are defined and described in the text. Typical radii (order of magnitude in cm) are shown in red on the right, corresponding to the progenitors and average jet properties we employ when computing our synthetic jet populations (note that the collimation-related radii zc, zcoll and $\widehat{z}$ change during the jet propagation).

In the text

Fig. 2. Radial density profile from our representative LGRB progenitor (model 16TI of Woosley & Heger 2006), which represents a low-metallicity, high-angular momentum, massive star prior to collapse. The low-density power-law tail, visible in the lower-right corner of the figure, has been added by ourselves to the original density profile of the model. It represents a stellar wind with typical Wolf-Rayet mass loss rate Ṁw = 10−5 M⊙ yr−1 and velocity vw = 1000 km s−1.

In the text

Fig. 3. Distribution of jet structures in our LGRB population. In each panel, at each fixed angle θ, the yellow (resp. orange) filled area encloses 90% (resp. 50%) of the structures in the population, while the red line shows the median of the distribution at that angle. Top panel: energy structures (kinetic energy per unit solid angle as a function of the angle θ from the jet axis). Bottom panel: velocity structures (average Lorentz factor diminished by one, as a function of θ).

In the text

	Fig. 4. LGRB luminosity distribution of our model population (thick solid red histogram) compared to the observed distribution (black points, from Pescalli et al. 2015). Thinner histograms show the luminosity distributions of three sub-classes corresponding to three viewing angle bins, reported in the legend.
In the text

	Fig. 5. Long GRB luminosity function dependence on the parameters of the jet at launch. Each panel shows how the fiducial luminosity function (thick grey solid line) is affected by a change in one of the parameters (coloured thinner lines). See text for an interpretation of the dependencies.
In the text

Fig. 6. Density and velocity profile of the SGRB ambient medium model. The red line shows a snapshot of the density profile (along the polar axis) of the neutron star merger ejecta model described in Xie et al. (2018) (also used in Kasliwal et al. 2017 and Gottlieb et al. 2018b). The blue dashed line shows the corresponding velocity profile (velocity axis on the right). The ejecta are assumed to expand homologously.

In the text

	Fig. 7. Distribution of structures in our SGRB population. The filled areas and the red lines have the same meaning as in Fig. 3. The black lines show the confidence regions (dotted: two sigma, dashed: one sigma, solid: best fit) for the jet structure of GRB 170817A from Ghirlanda et al. (2019).
In the text

Fig. 8. SGRB luminosity distribution of our model population (thick solid red histogram) compared to two luminosity function models based on observations (Wanderman & Piran 2015 – green solid line and shaded area – and Ghirlanda et al. 2016 – purple solid line and shaded area). Thinner histograms show the luminosity distributions of sub-classes of jets belonging to three viewing angle bins, reported in the legend. The grey line and filled horizontal band show the binary neutron star merger local rate based on gravitational wave observations (The LIGO Scientific Collaboration 2019).

In the text

	Fig. 9. Short GRB luminosity function dependence on the parameters of the jet at launch. In each panel, the thick grey line shows the synthetic short GRB luminosity function with fiducial parameters. The thinner, coloured lines show the impact of changing a single parameter. See text for a description of the various dependencies.
In the text

	Fig. 10. Short GRB luminosity function dependence on the delay Δtinj between merger and jet injection.
In the text

Fig. 11. GRB duration probability distribution. The blue histogram shows the estimated probability distribution of intrinsic duration (defined as T90/(1 + z)) of Swift-detected GRBs with a measured redshift. Red points represent the corresponding distribution for our simulated populations (assuming an observed rate ratio of 13.5 between LGRB and SGRB, i.e. the same rate ratio as in the Swift sample), while the grey points show the distribution for the LGRB synthetic population alone (downward arrows are upper limits that correspond to empty bins). Vertical error bars show the Poisson uncertainty associated to each bin.

In the text

Fig. A.1. Comparison between the height-averaged opening angle evolution from the numerical simulations in Nagakura et al. 2014 (right-hand panel, adapted from their Fig. 3) and those predicted by our semi-analytical model (left-hand panel). Each line refers to a different initial configuration (see text). Stars and upward triangles in the left-hand panel mark tbo and tbo + tdelay (see Sect. 3), respectively. The evolution for configuration M-th45 is not shown in the original figure.

In the text

Fig. A.2. Comparison between the jet structure predicted by our model and the results of numerical simulations by Lazzati et al. (2017b). Upper panel: distribution of kinetic energy per unit solid angle, while Lower panel: average Lorentz factor diminished by one. Grey dashed and dotted lines represent respectively the jet and the cocoon as computed by our model (Sects. 3, 4 and 4.2). Red lines show the combined structure. Blue thick lines show the results of the numerical simulations (the Lorentz factor profile is given in Lazzati et al. 2018).

In the text

	Fig. A.3. Same as Fig. A.2, but for the “narrow engine” case in Xie et al. (2018) (we choose the structure extracted from the 104 s snapshot as shown in their Fig. 2, right-hand panels).
In the text

	Fig. A.4. Same as Fig. A.2, but for the “wide engine” case in Xie et al. (2018) (we choose the structure extracted from the 104 s snapshot as shown in their Fig. 7, right-hand panels).
In the text