A THEORY OF MULTICOLOR BLACKBODY EMISSION FROM RELATIVISTICALLY EXPANDING PLASMAS

Relativistic outflows, often in the forms of jets, are a common phenomenon in many astronomical objects, such as microquasars (Mirabel & Rodriguez 1994; Hjellming & Rupen 1995), active galactic nuclei (AGNs; Lind & Blandford 1985; Gopal-Krishna et al. 2007), and gamma-ray bursts (GRBs; Paczyński 1986; Goodman 1986). In many of these objects, the density at the base of the flow is sufficiently high, so that the optical depth to Thomson scattering by the baryon-related electrons is much larger than unity. As a result, if a source of photons exists deep enough in the flow, the emerging spectrum is inevitably thermal or quasi-thermal (a Wien spectrum could also emerge if the number of photons is conserved by the radiative processes). These photons escape the flow once they decouple from the plasma, at the photosphere (e.g., Paczyński 1990).

As a result of the relativistic expansion of the source, even if the emitted spectrum (in the plasma comoving frame) is purely thermal—i.e., it is not modified by any additional, non-thermal radiative process—still the observed spectrum is not expected to follow the Planck distribution. This results from the aberration of light in the expanding plasma, and is of pure geometrical nature. It is therefore an inherent property of any relativistically expanding photon emitting source. The origin of this effect lies in the fact that at any given instance, an observer simultaneously sees photons that emerge the expanding plasma from a range of radii and angles. Therefore, each (thermal) photon has its own comoving energy, and is seen with a particular Doppler shift. The resulting, integrated spectrum is non-thermal.

Analysis of this effect begins with the non-trivial shape of the photosphere in a relativistically expanding plasma source. By definition, the photosphere is a surface in space that fulfills the following requirement: the optical depth to scattering a photon originating from a point on this surface and reaching the observer is equal to unity. Although it is mathematically a two-dimensional surface in space, for spherically symmetric wind, as is considered here, this surface is symmetric with respect to rotation around the axis to the line of sight. It is therefore appropriate to refer to it as the "photospheric radius," which is a function of the angle to the line of sight, r_ph = r_ph(θ).

Calculation of the photospheric radius for the scenario of a steady, spherically symmetric, relativistic wind characterized by a constant Lorentz factor Γ was first carried by Abramowicz et al. (1991) and later extended by Pe'er (2008). For relativistic winds, Γ ≫ 1, and small angles to the line of sight, θ ≪ 1, the calculation results in a simple, yet non-trivial dependence of the photospheric radius on the angle to the line of sight: r_ph(θ)∝(Γ⁻² + θ²/3), where the proportionality constant depends on the properties of the outflow. It is thus found that in relativistic expanding wind, the photospheric radius has a strong dependence on the angle to the line of sight, a fact which leads to several non-trivial consequences: for example, it was shown by Pe'er (2008) that for parameters characterizing GRBs, photospheric emission can be expected up to tens of seconds (albeit with a decreasing flux; see further discussion below).

By definition, the photospheric radius provides only a first-order approximation to the last scattering position (=decoupling position) of the thermal photons. This is because, in principle, photons can be scattered at any point in space in which electrons exist. Therefore, a full description of the last scattering position and scattering angle can only be done in terms of probability density function P(r, θ). As was shown by Pe'er (2008) and will be further discussed here, use of this function is essential in calculating the spectrum and flux. The shape of the observed spectrum, for example, depends on the individual Doppler shifts of the observed photons, but also on the number of photons that undergo the particular Doppler shift—information that is only held in the probability density function.

The probability density function can be considered as an extension of the standard use of the photospheric radius. Instead of considering a surface in space from which thermal photons emerge, one considers the entire space, weighted by the finite probability of a photon to emerge from an arbitrary radius r and arbitrary angle θ. Using this function, Pe'er (2008) calculated the expected temporal decay laws of the observed thermal photon flux and average temperature, following an abrupt termination of the inner engine. It was shown there that at late times the thermal flux is expected to decay as F∝t⁻² and the average temperature decays as 〈T〉 ∝t^−2/3. The agreement found between these theoretical predictions and the late time decay laws of the peak energy and flux observed during GRB prompt emission (Ryde 2004, 2005; Ryde & Pe'er 2009), is one of the key motivations in studying the properties of thermal emission in the context of the GRB prompt emission phase.

On a more general ground, thermal emission may be crucial in understanding the nature of GRB prompt emission. In spite of being studied for nearly two decades now, the origin of the prompt emission in GRBs is still puzzling. In recent years, it became clear that synchrotron emission, the leading emission model (Rees & Mészáros 1994; Tavani 1996; Cohen et al. 1997; Sari & Piran 1997; Panaitescu et al. 1999) cannot account for the steepness of the low-energy spectral slopes seen (Crider et al. 1997; Preece et al. 1998, 2002; Ghirlanda et al. 2003) (see, however Bošnjak et al. 2009). This motivated some alternative ideas, such as reprocessing through a heated cloud (Dermer & Böttcher 2000), jitter radiation (Medvedev 2000), or decaying magnetic field (Pe'er & Zhang 2006). Moreover, in order to account for the high efficiency of the prompt emission seen (Zhang et al. 2007; Nysewander et al. 2009) using the synchrotron model, highly efficient energy dissipation is required, which is difficult to be accounted for in the classical internal shocks scenario (Kobayashi et al. 1997; Daigne & Mochkovitch 1998; Lazzati et al. 1999; Guetta et al. 2001). Contribution from thermal emission thus seems a natural way of overcoming both these issues. First, photospheric emission is inherent to the fireball model (Eichler & Levinson 2000; Mészáros & Rees 2000; Mészáros et al. 2002; Daigne & Mochkovitch 2002; Rees & Mészáros 2005). Second, as it does not originate from any internal dissipation, the contribution from thermal photons reduces the efficiency requirement (Ryde & Pe'er 2009). Finally, as will be shown here, photospheric emission is able to produce low-energy spectral slopes that are consistent with those observed.

