{'fraction_non_alphanumeric': 0.03963999900276732, 'fraction_numerical': 0.023397571738425868, 'mean_word_length': 4.1003242418462715, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Photospheric emission from stratified jets

We explore photospheric emissions from stratified two-component jets, wherein a highly relativistic spine outflow is surrounded by a wider and less relativistic sheath outflow. Thermal photons are injected in regions of high optical depth and propagated until they escape at the photosphere. Due to the presence of shear in velocity (Lorentz factor) at the boundary of the spine and sheath region, a fraction of the injected photons are accelerated via a Fermi-like acceleration mechanism such that a high energy power-law tail is formed in the resultant spectrum. We show, in particular, that if a velocity shear with a considerable variance in the bulk Lorentz factor is present, the high energy part of observed Gamma-ray Bursts (GRBs) photon spectrum can be explained by this photon acceleration mechanism. We also show that the accelerated photons may also account for the origin of the extra hard power-law component above the bump of the thermal-like peak seen in some peculiar bursts (e.g., GRB 090510, 090902B, 090926A). It is demonstrated that time-integrated spectra can also reproduce the low energy spectrum of GRBs consistently due to a multi-temperature effect when time evolution of the outflow is considered. Finally, we show that the empirical Ep-Lp relation can be explained by differences in the outflow properties of individual sources.


INTRODUCTION
Gamma-ray Bursts (GRBs) are the most powerful explosions in the Universe. Prompt emission of GRBs are mainly observed in the energy range of 10 keV -MeV and show rapid time variability in their lightcurves. Their spectra are often modeled by an empirical 'Band' function (Band et al. 1993;Preece et al. 2000;Kaneko et al. 2006Kaneko et al. , 2008, which is a smoothly jointed broken powerlaw, whose physical origin is not yet identified. The typical low and high energy photon indices are distributed around α ph ∼ −1 and β ph ∼ −2.5, respectively, while the spectral peak (break) energy is distributed around E p ∼ a few 100 keV.
The most widely discussed model for the prompt emission mechanism is the internal shock model (Rees & Meszaros 1994;). In this model, it is assumed that relativistically moving shells emanate from the central engine with diverse velocities, and shocks form due to collisions of the shells. Shocks are accompanied by particle acceleration, and eventually gamma-rays are produced by relativistic electrons via synchrotron radiation. The internal shock model can naturally explain the observed rapid time variability in the lightcurve and non-thermal nature of the spectra. However, it is known that the model suffers from poor radiation efficiency, since only the kinetic energy associated with the relative motion of shells can be released (Kobayashi et al. 1997;Lazzati et al. 1999;Guetta et al. 2001;Kino et al. 2004). This contradicts with observations which show high efficiencies of a few tens of percent (Fan & Piran 2006;Zhang et al. 2007). Another difficulty is the low energy spectral index (α ph ). A nonnegligible fraction of observed GRBs show hard spectra at low energies which cannot be explained by synchrotron emission (Crider et al. 1997;Preece et al. 1998;Ghirlanda et al. 2003).
In this model, prompt gamma-rays are released at the photosphere when the fireball becomes optically thin. The peak energy of the spectra is determined by the temperature at the photosphere. This model has an advantage in that high emission efficiencies can be naturally achieved. Strong observational support for this scenario has been provided by the recent detection of quasi-thermal emission by the Fermi Large Area Telescope (LAT), notably GRB 090902B (Abdo et al. 2009;Ryde et al. 2010Ryde et al. , 2011Pe'er et al. 2012).
However, while there are many advantages, the pho-tospheric emission model must overcome several difficulties. The major flaw is reproduction of the observed broad non-thermal spectra, since the photons that originate from regions of very high optical depths are inevitably thermalized. Photospheric emission can, in principle, explain the hard low-energy spectra that cannot be accounted for by synchrotron emission, by considering a superposition of emission components which have different temperatures (multi-color Blackbody). This is because the Rayleigh-Jeans part of the emission component has a much harder photon index (α ph = 1) than those inferred in all observed GRBs. It is noted, however, that the existence of outflow models that can naturally account for the low energy spectra is still debated. Recent theoretical studies which explore peak energies of photospheric emissions based on spherical (one-dimensional) outflow models (e.g., Levinson 2012;Beloborodov 2013;Vurm et al. 2013) suggest that it is difficult to regulate the temperature of the emissions to be well below the typical observed peak energies E p , particularly when dissipation is present (but see, e.g., Levinson & Globus 2013, for a mechanism that may overcome this difficulty). Although dissipative processes are not considered, an outflow structure in twodimensions suggested by (Lundman et al. 2013) may be a key ingredient to resolve this difficulty. Based on hydrodynamical simulations of jet propagation (Zhang et al. 2003), Lundman et al. (2013) modeled an outflow with a Lorentz factor gradient in the lateral direction. By solving the transfer of photons within the jet, they indeed showed that the low energy spectra can be reproduced by the multi-color temperature effects due to the superposition of the photons released from different lateral positions that have various Lorentz factors. On the other hand, the high energy part of the observed spectra is difficult to reproduce by the superposition of thermal emissions, since an unreasonably high temperature must be realized within the outflow. To overcome this difficulty, dissipative processes around the photosphere have been suggested in previous studies. Various mechanisms such as shocks (Pe'er et al. 2005(Pe'er et al. , 2006Ioka et al. 2007;Lazzati & Begelman 2010), magnetic reconnection (Giannios 2006;Giannios & Spruit 2007;Giannios 2008) and proton-neutron nuclear collisions (Beloborodov 2010;Vurm et al. 2011) have been proposed. All models are accompanied by the generation of relativistic electrons (and positrons) 9 that upscatter the thermal photons to create the non-thermal spectra (Giannios 2012). However, it is quite uncertain whether dissipative processes can operate efficiently enough to deposit in relativistic electrons (and positrons) the copious amounts of energy needed to reproduce observations. Although the details are not well studied, Fermi-like acceleration of photons at the shocks formed below the photosphere due to bulk Compton scattering can be regarded as an alternative mechanism for producing high energy non-thermal spectra (Eichler 1994). Presence of relativistic electrons (and positrons) is not required in this scenario, since the energy of the background fluid can be directly transferred to the photons as they propagate. Shocks formed below the photosphere are inevitably mediated by radiation (Levinson & Bromberg 2008). The properties of relativistic radiation mediated shocks have been explored in various astrophysical contexts (e.g., shock breakouts) including GRBs. The most detailed and fully self-consistent analysis of relativistic radiation mediated shocks have been performed by Budnik et al. (2010) (see also Katz et al. 2010). Although not as sophisticated as the above study, Bromberg et al. (2011) explored the properties of photon acceleration in the context of shocks within GRB jets. These studies have shown that the existence of accelerated photons above the thermal peak energy is an inherent feature in these shocks. However, all the above studies assume a onedimensional stationary planar shock in their analysis, and the obtained spectra are far from those typically observed. Therefore, for a more accurate estimation, a global photon transfer calculation incorporating the multi-dimensional structure of the jet must be performed. Ioka et al. (2011) explored the properties of shocks that are formed in radiation dominated jets. They claimed that photons which are accelerated at the shock can generate a spectra which are close to the 'Band' function. However, since their argument is based on a simple analytical framework in which the transfer of photons is not solved, it is obvious that detailed calculations are required to form firmer conclusions.
While the previous studies focused on the shocks, any form of large velocity gradients that are present in the outflow can give rise to photon acceleration. For instance, a Fermi-like mechanism can also operate when velocity shear exists in the transverse direction of the outflow, just as in the case of the acceleration of the cosmic rays (e.g., Jokipii et al. 1989;Ostrowski 1990Ostrowski , 1998. Indeed, GRB jets are likely to posses rich internal velocity structures in the transverse directions. For example, numerical studies which explore the launching mechanism of the relativistic jets suggest that the produced jets have transversely stratified structures in which the bulk Lorentz factor increases towards the jet axis (McKinney 2006;Nagataki et al. 2007;Nagataki 2009;McKinney & Blandford 2009;Nagataki 2011). Since the highly magnetized jets seen in these simulations are cold, subsequent magnetic field dissipation must occur in order to produce a hot fireball which leads to an efficient photospheric emission. Although the flow structure after the dissipative process is highly uncertain, the stratified structure that originated at the central engine may remain even at this stage. Moreover, rich velocity structure can also develop after the launching phase as the jets drill through the progenitor star envelopes. A large number of hydrodynamic simulations regarding the jet propagation show that large velocity gradients appear within the jet (in radial as well as transverse directions) due to the interaction with the stellar envelope (Zhang et al. 2003;Mizuta et al. 2006;Morsony et al. 2007;Lazzati et al. 2009;Mizuta et al. 2011;Nagakura et al. 2011). Notably, large velocity gradient in the transverse directions produced by the strong recollimation shocks seen in these studies will likely provide an efficient acceleration site for photons as pointed out by Levinson (2012).
Motivated by this background, we explore the photon acceleration within a jet with velocity structures in the transverse direction and its effect on the resulting spectra of the photospheric emissions by solving the transfer of photons. As mentioned earlier, mainly focusing on the origin of the low energy spectra, Lundman et al. (2013) carried out a similar study. By performing a photon transfer calculation, they evaluated the photospheric emission from a relativistic jet that has a continuously decaying velocity profile in lateral direction Γ ∝ θ −p at the outer region. However, while they succeeded in reproducing the spectra below the peak energy, the imposed velocity gradient was not large enough to trigger efficient photon acceleration and, therefore, the high energy part was too soft compared to the Band spectra. With the aim to investigate the cases with efficient photon accelerations, we focus on jets that have a velocity shear in the transverse direction. This can be considered as a limiting case of p = ∞ in Lundman et al. (2013) work, a parameter region which was not explored. As a first step, this study considers a simple stratified two-component jet structure, wherein a highly relativistic spine outflow is surrounded by a wider and less relativistic sheath outflow. We show, in particular, that the high-energy part of the Band spectra (photon index β ph ) can indeed be reproduced by photons accelerated at the shear layer. In addition, we also demonstrate that the low-energy part of the spectra (photon index α ph ) can be reproduced consistently within the framework of the present study when time evolution of the outflow is considered.
The paper is organized as follows. In §2, we describe our model and numerical procedures. We present our main results in §3. Discussions of the implications of our results are given in §4. The summary of our main findings are given in §5.

MODEL AND METHODS
In the present study, we evaluate the photospheric emissions from an ultra-relativistic outflow which has a spine-sheath jet structure. The spine is defined as a region of conical outflow with a half-opening angle θ 0 , while the sheath is a region of outflow which surrounds the spine and extends up to an angle of θ 1 (> θ 0 ) (see Fig. 1). Each region has different fluid properties and a steady radial outflow is assumed. Under the assumption that the energy density in the magnetic field and the dissipated kinetic energy are subdominant, the fluid properties of the spine and sheath regions are described by the standard adiabatic "fireball" model (e.g., Piran 2004;Mészáros 2006). Here we give a brief review of the fireball model and show how the radial evolution of the Lorentz factor (velocity), Γ, and electron number density, n e , which are necessary background fluid information for evaluating the photon transfer, is determined.
The dynamics of the fireball can be characterized by three independent parameters, which are the initial fireball radius, r i , the kinetic luminosity, L, and the dimensionless entropy which characterizes the baryon loading, η ≡ L/Ṁ c 2 , whereṀ and c are the mass outflow rate ri rinj, rs1 θ1 rs0 spine sheath rph0 θ0 rph1 Fig. 1.-Schematic picture of the employed model. A fast spine jet is embedded in a slower sheath outflow. Spine and sheath start to accelerate at a radius r i , and the acceleration continues up to r s0 and r s1 in the spine and sheath region, respectively. Since the dimensionless entropy of the spine η 0 is larger than that of the sheath η 1 , the saturation radius and the terminal Lorentz factor of the spine (r s0 = η 0 r i and Γ 0 = η 0 ) are larger than that of the sheath (r s1 = η 1 r i and Γ 1 = η 1 ). The photospheric radius of the spine r ph0 is smaller than that of the sheath r s1 . In our model, thermal photons are injected at a the saturation radius of the sheath r inj = r s1 , and the transfer is solved up to the radius in which the optical depth is much lower than unity. Photons which cross the boundary layer of the spine and sheath multiple times gain energy and produce a non-thermal tail in the resultant spectrum. and the speed of light, respectively. Initially, the fireball accelerates as it expands by converting its internal energy into kinetic energy (acceleration phase). As a result, the bulk Lorentz factor of the flow increases with radius, r, as Γ(r) ≃ (r/r i ). If η is larger than the critical value given by η c = (σ T L/8πr i m p c 3 ) 1/4 (Piran 2004;Mészáros 2006), the fireball becomes transparent to radiation during the acceleration phase, where σ T and m p are the Thomson cross section and proton rest mass, respectively. On the other hand, if η < η c , the fireball becomes transparent after the acceleration has ceased and the flow asymptotically approaches a constant Lorentz factor Γ = η (coasting phase). In this case, the acceleration continues up to the saturation radius, r s = ηr i . In the present study, we focus on the latter case and consider the parameter space in which η < η c is satisfied in both the spine and sheath regions. Hence, the radial evolution of Lorentz factor is given by Γ(r) = r rs for r ≤ r s , η for r > r s . (1) In the relativistically expanding outflow, the electron number density in the comoving frame is given by where β is the velocity of the flow normalized by the speed of light. Here, we assumed that there are no strong dissipative processes in the outflow which may create copious pair plasma. Since we consider the case of η < η c , the electron number density decreases with radius as n e ∝ r −3 below the saturation radius (r ≤ r s ) and as n e ∝ r −2 at larger radii (r > r s ). Given the electron number density and bulk Lorentz factor of the flow, the optical depth to Thomson scatterings for the photons propagating in the radial direction to reach infinity can be evaluated as where we have assumed Γ ≫ 1. Here, r ph is the photospheric radius which corresponds to the radius where the optical depth becomes unity (τ = 1).

Spine-sheath structure
The radial profiles of the fluid quantities (n e (r) and Γ(r)) in the spine and sheath regions are different, since we impose different values on their fireball parameters (r i , L and η). In the present study, while we assume the same value for r i in both regions, the dimensionless entropy of the spine, η 0 , is taken to be larger than that of the sheath, η 1 . Hence, the saturation radius in the spine region, r s0 = η 0 r i , is larger than that of the sheath, r s1 = η 1 r i , and the terminal Lorentz factor of the former (η 0 ) is larger than that of the latter (η 1 ). Regarding the evolution of the Lorentz factor in the spine and sheath, they have equal values (Γ 0 (r) = Γ 1 (r) = r/r i ) up to r = r s1 . At larger radii (r > r s1 ), velocity shear begins to develop and the difference in the Lorentz factor increases with radius up to r s0 , since the spine is in the acceleration phase while the sheath is in the coasting phase (Γ 1 = η 1 ). Thereafter, the spine also enters the coasting phase (Γ 0 = η 0 ), and the difference in the Lorentz factor is constant. In determining the kinetic luminosities of the spine, L 0 , and sheath, L 1 , we assume that the mass outflow rate is equal in both regions (L 0 /η 0 = L 1 /η 1 ). Therefore, kinetic luminosity of the spine is larger by a factor of η 0 /η 1 . Under these assumptions, the photospheric radius in the sheath, r ph1 , is larger than that in the spine, r ph0 , by a factor of (η 0 /η 1 ) 2 (see equation (4)). A schematic picture of the employed model is given in Fig. 1.
Hereafter, the quantities corresponding to the spine and sheath regions are denoted by subscript 0 and 1, respectively.

Photon Transfer in a Spine-Sheath Jet
Having determined the background fluid properties (Γ and n e ), we evaluate the resultant photospheric emission by solving the propagation of photons which are injected far below the photosphere. The photon transfer is evaluated by performing a three-dimensional test particle Monte-Carlo simulation. In GRB jets, opacity of photons is strongly dominated by the scatterings with electrons. Therefore, we neglect the absorption process and only consider the scattering process by the electrons in our calculations. Furthermore we do not take into account the thermal motion of the electrons in evaluating the scattering for simplicity.

Initial condition
Initially, the photons are injected within the jet at the surface of a fixed radius where the velocity shear begins to develop r inj = r s1 . For the cases considered in this study, r inj is always located far below the photosphere (τ (r inj ) ≫ 1). Therefore, a tight coupling between the photons and matter is expected. For this reason, we can safely assume that the photons have an isotropic distribution with energy distribution given by a Planck distribution in the comoving frame. According to the fireball model, the radial evolution of the comoving temperature is given by where a is the radiation constant. Hence, we adopt the temperature at the corresponding radius given by above equation T ′ inj = T ′ (r inj ) for the comoving temperature of the injected photons. While the photons are isotropic in the comoving frame, they are strongly beamed in the laboratory frame due to the Doppler boosting effect. Due to this effect, the radiation intensity of the blackbody emission in the laboratory frame is given by where Γ inj = Γ(r inj ) is the bulk Lorentz factor of the flow at r = r inj determined from equation (1). Here, B ν (T ′ , ν) = 2hν 3 c −2 [exp(hν)/k B T ′ − 1] −1 is the Planck function, where h and k B are the Planck constant and the Boltzmann constant, respectively, and D(Γ, θ v ) = [Γ(1 − βcosθ v )] −1 is the Doppler factor, where θ v is the angle between the photon propagation direction and the fluid velocity direction (radial direction). In our calculations, the initial propagation direction and frequency of the injected photons are drawn from a source of photons given by the above equation.
It is emphasized that the results of our calculation are insensitive to the assumed position of the injection radius as long as r inj ≤ r s1 is satisfied. This is because, at a radius far below the photosphere (τ (r) ≫ 1), the photon energy distribution evaluated by solving the photon transfer does not deviate from the Planck distribution if velocity shear is not present (see next section), and its temperature evolution is well described by equation (5). The temperature of the injected photons in the spine, T ′ inj0 , is higher than that in the sheath, T ′ inj1 , by a factor of (η 0 /η 1 ) 1/2 , since the kinetic luminosity of the former is higher by a factor of (η 0 /η 1 ) 2 . Correspondingly, the luminosity of the injected photons in the laboratory frame (L inj ∝ r 2 inj Γ 2 inj T ′ 4 inj ) in the spine is higher than that of the sheath by a factor ∼ (T ′ inj0 /T ′ 4 inj1 ) ∼ (η 0 /η 1 ) 2 .

Boundary conditions
After the photons are injected, we track their path within the jet in three-dimensions by performing Monte-Carlo simulations ( §2.2.3) until they reach the outer or inner boundary of the calculation. Using spherical coordinates (r,θ,φ), the outer boundary r is set at a radius r out = 500r ph0 where the photons can be safely considered to have escaped since the optical depth is τ (r out ) = 2 × 10 −3 ≪ 1. While there is no boundary in the φ direction, outer boundary in the θ direction is set at θ out = θ 1 which corresponds to the edge of the whole jet. As for the inner boundary, we adopt a radius slightly below the injection radius r in = 0.5r inj . For photons which have reached the outer boundaries, we assume that they escape freely to r = ∞ without being scattered or absorbed. On the other hand, we assume that the photons are simply absorbed in the inner boundary. It is noted, however, that the fraction of absorbed photons is negligible, since most of the photons in ultra-relativistic outflows are strongly collimated due to the relativistic beaming effect and essentially streamed outward (e.g., Pe'er Beloborodov 2011).
The spectra of the emission are evaluated from the photons which have reached the outer boundaries. Due to the relativistic beaming effect, these photons are highly anisotropic and mostly concentrated within a cone of half-opening angle ∼ Γ −1 in the direction of the fluid velocity (radial direction) at the last scattering position. Hence, the observed emission spectra depend significantly on the angle between the direction to the observer and the jet axis, θ obs . In the present study, we evaluate the spectrum for observers located at direction θ obs by recording all photons which have reached the outer boundary that are propagating in direction within a cone of half-opening angle Γ −1 , which is small enough to regard that the emission is uniform within the cone. From the recorded photon flux, we calculate the isotropic equivalent luminosity by multiplying the photon flux by a factor 4π/dΩ, where dΩ = 2π[1 − cos(1/Γ)] is the solid angle of the cone.

Monte-Carlo simulation for solving photon transfer
Here we briefly describe the Monte-Carlo code used to solve the photon transfer. As noted previously, we neglect the thermal motions of the electrons and only take into account the scattering processes. Hence, the rest frame of the fluid is equivalent to that of the electrons. Under the above assumptions, the propagation of the photons is performed by directly tracking the path of the individual photons in the three-dimensional space of the calculation. Each photon has a specified position propagation direction and frequency, and these quantities are updated by using a uniform random number.
Within the jet, the photons travel along straight paths before they are scattered by the electrons. Firstly, the code determines the distance for the photons to travel before the scattering by drawing the corresponding optical depth δτ . The probability for the selected optical depth to be in the range of [δτ , δτ + dτ ] is given as exp(−δτ )dτ . Then, from the given optical depth δτ , the distance l to the scattering event is determined from the integration along the straight path of photons which can be expressed as where θ v is the angle between the direction of fluid velocity and photon. Here, σ sc is the total cross section for the electron scattering and is given as in our code, where σ KN is the total cross section for Compton scattering, and ν cmf is the frequency of the photon in electron (fluid) comoving frame. (The frequency ν cmf is evaluated by performing a Lorentz transformation using local fluid velocity). Given the distance l from the above equation, we update the position of the photons to the scattering position by shifting them from the initial position with the given distance in the initial direction of photon propagation. Note that, unlike the case of equation (3), the optical depth calculated by equation (7) is not limited to photons propagating in the radial direction. The path of integration is along the straight path of photons which can be in an arbitrary direction. For a given value of δτ , the distance l strongly depends on the propagation direction of the photons in the case of a relativistic flow (Γ ≫ 1). As is obvious from the above equation, the mean free path of photons Hence, a photon tends to travel a larger distance in the fluid velocity (radial) direction since the mean free path of the photon tends to be larger. Hereafter, quantities measured in the comoving frame of the fluid (electron) are denoted by tilde.
In evaluating the integration in equations (7), we employ two different methods depending on the frequency and position of the photon. For photons located above the saturation radius r ≥ r s (Γ = const) that satisfy hν cmf ≤ 100 keV (σ sc = σ T = const), analytical integration can be performed as shown by Pe'er (2008). Consider a photon path originating from a radius r I that has an angle θ v,I with respect to the fluid velocity (radial) direction at the original position. In this case, the optical depth to reach a radius r II along the straight photon path can be expressed as where θ v,II is the photon angle at the final position (r = r II ). In the second equality, we have used the relation r I sinθ v,I = r II sinθ v,II that holds for an arbitrary straight line. Hence, in this case, we determine the corresponding propagation length ∆l I−II = r II cosθ v,II − r I cosθ v,I from the drawn optical depth by solving equations (9). (Note that equation (9) is solved separately in the spine and sheath region, since the values of r ph and Γ (β) are different in each region).
On the other hand, when the photons have higher energies (hν cmf > 100 keV) or are located below the saturation radius (r < r s ), the integration is solved numerically. In this case, we divide the calculation region into a mesh in spherical coordinates (r,θ,φ). Then the integration is done by assuming that the physical quantities (velocity and number density) within the individual mesh are uniform and have the values corresponding to those at the position of the mesh center. We adopt 500 grid points which are logarithmically spaced for the mesh in the r coordinate between the inner boundary r in and outer boundary r out . As for the mesh in the θ coordinate, we adopt 800 uniformly spaced grid points in the range θ ≤ θ 1 . 1600 uniformly spaced grid points are adopted for the φ coordinate (0 ≤ φ < 2π). It is noted that the resolution of the grid is sufficiently high to reproduce the result obtained by the analytical solution given by equation (9) (corresponding to infinite resolution) in cases when r ≥ r s and hν cmf ≤ 100 keV are satisfied.
Given the position for the scattering from the above procedure, the four-momentum (the energy and propagation direction) of a photon after the scattering is determined based on the differential cross section for Thomson and Compton scattering. In our code, the scattering process is evaluated in the rest frame of the fluid (electron). First, the four-momentum of the photon before the scattering is Lorentz transformed into the fluid rest frame. For the photons that satisfy hν cmf ≤ 100 keV, the differential Thompson cross section is used, while the differential Compton cross section is used at higher energies (hν cmf > 100 keV). The scattering angle or equivalently the propagation direction of the outgoing photon in the fluid rest frame is drawn from the differential cross sections. Regarding the energy of the outgoing photons, we assume that it is conserved before and after the scattering (elastic scattering) in the case of Thomson scattering (hν cmf ≤ 100 keV). On the other hand, in the case of Compton scattering (hν cmf > 100 keV), energy loss due to the recoil effect is properly taken into account. The code then Lorentz transforms the outgoing photon fourmomentum back into the laboratory frame.
The above procedure is repeated until all the injected photons reach the boundary of the simulation grids.

RESULTS
In this section, we show the obtained photon spectra based on the model described in the previous section. We inject N = 2 × 10 8 photon packets in each calculation. In all cases, we employ a fixed value of θ 1 = 1 • for the half opening-angle of the jet, which is smaller than that of the typically observed values. It is emphasized, however, that the resulting spectra do not vary for wider jets (larger θ 1 ) as long as the observer angle stays in the range θ obs θ 1 − Γ −1 . This is simply because the emission from regions located at a angle |θ − θ obs | Γ −1 is negligible due to the relativistic beaming effect. Therefore, merely to reduce computational cost, we adopt the relatively small value in this study. We use r i = 10 8 cm for the initial radius of the fireball in all cases.

Uniform (Non-Stratified) Jet
Before we look at a stratified jet, we first present results for a one-component uniform jet that does not have structures in the θ direction (θ 0 = θ 1 = 1 • ). In this case, a thermal spectrum is expected, in contrast to a stratified jet as we will discuss later. The isotropic equivalent kinetic luminosity and the dimensionless entropy (terminal Lorentz factor) are set to be L 0 = 10 53 erg/s and η 0 = 400, respectively. As described in the previous section, we inject the photons at a radius of r inj = 4 × 10 10 η 0,400 r i,8 cm with intensity given by a blackbody for detail), where η 0,400 = η 0 /400, L 0,53 = L 0 /10 53 erg/s and r i,8 = r i /10 8 cm. The corresponding optical depth at the injection radius is τ (r inj ) ∼ 23. The results are insensitive to the value of r inj as long as τ ≫ 1 is satisfied as noted in §2. The advected photons lose energy adiabatically due to the expansion of the flow until the coupling with matter becomes weak near the photosphere (r ph0 ≃ 9.2 × 10 11 L 0,53 η −3 0,400 cm). The expected temperature at the photosphere can be calculated as k B T ′ ph0 = 0.38r keV. The observed peak energy of the photospheric emission is expected to be keV, since the photon energy is boosted by a factor D ∼ 2η 0 due to the Doppler effect. Also the luminosity of the emission can be estimated as Fig. 2 with a red solid line. Here we assume that the observer is aligned to the jet axis (θ obs = 0 • ). It is noted, however, that the result does not change if the observer angle stays in the range θ obs θ 1 − η −1 0 ∼ 0.86 • , since the fluid properties within the beaming cone (|θ − θ obs | 0.14 • η −1 0,400 ) do not change. In the figure, we also show an analytical solution for the expected emission when the photons are in complete thermal equilibrium up to a radius r which can be obtained as where θ is angle between the line of sight and fluid velocity (radial) direction. The green, blue and purple solid lines correspond to the solutions for r = r ph0 (τ = 1), r = r ph0 /2(τ = 2) and r = r ph0 /4(τ = 4), respectively. Regarding the case of r = r ph0 , the peak energy, E p , and luminosity, L p ∼ ν p L ν (ν p ), where ν p = E p /h, of the obtained spectrum agree well with the rough estimate given earlier. As for the cases of r = r ph0 /2 and r = r ph0 /4, both the peak energy and the luminosity are larger by a factor ∼ 2 2/3 ∼ 1.6 and ∼ 4 2/3 ∼ 2.5 than the case of r = r ph , respectively, since the temperatures at these radii are larger by the same factor. As shown in the figure, the peak luminosity of the numerical result is in good agreement with that of the analytical estimate for r = r ph . On the other hand, the spectrum extends up to higher energies and the peak energy is close to that for r = r ph /2. This is due to the fact that the coupling between the photon and matter is not complete near the photosphere τ ∼ 1 as shown in the previous studies (Pe 'er 2008;Beloborodov 2011;Bégué et al. 2013). As a result, photons which decouple with the matter at moderate optical depth (τ 5) are observed at higher energies.
Regarding the shape of the spectrum, while the emission is dominated by photons that escaped from the onaxis region (θ Γ −1 ) at energies near the peak energy and above, the low energy part (ν ≪ ν p ) is dominated by those from the off-axis region. The off-axis component becomes prominent at low energies because the Doppler factor is smaller which leads to a lower peak photon energy ∼ 4D(η 0 , θ)k B T ′ . Therefore, the low energy part of the spectrum can be expressed as a superposition of Blackbody spectra from the off-axis region which have different peak energies (multi-color Blackbody). As a result, the low energy slope of the spectra is somewhat softer than that expected from the Rayleigh-Jeans part of a single blackbody (νL ν ∝ ν 3 ) and can be roughly approximated as νL ν ∝ ν 2.4 .
-Observed luminosity spectrum in the case of a uniform jet (θ 0 = θ 1 = 1 • ) with parameters of the fireball given by L 0 = L 1 = 10 53 erg/s, η 0 = η 1 = 400 and r i = 10 8 cm. The spectrum corresponds to the case when the line of sight to the observer is aligned to the jet axis (θ obs = 0 • ). The red line shows the numerical result. The green, blue and purple lines show the analytical estimates of photospheric emission derived from equation (10) for r = r ph , r = r ph /2 and r = r ph /4, respectively.

Stratified Jet
Here we show the results for a two-component stratified jet. In all cases, the half opening angle of the spine is fixed at θ 0 = 0.5 • .
As mentioned in §3.1, when a uniform jet is assumed, spectra tend to be thermal-like with slight modifications from blackbodies originating from a sphere with radius r = r ph0 /5. On the other hand, the appearance of the spectrum can deviate significantly from a thermal one when a strong velocity shear is present in the outflow, since photons which cross the shear flow multiple times can gain energy through a Fermi-like acceleration mechanism. This can be understood as follows. Under the assumption of elastic scattering, the energy gain of photons in a single scattering event can be expressed as where ν in (ν sc ) and θ in (θ sc ) are the frequency and angle between the the fluid velocity and photon propagation direction before (after) the scattering, respectively. Hence, if θ sc < θ in , photons gain energy and vice versa. Photons which have crossed the boundary layer from the sheath to the spine region tend to gain energy when they are scattered there (upscatter). This is simply because the photons within the sheath region tend to have larger angle between their propagation direction and fluid velocity than those in the spine region. On the other hand, photon which have crossed the boundary layer from the spine to the sheath region tend to lose energy (downscatter) due to the same reason. Consequently, some fraction of photons which cross the boundary layer multiple times can gain energy, since the energy gain by the upscattering overcome the downscattering in average. This mechanism can give rise to a non-thermal spectrum at the high frequencies.
To obtain a rough estimation of the average energy gain and loss rate (ν sc /ν in ) for each process, let us approximate the radially expanding spine and sheath regions as a plane parallel flow. Under the above consideration, the typical angle between the photon propagation direction and the fluid velocity direction for the photons in the spine (sheath) region can be roughly estimated as Since the angle θ v is conserved along the photon's path in the case of a plane parallel flow, the typical energy gain rate by the upscattering in the spine region can be evaluated by substituting θ in = θ v 1 ∼ Γ −1 1 and θ sc = θ v 0 ∼ Γ −1 0 in equation (11) and is given as Similarly, the typical energy loss rate by the downscattering in the sheath region is given as From the above equations, it is clear that the energy gain by the upscattering overcomes the energy loss by the downscattering ( ν sc /ν in up ν sc /ν in down ∼ (1/4)[2 + (Γ 0 /Γ 1 ) 2 + (Γ 1 /Γ 0 ) 2 ] > 1). It is also clear that the efficiency of the acceleration per each cycle of crossing ν sc /ν in up ν sc /ν in down ∼ (1/4)[2+(Γ 0 /Γ 1 ) 2 +(Γ 1 /Γ 0 ) 2 ] is controlled by the ratio between the bulk Lorentz factor of the two regions Γ 0 /Γ 1 and increases as the ratio becomes larger. It is worth noting that, while the average value of the energy ratio ν sc /ν in roughly obeys the above equations, the dispersion around the average value is large, since it depends quite sensitively on the scattering angles (θ in and θ sc ; see Eq. (11)). When a photon from the sheath region that has an angle θ in = f 1 Γ −1 1 is scattered in the spine region with an angle θ sc = f 0 Γ −1 0 , the energy gain by the scattering can be written as . It is clear from the above equation that a small change in the scattering angles (θ in = θ v 0 ∼ Γ −1 0 and θ sc = θ v 1 ∼ Γ −1 1 ) leads to a quite large change in the energy ratio. For example, in the case of f 0 = 0 and f 1 = 2, the energy ratio due to the upscattering is larger than the typical value by a factor of (ν sc /ν in )( ν sc /ν in up ) −1 ∼ 2[1 + 4(Γ 0 /Γ 1 ) 2 ][1 + (Γ 0 /Γ 1 ) 2 ] −1 ∼ 8. Note also that, once the photon energy (evaluated in the electron rest frame) approaches close to the electron rest mass energy hν cmf ∼ m e c 2 , where m e is the electron rest mass, the scattering can no longer be approximated as elastic, since recoil effect becomes non-negligible (Klein-Nishina effect). In this case, the acceleration efficiency is significantly reduced.
In Fig. 3, we display the obtained result for the case of a stratified jet with η 0 = 400 and L 0 = 10 53 erg/s for the spine and η 1 = 200 and L 1 = (η 1 /η 0 )L 0 = 5 × 10 52 erg/s for the sheath. As mentioned in §2, the injection radius is set at a position where a velocity shear between the two regions develops (r inj = r s1 ). The corresponding optical depth is τ (r inj ) ∼ 100 for the spine and τ (r inj ) ∼ 180 for the sheath. The various lines in the figure show the cases for the observer angle with respect to the jet axis being θ obs = 0 • (red), 0.25 • (green), 0.4 • (blue), 0.5 • (purple), 0.6 • (light blue) and 0.75 • (black). As we can see, the spectrum varies quite sensitively with the observer angle. The spectrum for θ obs = 0 • is thermal-like and nearly identical to that obtained in the case of a uniform jet (Fig. 2). The reason for this is simple. Since most of the scattered photons propagate in a direction within a cone of half opening angle ∼ 1/Γ ∼ 0.14 • (Γ/400) −1 , the majority of the observed photons are from a region of θ 0.14 • . Hence, only a small fraction of photons from the sheath region and the boundary (θ ≥ θ 0 = 0.5 • ) can reach the observer, so that the spectrum does not deviate largely from the case of uniform jet. On the other hand, if the observer angle is larger, photons from the sheath and boundary layer become observable. As a result, a non-thermal component appears above the peak energy of the thermal spectrum due to the photon acceleration in the boundary layer. The non-thermal component is hardest when the observer angle is aligned to the boundary layer θ obs = θ 0 = 0.5 • and becomes softer as the deviation between θ obs and θ 0 becomes larger, simply because the boundary layer corresponds to the site of photon acceleration. As mentioned earlier, the photon acceleration becomes inefficient when the photon energy becomes large enough so that the recoil of electrons cannot be neglected (Klein-Nishina effect). Hence, in all cases, the spectrum does not extend up to energies higher than hν ∼ Γ 0 m e c 2 ∼ 200(Γ 0 /400) MeV.
Note also that the peak energy and the luminosity of the thermal component differs enormously for θ obs > θ 0 and for θ obs < θ 0 , due to the differences in the assumed parameters in the spine and sheath regions. For an observer at θ obs < θ 0 , the thermal component is determined mainly by photons which have propagated through the spine region. Therefore, the observed spectrum is nearly identical to the case of the uniform jet considered above in which a same set of parameters (η 0 , L 0 and r i ) is assumed. On the other hand, for an observer at θ obs > θ 0 , photons which have propagated through the sheath region dominate the thermal component. Accordingly, the peak energy and luminosity are lower by a factor ∼ (η 0 /η 1 ) 8/3 (L 0 /L 1 ) −5/12 ∼ 4.7 and ∼ (η 0 /η 1 ) 8/3 (L 0 /L 1 ) 1/3 ∼ 8, respectively.
In Fig. 4, we display the results obtained for η 0 = 200 and L 0 = 10 52 erg/s for the spine and η 1 = 100 and L 0 = 5 × 10 51 erg/s for the sheath. That is, the terminal Lorentz factor of the outflow is smaller by a factor of 2 for both the spine and the sheath than those assumed in the previous case. The optical depths at the injection radius are τ (r inj ) ∼ 170 for the spine and τ (r inj ) ∼ 290 for the sheath. As in the previous case, the non-thermal component is hardest when θ obs = θ 0 = 0.5 • and becomes softer as the deviation between θ obs and θ 0 becomes larger. However, the major difference with the previous case is that the observer dependence of the hardness is weaker. For example, as shown in Fig. 4, the non-thermal component can be prominent even for θ obs = 0 • . This is because beaming effect is weaker than the previous case due to the smaller values of the Lorentz factor, so that the photons can spread out in wider angles.
The dependence of the spectrum on the difference between the dimensionless entropies (terminal Lorentz factor) of the spine and sheath is displayed in Fig. 5. Each panel corresponds to the result for the observer angle fixed at θ obs = 0.25 • (top panel), θ obs = 0.4 • (middle left), θ obs = 0.5 • (middle right), θ obs = 0.6 • (bottom left) and θ obs = 0.75 • (bottom right). In all cases, the dimensionless entropy and kinetic luminosity of the spine are chosen to be η 0 = 400 and L 0 = 10 53 erg/s, respectively. The red line shows the case for a uniform jet, while the green, blue, purple, light blue and black lines show the cases for a sheath with dimensionless entropies of η 1 = 300, η 1 = 250, η 1 = 200, η 1 = 150 and η 1 = 100, respectively. The kinetic luminosity of the sheath is determined by L 1 = (η 0 /η 1 )L 0 .
As mentioned earlier, the peak energy and luminosity of the thermal-component depend on the dimensionless entropy and the kinetic luminosity as E p ∝ η 8/3 L −5/12 and luminosity L p ∝ η 8/3 L 1/3 . Hence, for observer mainly seeing the photons from the sheath region (θ obs > θ 0 ), these values show considerable decrease in models assuming smaller η 1 as is seen in the figure. The nonthermal component becomes significant as η 1 becomes smaller. This tendency is due to the following reasons. One is simply because the bulk Lorentz factor of the sheath becomes smaller for smaller η 1 . As a result, the ratio between the bulk Lorentz factor of the spine and sheath Γ 0 /Γ 1 becomes larger which in turn leads to an increase in the energy gain per each crossing as explained earlier in this section. In addition, the wider spreading of the photons propagating in the sheath region due to the increase in the beaming angle ∼ Γ −1 1 increases the chance for the photons to cross the boundary layer from the sheath to the spine region. Another reason is that, for smaller value of η 1 , the radius where the velocity shear begins to develop r s1 = η 1 r i becomes smaller. This also leads to an increase in the probability for the photons to be accelerated, since the optical depth of the acceleration region (r > r s1 ) increases (see equation (3)).

Relation Between Photon Energy and Number of
Crossings To demonstrate that photons accelerated via the multiple crossing of the spine and sheath boundary layer are indeed the origin of the non-thermal component, we analyzed the relation between the photon energy and the number of crossings that the corresponding photons have experienced. Here the number of crossings, n cr , is defined as the total number of events that the photon has crossed the boundary layer (either from spine to sheath or sheath to spine) before it reaches the outer boundary r out .
In Fig. 6, we show the distribution of the average number of crossings for a given observed photon energy, n cr ν . On the other hand, distribution of the average observed energy for a given number of crossings, ν ncr , is displayed in Fig. 7. The two cases of stratified jet that are displayed in the figures by the purple (Case I: η 0 = 400 and η 1 = 200) and black lines (Case II: η 0 = 400 and η 1 = 100) correspond to the analysis of photons shown in Fig. 5 using same colors.
From Fig. 6, it is confirmed that the photons at higher energies tend to have larger number of crossings. While the photons below the thermal peak energy (hν MeV) do not require multiple crossings, the prominent nonthermal component extending above ∼ 1 MeV is produced by the photons that cross the boundary layer ∼ 10 − 15 times in average. Comparing the two cases, the overall distribution of n cr ν . does not vary much. The difference in the average number of crossing is within ∼ 2 in all energies up to ∼ 100 MeV. Note that, however, that this does not imply that the average energy for a given number of crossings ν ncr does not vary much in the two cases. Conversely, the difference in ν ncr is quite large between the two cases as is seen in Fig. 7, since the average energy gain per crossing is quite sensitive to the ratio in the terminal Lorentz factor Γ 0 /Γ 1 (see Eqs. (12) and (13)). However, due to the large dispersion in the energy gain per crossing, the energy distribution of photons does not show a sharp peak at the average energy ν ncr but extends to energies below and above ν ncr by many orders of magnitude. As a result, the distribution of n cr ν do not directly reflect the distribution of ν ncr , since the photons which have crossed the boundary layer a certain number of times can dominate over the other population of photons in wide energy ranges. To clarify this, we show the energy distribution of the photon number count rate, νN ν [photons/s], for a given number of crossings in 8. From the figure, it is confirmed that, while there is a significant discrepancy in ν ncr , the photons with number of crossing n cr ∼ 10 tend to dominate the population in the energy range in which nonthermal component becomes prominent (hν few MeV) in both cases. For this reason, the resultant distribution of n cr ncr does not vary much in the two cases.
Regarding the distribution of the average energy, ν ncr tends to increase with the increasing number of crossing n cr initially and then approaches a constant value. This asymptotic behaviour is due to the Klein-Nishina effect. As the photon energy becomes large and exceeds 100 keV in the comoving frame, acceleration efficiency is reduced by the effect (see §3.2 for detail), and the average energy can no longer increase. As in seen in Fig. 7, the dependence of ν ncr on n cr is not smooth but rather bumpy. This reflects the fact that the photons tend to be upscattered when crossing from the sheath to spine region occurs (see Eq. (12)), while, on the other hand, tend to be downscattered when crossing in the opposite direction occurs (see Eq. (13)). Therefore, the energy of individual photons is not a monotonically increasing function of n cr but shows a bumpy dependence, since upscattering and downscattering occurs alternately in each crossing events. Although somewhat reduced, this feature remains even after averaging up and leads to the appearance wiggles in the distribution of ν ncr .
To quantify the energy gain and loss in each crossing, we display the ratio of the average energy for photons with n cr crossings to that for photons with n cr − 1 crossings ( ν ncr / ν ncr−1 ) in Fig. 9. In addition to the analysis of the total photons, we also display the results of the analysis for the two populations of photons that were initially injected in the spine region (θ inj ≤ θ 0 ; blue line) and sheath regions (θ inj > θ 0 ; green line). In the former case (θ inj ≤ θ 0 ), if n cr is an odd (even) number, the number of crossings from the spine (sheath) to sheath (spine) region is greater by 1 than that for n cr − 1 crossings. Hence, the number of downscattering (upscat-tering) event is greater by 1 for the photons with an odd (even) number of n cr than for those with n cr − 1 crossings. As a result, ν ncr / ν ncr−1 is less (greater) than unity for the odd (even) number of n cr . On the other hand, as is obvious, the opposite is true for the latter case (θ inj > θ 0 ). Jets with a larger difference in the terminal Lorentz factor (Case II; right panel) show a larger range of energy ratio ν ncr / ν ncr−1 than those with a smaller difference (Case I; left panel), since the efficiency of the single upscattering (downscattering) increases (decreases) as the relative difference in the Lorentz factor becomes larger.
Regarding the upscatterings, the energy ratio ν ncr / ν ncr−1 for photons that cross the boundary from the sheath to spine region (green line) only once (n cr = 1) is relatively small because a large fraction of these photons experience the crossing when the velocity shear is not fully developed (r < r s0 ). For photons with a larger number of crossings (n cr 2), a large fraction of the photons experience the last crossing at r > r s0 , where the velocity shear is fully developed. Therefore, the energy ratio is larger than that for n cr = 1 and is roughly constant as long as the Klein-Nishina effect is negligible. When the Klein-Nishina effect becomes important (n cr 10 for Case I and n cr 5 for Case II), ν ncr / ν ncr−1 decreases as n cr increases, and again asymptotically approaches a constant value. On the other hand, regarding the case of downscatterings, the energy ratio is relatively insensitive to the number of crossings in both cases and is roughly in the range ∼ 0.25 − 0.6. It is worth noting that the values of energy ratio are consistent within a factor of ∼ 2 with the rough estimations given by Eqs. (12) and (13). Above the saturation radius (r ≥ r s0 ), the equations predict an upscattering (downscattering) with energy ratio of ∼ 2.5 (∼ 0.6) and ∼ 8.5 (∼ 0.5) for Cases I and II, respectively. The energy ratio of the total photons is larger and smaller when the photons from the spine correspond to the upscattering (even n cr ) and downscattering (even n cr ), respectively, since more photons originate from the spine region than from the sheath region.
Lastly, to obtain a further insight into the relation between the photon acceleration and number of crossing, we show the evolution of the average photon energy evaluated in the comoving frame, ν cmf , with radius in Fig 10. The red, green, blue, purple, light blue, yellow and black lines display the photons that have experienced n cr = 0, 3, 5, 10, 15, 20 and 25 crossings, respectively. For comparison, we also plot the curve of ν cmf ∝ r −2/3 with a thin black line, which corresponds to the adiabatic cooling expected above the saturation radius (r ≥ r s0 ).
Regarding the photons that do not experience any crossings (n cr = 0), the overall evolution of the energy is determined solely by the adiabatic cooling due to the expansion of the jet, since photon acceleration does not take place. Below the saturation radius of the spine r ≤ r s0 , the cooling rate in the spine region (ν cmf ∝ r −1 ) is higher than that in the sheath region (ν cmf ∝ r −2/3 ), since the spine region is in the acceleration phase, while the sheath region is in the coasting phase (r ≥ r s1 ). Hence, the energy evolution of the average comoving energy is in between the two cooling rates. At larger radii (r ≥ r s0 ), since both regions are in the coasting phase, the average energy evolves as ν cmf ∝ r −2/3 un-til the coupling between the photon and matter becomes weak (r 0.1r ph0 ). The cooling rate gradually reduces at r 0.1r ph0 due to the weak coupling, and the comoving energy approaches a constant value. It is worth noting that the behaviour of the average comoving energy for photons above the saturation radius is consistent with that found in the previous studies (e.g., Pe'er 2008;Beloborodov 2011;Bégué et al. 2013).
Regarding the photons with at least one crossing (n cr ≥ 1), the evolution of the comoving energy cannot be described only by the adiabatic cooling due to the presence of the photon acceleration. As the number of the crossings increases, the departure from the simple adiabatic cooling becomes significant and the comoving energies tend to be larger. For a given number of crossings, the departure is more prominent in Case II than in Case I due to the increase in the acceleration efficiencies. When the number of the crossings is sufficiently large so that the effect of Klein-Nishina becomes non-negligible, the evolution of the average comoving energy asymptotically approaches a single curve since the acceleration saturates. This tendency is clearly seen in Fig. 10 (for example, see lines that display the evolution of photons with n cr ≥ 10 in Case II).

Comparison with Observations
So far, we have shown that in the presence of velocity shear, photons can be accelerated at the shear region producing a high energy power-law component above the thermal peak energy. From observations, the observed spectra of GRB prompt emission can be often modeled by a Band function, which is a smoothly joined broken power-law that peaks at ∼ a few 100 keV. The photon indices below and above the peak (α ph and β ph ) vary from source to source (νL ν ∝ ν α ph +2 for ν < ν p and νL ν ∝ ν β ph +2 for ν > ν p ). Focusing on the high-energy slope β ph , the measured values are roughly in the range between −2 and −3, with a typical value at ∼ −2.5. The power-law component predicted in our model can accommodate various spectral indices depending on the flow parameters as well as observer angle. Hence, by adopting appropriate values for our parameter set, the observed high energy spectral slope (νL ν ∝ ν −0.5±0.5 ) as well as the harder and softer slopes observed in some GRBs can be reproduced as shown in Figs. 3-5. Interestingly, our model can also reproduce spectral features seen in some peculiar GRBs such as GRB 090510 (Ackermann et al. 2010), 090902B (Abdo et al. 2009) and 090926A (Ackermann et al. 2011). In these bursts, in addition to a Band-type (or thermal-like) component, an extra hard power-law component (photon index larger than −2) is required to model the overall spectrum. As shown in Fig. 5, in some cases, especially when the difference between η 0 and η 1 is large, the non-thermal power-law component does not join the thermal component smoothly at the peak energy, but appears as a hard power-law component which becomes prominent at energies above the bump of the thermal component. These spectral features resemble those found in these bursts. Hence, we emphasize that the spectra of peculiar bursts may be due to the relatively large difference in η (or bulk Lorentz factor) in the shear region. It is noted, however, that the extra power-law component that also extends to energies below the thermal-like component found in GRB 090902B and 090510 is hard to explain. If present, different emission components such as synchrotron emission as discussed in Pe'er et al. (2012) may be required in order to explain the low energy end of the power-law component. We also note that, in order to extend the powerlaw component up to GeV energies, quite large Lorentz factors are required, since the non-thermal component extends only up to energy of hν 200(η/400) MeV due to the Klein-Nishina effect. A detailed discussion on this issue is given in §4.6.
On the other hand, the spectral slope below the peak The red line shows the case of a uniform jet, while the green, blue, purple, light blue and black show the cases for the dimensionless entropy of the sheath given by η 1 = 300, η 1 = 250, η 1 = 200, η 1 = 150 and η 1 = 100, respectively. In each case, the kinetic luminosity of the sheath is given by L 1 = (η 0 /η 1 )L 0 . energy is only moderately sensitive to the values of the parameters and can be well approximated by a blackbody emission from the off-axis region (|θ − θ obs | Γ −1 ) of a sphere at a radius r ∼ r ph /5 as in the case of a uniform jet. In all cases, the low energy photon indices are roughly α ph ∼ 0.5 which is harder than the typical observed value (α ph ∼ −1). It is noted, however, that this does not imply that the low energy slope cannot be reproduced within this scenario. For instance, while the instantaneous spectrum is hard, the time integration can lead to a softer spectrum if the evolution history of the outflow is considered. To demonstrate this, let us con-sider the case when the dimensionless entropies of the spine-sheath jet have evolved from η 0 = 400 and η 1 = 200 (initial stage) to η 0 = 200 and η 1 = 100 (later stage), while the kinetic luminosity is fixed as L 0 = 10 53 erg/s and L 1 = 5 × 10 52 erg/s. We illustrate the resultant spectrum in Fig. 11 under the assumption that a nearly equal time period is spent in the two stages. The black solid line shows the overall spectrum, and the red and green solid lines show the contribution from the initial and later stages, respectively. Since the peak energy and the luminosity vary with the dimensionless entropy of the flow roughly as ∝ η 8/3 , they are smaller by a factor ∼ 6, θ 0 =0.5° θ 1 =1°L 0 =10 53 erg/s L 1 =(η 1 /η 0 )L 0 r i =10 8 cm Fig. 6.-Average number of crossing events of the spine-sheath boundary layer as a function of observed photon energy for a jet in which the spine region with half opening angle of θ 0 = 0.5 • is embedded in a wider sheath outflow with half opening angle of θ 1 = 1 • . The employed values for dimensionless entropy (terminal Lorentz factor) and kinetic luminosity are chosen as η 0 = 400 and L 0 = 10 53 erg/s for the spine. The purple and black lines show the cases of the dimensionless entropy of the sheath given by η 1 = 200 and η 1 = 100, respectively. In each case, the kinetic luminosity of the sheath is given by L 1 = (η 0 /η 1 )L 0 .  Fig. 7.-Average energy of the photons as a function of the number of crossings of the spine-sheath boundary layer ncr for a jet in which the spine region with half opening angle of θ 0 = 0.5 • is embedded in a wider sheath outflow with half opening angle of θ 1 = 1 • . The employed values for dimensionless entropy (terminal Lorentz factor) and kinetic luminosity are chosen as η 0 = 400 and L 0 = 10 53 erg/s for the spine. The purple and black lines show the cases of the dimensionless entropy of the sheath given by η 1 = 200 and η 1 = 100, respectively. In each case, the kinetic luminosity of the sheath is given by L 1 = (η 0 /η 1 )L 0 . in the later stage. As a result, while the peak energy and luminosity of the overall spectra is determined by the initial stage, at lower energies, contribution from the later stage becomes dominant. Due to the superposition of the two components, the spectrum below the peak becomes softer than that of the individual component, so that a typical observed spectrum (νL ν ∝ ν) is reproduced. Note also that the high energy part of the overall spectrum largely resembles typical ones from observations. Hence, we conclude that typical observed spectra can be successfully reproduced when a time evolution is considered. 10 4.2. On the Observer Angle Dependence and Structure of the Jet As shown in the previous section, our results have a quite strong dependence on the observer angle θ obs . While strong non-thermal emission can be seen when the observer angle is close to the angle in which a velocity shear is present (|θ 0 − θ obs | Γ −1 ), the non-thermal feature tends to be weaker for observer angles far from θ 0 . This tendency seems to contradict the fact that most of the observed GRBs have non-thermal features in their spectra. However, we emphasize that this difficulty can be overcome if the jet possesses a more complex structure. For example, if the jet has more than two components and velocity shear is present at multiple angles more closely spaced than Γ −1 , photons from the acceleration regions will be prominent for all observers lying within the opening angle of the jet (θ obs ≤ θ 1 ). Even in the case of a simple spine-sheath jet, if the angular boundary of the spine and sheath (θ 0 ) varied rapidly with time, the accelerated non-thermal photons would be observable across a broad range of angles. It should be also noted that the structure need not be in the direction of θ. Strong velocity shear in the azimuthal (φ) direction and/or radial direction due to the presence of turbulence or shocks can also provide an acceleration sites for the photons (Bromberg et al. 2011;Ioka et al. 2011). Any structure showing strong velocity shear within an angle Γ −1 from the line-of-sight can give rise to a non-thermal component above the thermal peak. Hence, it is expected that jet having a rich structure and/or rapid time variability will be naturally accompanied by non-thermal emission, irrespective of the observer angle. The origin of the structure and variability could be due to the nature of the central engine (McKinney 2006;Nagataki et al. 2007;Nagataki 2009;McKinney & Blandford 2009;Nagataki 2011) and/or the propagation of jet through the envelope of the progenitor star (Zhang et al. 2003;Mizuta et al. 2006;Morsony et al. 2007;Lazzati et al. 2009;Mizuta et al. 2011;Nagakura et al. 2011).

4.3.
On the E p -L p relation Our calculations have shown that the peak energy of the observed spectra can be roughly approximated as the thermal peak of a blackbody emission from the surface of optical depth τ (r) ∼ 2 (r ∼ r ph /2). On the other hand, the peak luminosity shows rough agreement with that of the emission from τ (r) ∼ 1 (r ∼ r ph ). This result is valid for an adiabatic fireball which has a photosphere above the saturation radius (r ph > r s ) (see §3.1 for detail). As a result, the peak energy and luminosity are given by θ 0 =0.5° θ 1 =1° η 0 =400 η 1 =100 r i =10 8 cm L 0 =10 53 erg/s L 1 =2.5×10 52 erg/s n cr =0 n cr =3 n cr =5 n cr =10 n cr =15 n cr =20 n cr =25 -Ratio of the average observed energy for photons which cross the spine-sheath boundary layer (θ = θ 0 ) ncr times to those which experience ncr − 1 crossings for a jet in which a spine region with half opening angle of θ 0 = 0.5 • is embedded in a wider sheath outflow with half opening angle of θ 1 = 1 • . The employed values for dimensionless entropy (terminal Lorentz factor) and kinetic luminosity are chosen as η 0 = 400 and L 0 = 10 53 erg/s for the spine. The left and right panels show the cases of the dimensionless entropy of the sheath given by η 1 = 200 and η 1 = 100, respectively. In each case, the kinetic luminosity of the sheath is given by L 1 = (η 0 /η 1 )L 0 . The energy ratio for total photons are displayed in the left and right panels with the purple and black lines, respectively. In addition, also the energy ratio for the photons that were initially injected in the spine region (θ inj ≤ θ 0 ) and sheath region (θ inj > θ 0 ) are shown by blue and green lines, respectively, in both panels. respectively 11 . Derived from observations, there is an empirical relation between the peak energy and luminosity (Yonetoku et al. 2004;Kodama et al. 2008;Nava et al. 2011;Amati et al. 2002;Wei & Gao 2003) that is roughly given by Therefore, in order to reproduce the empirical relation from the photospheric emission, the parameters of the 11 For cases where efficient energy dissipation is present within the flow, the dependence of Ep and Lp on the fireball parameters can be significantly different (e.g., Giannios 2012;Lazzati et al. 2013). It is worthy to note that comparison of Eqs. (14) and (15) Another important ingredient for comparison with observations is the emission efficiency η γ = L p /L. Observations of the afterglows suggest a quite high GRB efficiency in the range 0.01 η γ 1, with a typical value at η γ ∼ 0.1 − 0.2 (e.g., Fan & Piran 2006;Zhang et al. 2007). This gives another constraint on the fireball parameters which can be written as 0.01 0.1r To sum up, the observed E p -L p relation as well as the high efficiency can be reproduced when equations (17) and (18) are satisfied. In Fig. 12, we show the calculated spectra for a uniform jet (θ 1 = θ 0 = 1 • ) in which both conditions are  η 0 =400 η 1 =100 comoving frame L 0 =10 53 erg/s L 1 =2.5×10 52 erg/s n cr =0 n cr =3 n cr =5 n cr =10 n cr =15 n cr =20 n cr =25 Fig. 10.-Average comoving energy of photons as a function of radius for a jet in which a spine region with half opening angle of θ 0 = 0.5 • is embedded in a wider sheath outflow with half opening angle of θ 1 = 1 • . The red, green, blue, purple, light blue, yellow and black lines display the photons that have experienced ncr = 0, 3, 5, 10, 15, 20 and 25 crossings, respectively. The thin black solid line plots the curve of ∝ r −2/3 , which corresponds to the adiabatic cooling expected above the saturation radius (r ≥ r s0 ). The employed values for dimensionless entropy (terminal Lorentz factor) and kinetic luminosity are chosen as η 0 = 400 and L 0 = 10 53 erg/s for the spine. The left and right panels show the cases of the dimensionless entropy of the sheath given by η 1 = 200 and η 1 = 100, respectively. In each case, the kinetic luminosity of the sheath is given by L 1 = (η 0 /η 1 )L 0 .
total Fig. 11.-Observed luminosity spectrum in the case of spinesheath jet in which a spine jet with half opening angle of θ 0 = 0.5 • is embedded in a wider sheath outflow with half opening angle of θ 1 = 1 • . The observer angle is fixed as θ obs = 0.4 • . The black solid line shows the overall spectrum for the case when the dimensionless entropy of the spine-sheath jet have evolved from η 0 = 400 and η 1 = 200 (initial stage) to η 0 = 200 and η 1 = 100 (later stage), while the kinetic luminosity is fixed as L 0 = 10 53 erg/s and L 1 = 5 × 10 52 erg/s. The red and green solid lines show the contribution from the initial and later stages, respectively.
fulfilled. The adopted values of the parameters are summarized in the figure. The red, green, blue and purple lines correspond to the cases of η 0 = 100, η 0 = 200, η 0 = 400 and η 0 = 500, respectively. The remaining parameters (L 0 and r i ) are determined so that an efficiency of η γ ∼ 0.1r

4.4.
Dependence on the width of the spine-sheath boundary layer In the present study, we have assumed an infinitesimal width for the boundary layer of the spine and sheath region. However, in reality, the boundary layer is expected to possess a finite width due to the interaction between the two regions. Firstly, by definition, photons are closely coupled to the matter below the photosphere, and, therefore, will couple the two regions due to Compton friction (radiative viscosity). In particular, this effect will be important in the regions where the energy is dominated by radiation (r r s ) (e.g., Arav & Begelman 1992). Secondly, even in the absence of radiation coupling, the Kelvin-Helmholtz instability should grow whenever velocity shear is present (e.g., Turland & Scheuer 1976;Bodo et al. 2004). These effects will relax the discontinuous change in the velocity and lead to broadening of the boundary layer. Although the detail analysis of the resultant jet structure is beyond the scope of the present study, these effects will reduce the acceleration of the photons. To quantify the effect of the broadening of the boundary layer on the spectra, here we compute the photon propagation within a spine-sheath jet having a boundary layer with finite width in transverse direction. We explore the dependence on the width of the boundary layer by considering cases with various widths.
In modeling the jet structure, we added slight modifications in the original spine-sheath jet model considered in the present study. The boundary layer is defined as a region having finite transverse width dθ B that is located in the range θ i < θ < θ ii , where θ i = θ 0 − dθ B /2 and θ ii = θ 0 + dθ B /2. Correspondingly, the spine and sheath regions are limited to the ranges 0 ≤ θ ≤ θ i and θ ii ≤ θ ≤ θ 1 , respectively. While the fluid properties (fireball parameters) of the spine and sheath region are determined in the same way as in the cases of infinitesimal boundary layer ( §2.1.2), the properties of the boundary layer are determined by simply imposing a linear interpolation of the fireball parameters from the two regions. Hence, the initial radius of the fireball is fixed -Observed luminosity spectrum for a spine-sheath jet in which the width of the boundary layer dθ B is taken to be nonzero. The the spine and sheath regions are located in the transverse range 0 ≤ θ ≤ θ i = θ 0 − dθ B /2 and θ ii = θ 0 + dθ B /2 ≤ θ ≤ θ 1 , respectively. In all cases, fixed values of θ 0 = 0.5 • and θ 1 = 1 • are employed and the observer angle is fixed at the midpoint of the boundary layer θ obs = θ 0 = 0.5 • . The employed values for dimensionless entropy (terminal Lorentz factor) and kinetic luminosity are chosen as η 0 = 400 and L 0 = 10 53 erg/s for the spine and η 0 = 200 and L 0 = 5 × 10 52 erg/s for the sheath. While the thick red line displays the case for an infinitesimal width (dθ B = 0), the thin green, blue, purple, light blue, yellow, black and red lines display the cases for widths of in all regions (r i = 10 8 cm), and the transverse distribution of the dimensionless entropy and kinetic luminosity within is given by and respectively. Hence, the profile of the velocity (Lorentz factor) and density within the jet is continuous in all regions. Having determined the background fluid properties, the propagation of photons is calculated in the same way as in the case of an infinitesimal boundary layer. Black body emission is injected at the surface of a fixed radius r in = r s1 with a comoving temperature determined from the fireball parameters ( §2.2.1). Then, photons are propagated until they reach the outer boundary of the calculation range (see §2.2 for detail).
In Figs. 13 and 14, we show the obtained spectra for a stratified jet in which the fireball parameters for the spine and sheath regions are identical to those employed in Case I (η 0 = 400 and η 1 = 200) and Case II (η 0 = 400 and η 1 = 100) discussed in §3.3. The various lines reflect the width of the boundary layer dθ B . While the thick red line displays the case for the infinitesimal width (dθ B = 0), the thin green, blue, purple, light blue, yellow, black and red lines display the cases for the finite widths of dθ B = (100η 0 ) −1 , (20η 0 ) −1 , (10η 0 ) −1 , (5η 0 ) −1 , (2η 0 ) −1 , η −1 0 and (η 0 /2) −1 , respectively. In all cases, the observer angle is fixed at θ obs = θ 0 = 0.5 • . As expected, the non-thermal component becomes softer as the width of the boundary layer broadens due to the reduction in efficiency of photon acceleration. Note that, while the broadening of the boundary layer causes the spectra to depart from the high energy tail of the Band function for Case I, good agreement is found at dθ B ∼ (5η 0 ) −1 for Case II. Therefore, although a large gradient in the Lorentz factor is required, we emphasize that the broadening of the boundary layer does not rule out the photon acceleration mechanism in a stratified jet as the origin of the high energy spectra in the prompt emission.
From the figures, it is seen that, in order to provide an efficient acceleration site comparable to the case of boundary layer with infinitesimal with, the bulk Lorentz factor must vary by a factor of few roughly within an angle dθ B ∼ (100η 0 ) −1 ≪ Γ −1 . As is expected, this condition is roughly equivalent to the condition for the boundary layer to be optically thin above the radius where the velocity shear begins to develop (r ≥ r s1 ). This can be shown as follows. The typical angle between the photon propagation direction and the fluid velocity (radial) direction is roughly θ v ∼ Γ −1 . Hence, at a given radius r, the typical distance that the photons must propagate to cross the boundary layer is roughly given as l cr ∼ r(dθ B / θ v ) ∼ rΓdθ B . The above estimate is valid as long as dθ B ≪ Γ −1 . On the other hand, the mean free path for the typical photon at a given radius r is roughly given by l mfp ∼ D( θ v ) −1 σ T n e (r) ∼ r/τ (r). Thus, the typical optical depth for a photons to cross a boundary layer of finite width can be estimated as τ cr ∼ l cr /l mfp ∼ τ (r)Γdθ B . Since τ (r s1 ) and Γ(r s1 ) are ∼ 100 (∼ 600) and 1/2η 0 (1/4η 0 ) for Case I (Case II), respectively, τ cr 1 is obtained at r = r s1 for dθ B (100η 0 ) −1 . As a result, since τ cr decreases as the radius increases, when the condition dθ B (100η 0 ) −1 is satisfied, the boundary layer is optically thin to crossings at r ≥ r s1 .
Lastly, let us briefly comment on a similar calculation performed by Lundman et al. (2013). In their study, photon propagation within a jet having a continuously decaying velocity profile Γ ∝ θ −p at the outer region (θ θ j ) is explored. Since the efficiency of photon acceleration increases as the difference in velocity (Lorentz factor) between the spine and sheath region increases, their calculation should also show non-thermal high energy photons when a large velocity gradient is considered. Indeed, although not so prominent as is shown by our results, signs of photon acceleration are also reported in their study. While only a thermal component is present when a relatively small velocity gradient (p = 2) is assumed, power-law excess above the thermal peak energy appears in the calculated spectra when a large velocity gradient (p = 4) is assumed. Therefore, as the velocity gradient enlarges, we expect to find more efficient acceleration as is seen in the present study.

On the Absorption Processes
Since the photon number density is significant inside the jet (r r ph ), photons are subject to γγ attenuation once they exceed the threshold energy for the process. The threshold energy is given by E th = 2(m e c 2 ) 2 /[(1 − cosθ ph )ǫ], where ǫ and θ ph are the energy of the target photons and angle between the propagation direction of the target and incident photons, respectively. Since the photons are basically advected in the radial direction within a small angle Γ −1 , the collision angle of the photons is expected to be in a range θ ph Γ −1 . Hence, for photons confined in an angle Γ −1 , the minimum value of the threshold energy (θ ph ∼ Γ −1 ) can be roughly esti-mated as E th,min ∼ 100 (Γ/100) 2 (ǫ/100 MeV) −1 MeV. For the cases considered in the present study (Γ ≥ 100), the photon energies do not exceed 100 MeV by much and, therefore, the pair creation process is negligible for the majority of the photons. However, for the small fraction of photons which propagate with a large angle ( Γ −1 ) with respect to the radial direction, this effect can become non-negligible. In particular, pair annihilation may play an important role for the high energy photons which are being accelerated, since scattering with a large angle with respect to the fluid velocity (radial) direction is favored for gaining photon energy (see §2 for detail). To check this, we have compared our results with those obtained by discarding the photons which have exceeded the threshold energy obtained by substituting the angle between the photon propagation direction and radial direction for θ ph and the highest photon energy appear in the calculation ∼ Γ 0 m e c 2 for ǫ. Indeed, we find that only a small fraction of photons are absorbed and the change in the spectrum is negligible. Therefore, we conclude that γγ attenuation effect does not affect the obtained results.
The flow within the fireball is expected to be fully ionized, and the effect of free-free absorption should also be discussed. For an electron-proton plasma, the frequency averaged free-free opacity can be roughly written as α ff (r) ≈ 1.7 × 10 −25 T ′ (r) −7/2 n e (r) 2 cm −1 (Rybicki & Lightman 1979) in the fluid rest frame. In the laboratory frame, by using a Doppler factor, the opacity can be expressed as α ff,lab = D −1 α ff . Since photons advect in the radial direction within an angle Γ −1 (D ∼ Γ), the optical depth for photons which have propagated from r = r inj to the observer (infinity) can be roughly estimated as τ ff ∼ r inj α ff,lab Γ −1 ∼ 8.3 × 10 −6 (L/10 53 erg s −1 ) 9/8 (η/400) −2 (r i /10 8 cm) 5/4 (r inj /10 10 cm) −5/2 , we have assumed r inj ≤ r s and used equations (2) and (5) in the last equality. Therefore, we can conclude that free-free absorption is also negligible. 4.6. On the GeV Gamma-ray Fermi observations have shown that a fraction (∼ 8%) of GRBs are accompanied by significant emissions at energies well above ∼ 100MeV (e.g., Zhang et al. 2011;The Fermi Large Area Telescope Team et al. 2012). Within the framework of our model, the energy of the photons is limited by the bulk Lorentz factor of the flow as hν obs 50(Γ/100) MeV ( §3). Hence, in order to generate emissions at hν obs ∼ GeV, large bulk Lorentz factors such as Γ > 2000 are required. Alternatively, if we consider the presence of relativistic electrons due to some kind of dissipative processes within the flow (e.g., Ioka et al. 2007;Lazzati & Begelman 2010;Giannios 2012;Beloborodov 2010;Vurm et al. 2011), lower values are allowed for the bulk Lorentz factor. By denoting the maximum Lorentz factor of the electrons measured in the rest frame of the fluid as γ max , the lower limit on the bulk Lorentz factor to produce GeV photons decreases as Γ > 200(γ max /10) −1 . Therefore, within this picture, GRBs with intense GeV emissions may imply the presence of fluid components with very high bulk Lorentz factors (Γ > 2000) or a dissipative process producing relativistic electrons within the flow. In either case, pair cascades due to γγ attenuation may play an important role (Ioka et al. 2011), and the details of the process and its effect on the observed photon spectra are beyond the scope of the present study.
It should be also noted that photons with energies above and below ∼ 100 MeV may have distinct origins. For example, it is widely discussed that ∼ GeV emission can result from energy dissipation by external shocks similar to afterglow emission (Kumar & Barniol Duran 2009, 2010Ghisellini et al. 2010), due to the smooth temporal decay seen in lightcurve of many LAT detected GRBs. Therefore, one possible interpretation is that, while emission at 100 MeV has a photospheric origin as discussed in the present study, higher energy photons are independently produced at the external shock.