In spite of the success of the thermal emission model in reproducing the late time decay of the peak energy and flux seen in many GRBs (Ryde & Pe'er 2009), the simple version of the model suffers several drawbacks. One issue that is often raised, is that the low-energy spectral slopes seen during the prompt phase of many GRBs are too shallow to be accounted for by the Rayleigh–Jeans tail of the thermal spectrum (e.g., Bellm 2010). Indeed, the Rayleigh–Jeans tail implies a spectral slope F_ν∝ν², while GRB observations show that on the average, the low-energy photon index in the "Band" function fits is α ≃ −1 (i.e., F_ν∝ν⁰; see Preece et al. 2000; Kaneko et al. 2006, 2008).

As discussed above, the observed spectrum resulting from a photospheric emission in relativistically expanding plasma does not necessarily need to be a pure blackbody, but should, in general, be modified. Therefore, a pure blackbody spectrum will in many cases not be able to fit the observed spectrum. The main purpose of this paper is to calculate the observed spectrum resulting from photospheric emission in a scenario of a steady, relativistic outflow. In fact, as we show below, at late times following the termination of the inner engine, the resulting spectrum is expected to be close to a power law below the thermal peak, with power-law index F_ν∝ν⁰. We point out that this result is remarkably similar to the average value of the low-energy spectral slopes seen in large samples of GRBs (Kaneko et al. 2006, 2008). While this is very encouraging, we point below to several practical limitations that limit the validity of our result. While therefore we cannot claim that the problem is fully solved, the results presented here, which show the important role of the geometrical aberration of the light in determining the spectrum, may be significant in understanding the low-energy spectral slopes of GRBs.

This paper is organized as follows. In Section 2, we briefly discuss the properties of the photosphere, and the characteristic time scales up to which thermal emission is expected. In Section 3, we describe the construction of the probability density function. These sections closely follow the treatment by Pe'er (2008), and are given here for completeness. We then calculate in Section 3.1 the observed spectrum at late times, and show that the energy spectrum can be approximated as F_ν(t)∝ν⁰t⁻². We compare the analytical predictions to the numerical results in Section 4. We then discuss the implications of our results on the observed GRB prompt emission spectra in Section 5.

Consider the ejection of a spherically symmetric plasma wind from a progenitor characterized by constant mass-loss rate $\dot{M}$ that expands with time independent velocity v = βc. The ejection begins at t = 0 from radius r = 0, thus at time t the plasma outer edge is at radius r_out(t) = βct from the center. However, here we assume that the plasma wind occupies the entire space, i.e., r_out(t) → ∞.⁵ For constant ${\dot{M}}$ and Γ, at r < r_out the comoving plasma density is given by $n^{\prime }(r) = {\dot{M}}/(4\pi m_p v \Gamma r^2)$, where Γ = (1 − β²)^−1/2. We further assume that emission of photons occurs deep inside the flow where the optical depth is τ ≫ 1, as a result of unspecified radiative processes. The emitted photons are coupled to the flow (e.g., via Compton scattering), and are assumed to thermalize before escaping the plasma once the optical depth becomes low enough.

Under these assumptions, it was shown by Pe'er (2008) that the optical depth of a photon emitted at radius r, angle to the line of sight θ and propagates toward the observer is

$$\begin{equation} \tau (r,\theta) = {R_d \over \pi r}\left[{\theta \over \sin (\theta)} - \beta \right] \simeq {R_d \over 2 \pi r} \left({1\over \Gamma ^2} +{\theta ^2 \over 3}\right). \end{equation} \tag{ 1 }$$

Here,

$$\begin{equation} R_d \equiv {{\dot{M}} \sigma _T \over 4 m_p \beta c}, \end{equation} \tag{ 2 }$$

m_p is the proton rest mass and σ_T is Thomson cross section. The last equality in Equation (1) holds for Γ ≫ 1 and small angle to the line of sight, θ ≪ π/2, which allows the expansion sin (θ) ≃ θ − θ³/6.

The photospheric radius is obtained by setting τ(r_ph, θ) = 1,

$$\begin{equation} r_{{\rm ph}}(\theta) \simeq {R_d \over 2 \pi } \left({1\over \Gamma ^2} +{\theta ^2 \over 3}\right). \end{equation} \tag{ 3 }$$

Equation (3) implies that for small viewing angle, θ ≪ Γ⁻¹ the photospheric radius is angle independent, r_ph ≃ R_d/2πΓ², while for large angles θ ≫ Γ⁻¹, the photospheric radius is r_ph(θ) ≃ R_dθ²/6π.

Characteristic observed times. Below the photosphere, the photons are coupled to the flow, therefore their effective propagation velocity in the radial direction is similar to the outflow velocity, ≈βc. Assuming that a photon is emitted at t = 0, r = 0, it decouples the plasma at time t = r/βc.⁶ Consider photons that propagate toward the observer at angle to the line of sight θ. These photons are observed at a time delay with respect to a hypothetical photon that was emitted at t = 0, r = 0 and did not suffer any time delay ("trigger" photon), which is given by

$$\begin{equation} \Delta t^{{\rm ob}.}(r,\theta) = {r \over \beta c}\times [1 - \beta \cos (\theta)]. \end{equation} \tag{ 4 }$$

For relativistic outflows, Γ ≫ 1, thermal photons emitted from the photospheric radius r_ph(θ) on the line of sight (θ = 0), are seen at a very short time delay with respect to the trigger photon, which is given by

$$\begin{equation} \Delta t^{{\rm ob}.}(r_{{\rm ph}},\theta =0)\equiv t_N \simeq {R_d \over 4 \pi \Gamma ^4 c} \simeq 10^{-2} \;L_{52}\; \Gamma _{2}^{-5} \; {\rm s}. \end{equation} \tag{ 5 }$$

Here, ${\dot{M}} = L/\Gamma c^2$, and typical parameters characterizing GRBs, L = 10⁵² L₅₂  ergs⁻¹ and Γ = 100 Γ₂ were used.

On the other hand, due to the strong angular dependence of the photospheric radius, thermal photons emitted from the photosphere with high angles to the line of sight, θ ≫ Γ⁻¹ (and θ ≪ 1) are observed at a much longer time delay,

$$\begin{equation} t^{{\rm ob}.}(\theta \gg \Gamma ^{-1}) \simeq {R_d \over 3 \pi c}\left({\theta ^2 \over 2}\right)^2 \simeq 30 \; L_{52}\; \Gamma _{2}^{-1} \theta _{-1}^4 {\rm \; s}, \end{equation} \tag{ 6 }$$

where θ = 0.1 θ₋₁ rad. The time scale derived on the right-hand side of Equation (6) is based on the estimate of the jet opening angle in GRB outflow, θ ⩽ θ_j ≃ 0.1 rad (e.g., Berger et al. 2003). Equation (6) therefore shows that in a relativistically expanding wind with parameters that can characterize GRBs, thermal emission is expected up to tens of seconds following the decay of the inner engine.

Any attempt of describing the photospheric emission must consider the fact that photons have a finite probability of being scattered at every point in space in which electrons exist (i.e., in the entire space, and not only on the photospheric surface). This led Pe'er (2008) to introduce the concept of a probability density function as a mathematical tool which is needed in calculating the observed properties of the photospheric emission, for example, the temporal decay of the temperature and flux. Since, as we show here, this is a fundamental concept that is necessary in an analysis of the properties of the photospheric emission, and in particular calculation of the expected spectrum, we briefly repeat in this section the basic definition and use of this function, before calculating the spectrum in Section 3.1.

The thermal photons are advected with the flow below the photosphere, until the last scattering event (the decoupling) takes place, at some radius r. For every radius r, there is an associated probability that the last scattering event occurs at that particular radius (see below). During the last scattering event, a photon is scattered into angle θ. An underlying assumption is that photons that decouple from the plasma at radius r and scattered into angle θ are observed at a delay given by Equation (4). This assumption implies that (1) the delay time of a photon is solely determined by two parameters, the last scattering radius and scattering angle (for constant outflow velocity), and (2) at any given instance, an observer simultaneously sees photons emitted from a range of radii and angles, all fulfilling the requirement set by Equation (4).

The θ-dependence of the photospheric radius implies that the probability of a photon to be scattered into angle θ depends on the radius at which the scattering event takes place (or vice versa). However, here we assume that the probabilities are independent, i.e., P(r, θ) = P(r) × P(θ). This assumption is made in order to simplify the calculation, and is tested against the numerical results (see Section 4). While clearly this approximation has only a limited validity, the results obtained are in good agreement with the precise calculation done numerically. We further discuss this approximation, as well as its limitations in Section 5.

Equation (1) implies that the optical depth to scattering depends on the radius as τ(r)∝r⁻¹. This optical depth is the integral over the probability of a photon propagating from radius r to +∞ to be scattered, i.e., τ(r) = ∫^∞_r(dτ/dr)dr, from which it is readily found that (dτ/dr)|_r∝r⁻². As a photon propagates in the radial direction from radius r to r + δr, the optical depth in the plasma changes by δτ = (dτ/dr)|_rδr. Therefore, the probability of a photon to be scattered as it propagates from radius r to r + δr is given by

$$\begin{equation} P_{{\rm sc}.}(r \cdots r + \delta r) = 1 - e^{-\delta \tau } \approx \delta \tau \propto {\delta r \over r^2}. \end{equation} \tag{ 7 }$$

For the last scattering event to take place at r⋅⋅⋅r + δr, it is required that the photon does not undergo any additional scattering before it reaches the observer. The probability that no additional scattering occurs from radius r to the observer is given by exp (− τ[r]). The probability density function P(r) for the last scattering event to occur at radius r is therefore written as

$$\begin{equation} P(r) = {r_0 \over r^2} e^{-(r_0/r)}. \end{equation} \tag{ 8 }$$

The function P(r) in Equation (8) is normalized, ∫^∞₀P(r)dr = 1. Comparison to Equation (1) gives the proportionality constant, r₀ ≡ r_ph(θ = 0) = R_d/2πΓ².

The probability of a photon to be scattered into angle θ is calculated assuming isotropic scattering in the comoving frame, i.e., dσ/dΩ' = Const.⁷ The comoving spatial angle is dΩ' = sin θ'dθ'dϕ', and therefore the probability of a photon to be scattered to angle θ' (in the comoving frame) is dP/dθ'∝sin θ'. The proportionality constant is obtained by integrating over the range 0 ⩽ θ' ⩽ π, and is equal to 1/2. Thus, the isotropic scattering approximation leads to P(θ') = (sin θ')/2.

Assuming that on the average, photons propagate in the radial direction, by making Lorentz transformation to the observer frame one obtains the probability of scattering into angle θ with respect to the flow direction. This angle is equal to the observed angle to the line of sight, and is given by

$$\begin{equation} P(\theta) = P(\theta ^{\prime }) {d\theta ^{\prime } \over d\cos \theta ^{\prime }} {d \cos \theta ^{\prime } \over d \cos \theta } {d \cos \theta \over d\theta } = {\sin \theta \over 2 \Gamma ^2 (1-\beta \cos \theta)^2}. \end{equation} \tag{ 9 }$$

Defining u ≡ 1 − βcos θ, Equation (9) becomes

$$\begin{equation} P(u) = { 1 \over 2 \Gamma ^2 \beta u^2}. \end{equation} \tag{ 10 }$$

Note that 1 − β ⩽ u ⩽ 1 + β, and the function P(u) in Equation (10) is normalized, ∫^{1 + β}_{1 − β}P(u)du = 1.

3.1. Spectrum and Decay Law of the Thermal Flux at Late Times

As long as the radiative processes that produce the thermal photons deep inside the flow are active, the observed thermal radiation is dominated by photons emitted on the line of sight toward the observer. Once these radiative processes are terminated, the radiation becomes dominated by photons emitted off-axis and from larger radii, which determine the late time behavior of the spectrum and flux. Therefore, the limiting case of a δ-function injection, both in time and radius (t = 0, r = 0) is expected to closely describe the late time behavior of the thermal spectrum. We calculate here the observed spectrum and flux of the thermal emission at late times, under these assumptions.

Denote by T'(r) the photon comoving temperature, its observed temperature is $T^{{\rm ob}} = T^{\prime }(r) {\mathcal {D}}$, where ${\mathcal {D}}= [\Gamma (1-\beta \cos \theta)]^{-1} = (\Gamma u)^{-1}$ is the Doppler factor. The observed photon temperature therefore depends on the viewing angle as well as on the radius of decoupling. Below the photospheric radius, photons lose their energy adiabatically, hence T'(r)∝r^−2/3. It was shown by Pe'er (2008) that a similar decay law for the photon temperature exists even if the energy density in the photon field is much smaller than the energy density in the electrons (which can in principle be non-relativistic in the comoving frame, hence have a different temperature decay law), resulting from the spatial three-dimensional expansion of the plasma. The adiabatic losses take place only as long as the photons propagate at radii smaller than ∼few × r₀. While photons that propagate at high angles decouple the plasma at much larger radii than r₀ (see Equation (3)), above few × r₀, the number of scattering is small, and hence the photons maintain their energy. As here we are interested in the late time evolution of the spectrum and flux, where late time imply t^ob ≫ t_N,⁸ we can safely assume that the comoving temperature of photons that dominate the flux at late times is (on the average) constant, T'(r) = Const.

Assume that N₀ photons are emitted instantaneously (a δ-function injection in time) at the center of the expanding plasma. Each photon is advected with the flow, until the last scattering event takes place at radius r and into angle θ, after which it propagates freely. The observed flux density (or differential energy flux) from a source at luminosity distance d_L is therefore given by

$$\begin{eqnarray} F_{\nu }(t^{{\rm ob}}) &\equiv& {d^2 F \over dt d\nu } = {N_0 \over 4 \pi d_L^2} \int P(r) dr\int P(u) du k_B T^{{\rm ob}}(r,u) \nonumber\\ &&\times \delta \left(t^{{\rm ob}} = {r u \over \beta c} \right) \delta \left(\nu = {k_B T^{\prime }(r) \over h \Gamma u} \right). \end{eqnarray} \tag{ 11 }$$

Here, F is the total fluence, h is Planck's constant, and k_B is Boltzmann's constant. Note that we assumed here that the photon observed frequency ν directly corresponds to the observed temperature, ν = k_BT^ob/h.

As discussed above, at late times, t^ob > t_N(ν_max/ν), one can write ν = k_BT'₀/hΓu, where T'₀ is r-independent. Using P(r) and P(u) from Equations (8) and (10), and the identities δ(t^ob = ru/βc) = δ(r = βct^ob/u) × βc/u and δ(ν = k_BT'₀/hΓu) = δ(u = k_BT'₀/Γhν) × hΓu²/k_BT'₀, Equation (11) becomes

$$\begin{eqnarray} F_{\nu }(t^{{\rm ob}} > t_N (\nu _{\max }/\nu)) & = & {N_0 \over 4 \pi d_L^2}\!\! \int {r_0 \over r^2} e^{-{r_0 \over r}} dr\!\! \int {1 \over 2 \beta \Gamma ^2 u^2} du {k_B T^{\prime }_0 \over \Gamma u}\nonumber\\ && \times {\beta c \over u} \delta \left(r = {\beta c t^{{\rm ob}} \over u} \right) {h \Gamma u^2 \over k_B T^{\prime }_0} \delta\nonumber\\ &&\times \left(u = {k_B T^{\prime }_0 \over \Gamma h \nu } \right) \nonumber \\ & = & {h N_0 c \over 8 \pi d_L^2 \Gamma ^2} \int {r_0 \over r^2} e^{-{r_0 \over r}} dr \delta\nonumber\\ &&\times \left(r = {\beta c t^{{\rm ob}} \Gamma h \nu \over k_B T^{\prime }_0} \right) {\Gamma ^2 {(h \nu)}^2 \over {(k_B T^{\prime }_0)}^2} \nonumber \\ & = & {h N_0 c \over 8 \pi d_L^2 r_0} \left({r_0 \over \Gamma \beta c t^{{\rm ob}} }\right)^2 e^{-{r_0 k_B T^{\prime }_0 \over \Gamma \beta c t^{{\rm ob}} h \nu }}. \end{eqnarray} \tag{ 12 }$$

This equation can further be simplified by using the definition of t_N from Equation (5) and noting that at any given instance, the maximum observed temperature is given by T^ob_max = T'₀/Γu_min = T'₀/Γ(1 − β) ≃ 2ΓT'₀. Using these in Equation (12) leads to the final form

$$\begin{equation} F_{\nu } (t^{{\rm ob}} > t_N (\nu _{\max }/\nu)) = {h N_0 c \over 2 \pi d_L^2 r_0} \left({\Gamma t_N \over t^{ob} }\right)^2 e^{-{t_N \over t^{{\rm ob}}}{\nu _{\max } \over \nu }}, \end{equation} \tag{ 13 }$$

where ν_max = T^ob_max/h.

Equation (13) is the key finding of this paper. It implies that at late times, t^ob ≫ t_N, for a wide frequency range ν_max(t_N/t^ob) < ν < ν_max, the exponent is close to unity, and therefore a flat energy spectrum F_ν∝ν⁰ is expected below the thermal peak. This spectrum results from a simultaneous observation of thermal photons emitted from a large range of radii and angles to the line of sight. We emphasis again that this is a purely geometrical effect, as no additional radiative processes are considered. Furthermore, this aberration of the Planck spectrum results from the fact that the "classical" definition of the photosphere is extended here to include the entire space. We further note that Equation (13) provides, in addition, the decay law of the energy flux at late times, F_ν(t^ob ≫ t_N)∝t^ob−2.

While Equation (13) is derived under general conditions, we note that in practice its use in the context of GRBs is somewhat limited. First, the result is derived under the assumption that the inner engine already terminated. As long as the inner engine is active, the flux is dominated by photons emitted at angles θ < Γ⁻¹, and hence the spectrum is steeper than derived here.⁹ Second, the rapid decay of the flux, F_ν∝t^ob−2, implies that practically, the flux can be seen only up to t^ob ∼ 30 t_N, as at longer times, the flux may become too weak to be observed. Moreover, for a physical jet with opening angle θ_j, jet edge effects will become important once the contribution from photons at these angles become important. We further discuss these limitations in the context of GRBs in Section 5.

The analytical calculations presented above were checked with a numerical code. The code is a Monte Carlo simulation, based on earlier code developed for the study of photon propagation in relativistically expanding plasma (Pe'er & Waxman 2004; Pe'er et al. 2006b). This code is essentially identical to the code used for the numerical calculations that appear in Pe'er (2008), and a description of it appears there. We give here only a basic description of the code, for completeness, before presenting the numerical results and a comparison to the analytical approximation presented above.

The code is essentially a Monte Carlo simulation of Compton scattering between photons and electrons. The uniqueness of it lies in the fact that it calculates the interactions during a relativistic, three-dimensional expansion of the plasma. As a result, the probability of a photon to be scattered at any given instance, which is translated to the distance traveled by a photon between two consecutive scattering events, depends on the instantaneous radius and propagation direction of the photon. In every interaction, the full Klein–Nishina cross section is used in calculating the outgoing photon energy and propagation direction. Since a scattering event is calculated in the electron's rest frame, before and after every scattering, the photon 4-vector is being Lorentz transformed twice. First, into the (local) bulk motion rest frame of the flow which assumes an expansion at constant Lorentz factor Γ in the radial direction. A second Lorentz transformation is made into the electrons rest frame: in the bulk frame, the electrons assume a random velocity direction, with velocity drawn from a Maxwellian distribution with temperature T'_el(r)∝r^−2/3. The proportionality constant is determined using parameters characterizing the prompt emission in GRBs.

Initially, photons are injected into the plasma in a random position on the surface of a sphere at radius r_inj = R_d/(2πΓ²d). The depth d is taken as d = 20 in order to ensure that the probability of a photon to escape without being scattered is smaller than exp (− 20), i.e., negligible.¹⁰ The initial photon propagation direction is random, and its (local) comoving energy at the injection radius is equal to the plasma comoving temperature at that radius.

4.1. Numerical Results

We show in Figure 1 the positions of the last scattering events for N = 3 × 10⁶ photons in the r–θ plane. The three contour lines (thin black) are added to the plot in order to help demonstrating the probabilities of photons to be scattered at a given radius and angle. The solid (green) line is the photospheric radius, as calculated in Equation (3). Clearly, the photospheric radius provides only a first-order approximation to the last scattering events positions. It is obvious from the figure that photons decouple from the plasma at a range of radii and angles, necessitating the use of the probability density functions introduced in Section 3. We further added to the plot (dashed, blue lines), three equal arrival times contours (see Equation (5)), for three different values of the observed time. These contours demonstrate one of the key results: while at early times, t^ob ≲ few × t_N, the emission is dominated by photons emitted at angles smaller than Γ⁻¹, at later times, most of the contribution to the emission is from photons emitted from higher angles, θ > Γ⁻¹. This is the origin of the emerging power-law spectrum.

Zoom In Zoom Out Reset image size

In preparing the plot, we used parameters that can characterized GRBs, such as initial expansion radius r_i = 10⁸ cm, luminosity L = 10⁵² erg s⁻¹, and Lorentz factor Γ = 400. In accordance to Equation (3), the photospheric radius is angle independent at angles smaller than Γ⁻¹ = 2.5 × 10⁻³ rad, and is ∝θ² at higher angles. We point out that in Figure 1, the radius of the last scattering event is normalized to r₀ = R_d/2πΓ² (r₀ = 9 × 10¹⁰ cm for the parameters chosen here). Therefore, results obtained for arbitrary values of the free model parameters (L and Γ) that characterize astrophysical transients other than GRBs, such as AGNs or microquasars, are similar to the ones presented.

In Figure 2, we present the resulting observed spectrum, F_ν calculated at different times. The solid lines show the full numerical results at observed times t^ob/t_N = 10, 10², 10³, 10⁴. The thin dashed lines are the analytical approximation given by Equation (13). In calculating the maximum observed temperature, T^ob_max, we used the standard "fireball" model to determine the temperature evolution below the photosphere (see, e.g., Piran 2005; Mészáros 2006). Thus, the plasma assumes to accelerate between r_i and the saturation radius r_s = Γr_i, and continues with constant outflow velocity at larger radii. For the values of the parameters chosen, the plasma comoving temperature at the photospheric radius r₀ is T'(r₀) = (L/4πr²_ica)^1/4(r₀/r_s)^−2/3 = 0.54 L^−5/12₅₂r^1/6_{i, 8}(Γ/400)^5/3 keV, and the maximum expected temperature is T^ob_max = 2ΓT'(r₀) = 430 keV.

Zoom In Zoom Out Reset image size

The numerical results are indeed in very good agreement with the analytical approximation. We note that for the parameters chosen here, which can characterize many GRBs, t_N ≈ 10⁻⁵ s, and therefore the observed characteristic times of ≲second (in a single pulse) are translated into t^ob/t_N ⩽ 10⁴. At early times, the exponential decay occurs at (relatively) high energies, and thus a flat spectrum is not seen. However, at later times, the spectrum becomes flat over a wide energy range, consistent with the analytical approximations.

At intermediate times, there is a slight discrepancy between the numerical results and the analytical approximation, both at the high and low energies. This discrepancy demonstrates the limitation of the delta-function approximation used in Equation (13). At high energies, the numerical results show "bending" at early times, which is not captured by the analytical approximation. This results from the δ-function approximation in time used (Equation (5)). The underlying assumption used in deriving Equation (5) is that the observed delay time is solely determined by the radius and angle of the last scattering event. In reality, of course, as the photons diffuse below the photosphere, the time delay is determined by the full history of the photon propagation below the photosphere. At high energies, the spectrum is dominated by photons emitted with high Doppler shift, i.e., at small viewing angles; as at late times there are relatively very few such photons, dispersion in the last scattering position and angle leads to a high-energy decay.

At low energies, the discrepancy between the analytical prediction and the numerical results is explained as being due to a combination of two phenomena. First, the analytical approximation used does not consider the coupling between the probability density functions P(r) and P(θ). As is clear from Figure 1, for any given radius r, the observed viewing angle θ is limited from above. A second source of spread lies in the fact that the comoving photons energy distribution also have an internal spread (although, the average comoving energy is constant above few × r₀). As a result of these two effects, the low-energy spectrum is also slightly bent with respect to the simple power-law prediction of the analytical approximation. For a full treatment of these effects, one needs to obtain a full solution of the diffusion equation, which will be presented elsewhere.

In spite of being only a first-order approximation, it is clear from Figure 2 that the analytical predictions are in very good agreement with the numerical results of the observed spectrum. At intermediate times, the spectrum approaches a power-law distribution F_ν∝ν⁰ over a relatively wide energy band. Even at very early times, t ≳ t_N, it is clear that the observed spectrum deviates significantly from the classical Planck function.

In this manuscript, we have studied both analytically and numerically the spectrum resulting from a photospheric emission in relativistically expanding plasma. We showed that even for a purely thermal distribution of the comoving photon spectrum, aberration of the light results in an observed spectrum at late times that is very different from a Planck distribution. We showed that at late times, t ≫ t_N, the observed spectrum approaches a power law, with power-law index F_ν∝ν⁰ over a wide spectral range. This result is remarkably similar to the observed low-energy spectral index of many GRBs (Preece et al. 2000; Kaneko et al. 2006, 2008).

The origin of this non-trivial result lies in two facts: First, at late times, the emission is dominated by photons emitted off-axis (from angles θ ≫ Γ⁻¹, see Figure 1). Due to the strong dependence of the Doppler shift on the angle to the line of sight, these photons are seen at much lower energies than photons emitted on-axis. Second, at any given instance, an observer simultaneously sees photons that are emitted from a range of radii and angles (see Equation (5)). As each photon has a finite probability of being emitted from a given radius, and is seen at a particular Doppler shift, the observed spectrum and flux can only be described in terms of probability density functions (see Section 3). These functions provide a mathematical tool to describe the probability of photons to be emitted from radius r and into angle θ. Using these functions, we calculated an analytical approximation to the observed spectrum (Equation (13)), which is the main result of this work. The analytical approximation was tested with a Monte Carlo simulation that tracks the evolution of thermal photons in relativistically expanding plasma (Section 4). The simplified analytical calculations are found to be in very good agreement with the accurate numerical results (see Figure 2).

In spite of the success in reproducing the low-energy spectral slopes at late times, the theory is still not complete. First, we note that the theory gives the spectrum following an abrupt termination of the inner engine. As long as the inner engine is active, the spectrum is dominated by photons emitted at angles θ ≲ Γ⁻¹, and the modification to the Planck spectrum is not as significant (note that nonetheless, some modification must exist due to contribution of photons from higher angles. It is easy to show¹¹ that in this case the Rayleigh–Jeans part of the spectrum is F_ν∝ν). Therefore, comparing our result with the data must be done with caution, as it should be relevant only during the decay phase of a GRB pulse. Moreover, we consider here a simplified scenario in which the Lorentz factor is constant. In practice, of course, the flow fluctuates, which implies that (1) the parameters of the outflow may vary in time, and (2) there may exist sub-photospheric dissipation, which will modify the properties of the photospheric emission.

The most noticeable drawback of the theory as stated here is the very short characteristic time scale. For parameters characterizing GRBs, the thermal peak (T^ob_max) naturally falls at the sub-MeV range (see Section 4.1), hence the observed peak can naturally be explained as having a thermal origin. However, using the same parameters leads to very short characteristic time scale, t_N ≈ 10⁻⁵ s. Since above t_N the flux decays rapidly, F_ν(t)∝t⁻² (see Equation (13)), one expects a relatively weak thermal signal at t ≫ t_N, which can still be very short time.

There are several possible solutions to this problem, although, admittedly none is yet complete and further research is being done. One possibility is that in many GRBs the Lorentz factor is not few hundreds, but only few tens. While this cannot be the case for the Large Area Telescope (LAT) GRBs, we note that these are only very small fraction, ∼5% of the total GRBs seen by Fermi. Without clear detection of >MeV emission, there is no certainty that the bulk Lorentz factor is indeed few hundreds, as is typically assumed. Another simplified assumption is that the calculation of t_N in Equation (5) is based on the assumption of constant outflow velocity. This assumption is too simplified: in order to reach high Lorentz factor, the plasma needs to undergo an acceleration phase, and so the average Lorentz factor below the photosphere is less than the terminal Lorentz factor. As a result, we expect that in practice the characteristic time scale relevant for thermal emission in GRBs is longer than the one considered in Equation (5). Other simplified assumptions used are the assumption of spherical explosion (as opposed to the more realistic jet-like explosion, in which Γ = Γ(θ)), and the neglection of a second photosphere, produced by a large number of pairs that may exist in the plasma. Relaxing these assumptions would also change the characteristic time scale of photon escape. We expect these effects to lead to characteristic time scale that is longer than that considered in Equation (5). However, a full treatment of these effects is beyond the scope of this paper. Here, for the first time we introduce the concept of geometrical broadening of the photospheric emission, which is done in the most simplified scenario.

The exact delay time depends on several uncertain conditions. One is the content of the fireball: for example, a Poynting-flux-dominated fireball is expected to have slower acceleration than matter-dominated fireball (Drenkhahn 2002; Drenkhahn & Spruit 2002; Giannios & Spruit 2005, 2006). Another is the baryon load (Ioka 2010), and in particular the baryon distribution along the jet (Morsony et al. 2007; Lazzati et al. 2009; Mizuta et al. 2010): although in this manuscript we considered a steady outflow, clearly the outflow in GRBs is characterized by regions of higher and lower densities, and is thus not steady (edge effects due to the finite opening angle may also play a role at late times).

An additional source of discrepancy between the theoretical predictions developed in this paper and the observed spectrum, lies in the fact that the theory here does not consider any additional, non-thermal radiative processes. As shown here, photospheric emission is capable of reproducing the peak energy and the low-energy spectral slope (α in the "Band" function) seen in GRBs at late times. However, photospheric emission is not capable of producing high-energy photons (above T^ob_max ≲ MeV), as are seen in some GRBs by the LAT detector on board the Fermi satellite. The inclusion of high-energy, non-thermal photons, necessitates additional radiative mechanisms that must take place following dissipation processes that occur above the photosphere (e.g., in GRB080916C analysis of the high-energy emission imply Poynting-dominated outflow; see Zhang & Pe'er 2009). Additional radiative mechanisms, such as synchrotron emission or Compton scattering, naturally produce a broadband energy spectrum, and thus may contribute not only to the high-energy spectrum but to the low-energy part (below the thermal peak) as well. The overall observed spectra below the thermal peak is thus generally expected to be hybrid—i.e., composed of both thermal and non-thermal parts. The exact contribution of the non-thermal part may vary from burst to burst, depending on the values of the free model parameters, such as the radius of the photosphere, the dissipation radius, the strength of the magnetic field, etc.

In addition, the inclusion of a non-thermal part resulting from dissipation above the photosphere leads to masking the effect: the non-thermal photons, produced above the photosphere are naturally delayed with respect to the photospheric photons. Thus, termination of the inner engine, which leads to gradual change in the photospheric spectrum, is expected to affect the non-thermal spectrum at later times. Hence, during the decay phase of the photospheric component, the non-thermal part of the spectrum is expected to be more pronounced, and the signal may be smeared. We stress, though, that the theory does not give predictions for the dissipation radius that creates the non-thermal photons, hence it is impossible to calculate the exact effect without additional assumptions.

In principle, the inclusion of thermal photons contributes to the high-energy, non-thermal part of the spectrum as well, as these photons serve as seed photons for Compton scattering by the non-thermal electrons and hot pairs (Rees & Mészáros 2005; Pe'er et al. 2005, 2006a; Lazzati & Begelman 2010; Beloborodov 2010). The exact contribution of the thermal photons depends on the optical depth at the dissipation radius. If the dissipation radius is close to the photosphere, the resulting spectrum has a complex shape that is very different than either a thermal spectrum or the optically thin synchrotron—synchrotron self-Compton (SSC) model predictions (Pe'er et al. 2005, 2006a). On the other hand, if the dissipation occurs at large radii, the two components, thermal and non-thermal, can be decoupled, a fact that can lead to a clear identification of the thermal component, as is in the case of GRB090902B (Ryde et al. 2010; Pe'er et al. 2011).

One clear prediction of the results presented here is that if a thermal component contributed significantly to the observed spectrum, then the low-energy spectral slope (α) is expected to vary with time: at early times, when the inner engine is active, α is expected to be steep (F_ν∝ν¹), while at later times, F_ν∝ν⁰. This result is expected to be correlated with a decay of the peak energy and the thermal flux, and was possibly observed in several bursts (e.g., Crider et al. 1997; Ghirlanda et al. 2007; Ryde & Pe'er 2009; Page et al. 2009). We note though that in many GRBs, the data are not easily compared to the theoretical prediction, since the outflows in GRBs are generally not smooth, but fluctuate significantly. Therefore, photons originating from different emission episodes (each characterized by its own mass-loss ejection rate and its own Lorentz factor) are often superimposed. Thus, in fact, thermal photons may only be identified using detailed time resolved spectral analysis of separate flares. This task was carried by Ryde & Pe'er (2009), which were indeed able to identify the thermal component in a large sample of bursts. Once done, this analysis method can be further used to deduce the properties of the GRB outflow (Pe'er et al. 2007).

In addition to the prompt emission phase in GRBs, thermal activity may occur as part of the flaring activity observed in the early afterglow phase of many GRBs (Burrows et al. 2005; Falcone et al. 2007; Chincarini et al. 2010; Margutti et al. 2010). The exact nature of these flares is currently not yet clear. As it is plausible that the flares result from renewed emission from the inner core, a renewed thermal emission may occur. The analysis presented here may therefore apply in the study of the late time flares as well.

While analysis of a GRB prompt emission spectrum serves as a major motivation to this work, we note that the analysis presented here is general. With modified parameters values, the analysis applies to emission from any astronomical transient, which is characterized by relativistic outflow speeds and a high density core. Thus, we expect the analysis carried here to be valid for objects such as AGNs and microquasars. This is valid because the exact nature of the radiative process that produces the thermal photons is of no importance, as long as the emission occurs deep enough in the flow, where the optical depth is high, τ ≫ 1.

We would like to thank Mario Livio, Željka Bošnjak, Chryssa Kouveliotou, Charles Dermer, Ralph A. M. J. Wijers, Peter Mészáros, Martin Rees, and Bing Zhang for many useful discussions and comments. A.P. is supported by the Riccardo Giacconi Fellowship award of the Space Telescope Science Institute.