VISCOUS BOUNDARY LAYERS OF RADIATION-DOMINATED, RELATIVISTIC JETS. I. THE TWO-STREAM MODEL

Astrophysical jets almost certainly exist as aggregates of massive particles, magnetic fields, and radiation. In certain scenarios, however, the contribution of radiation to the energetics of the outflow far outweighs those of the particles and magnetic fields, meaning that one can essentially neglect the presence of the latter two entities.

One such scenario occurs in the collapsar model of long gamma-ray bursts (GRBs; Woosley 1993; MacFadyen & Woosley 1999). In this model, the core of a massive, evolved star collapses directly (or with a short-lived neutron star phase) to a black hole during the infall stage of an SN II. The energy released by the material accreting onto the black hole, and ultimately observed as the gamma-ray emission, is collimated into bipolar jets—the jet formation confirmed by energetics arguments (Waxman et al. 1998; Fruchter et al. 1999; Frail et al. 2001) and the observations of breaks in the X-ray afterglow light curves (Panaitescu 2007; Dado et al. 2008; Racusin et al. 2009)—and is often sufficient to unbind the stellar envelope, resulting in a supernova (Galama et al. 1998; Bersier et al. 2004; Kamble et al. 2009; Levan et al. 2014, but see Fynbo et al. 2006). If one assumes that the mass of the remnant black hole is on the order of a few solar masses, its accretion luminosity exceeds the Eddington limit by roughly 10 orders of magnitude, meaning that radiation pressure, even if the flux is nearly isotropic, is likely an important mechanism for driving and sustaining the outflow (the fireball model; Rees & Mészáros 1992). Even if the jet is launched by magnetohydrodynamical mechanisms (Blandford & Znajek 1977; Blandford & Payne 1982), radiation could still play a prominent role in determining the dynamics of the jet. Arguments concerning the time necessary for the jet to break through the stellar envelope also seem to disfavor Poynting-dominated jets (Bromberg et al. 2015; but see Mundell et al. 2013).

Jets produced during tidal disruption events (TDEs; Giannios & Metzger 2011)—when a star is destroyed by the tidal force of a supermassive black hole—could provide another class of radiation-dominated outflow. After the star is tidally disrupted, roughly half of the shredded debris remains bound to the black hole and returns to the tidal disruption radius. If the black hole has a mass less than roughly ${10}^{7}{M}_{\odot }$, that rate of return can exceed the Eddington limit of the black hole by orders of magnitude for a significant amount of time (on the order of days to months; Evans & Kochanek 1989; Strubbe & Quataert 2009). Provided that this material can rapidly accrete onto the black hole, which is likely the case due to the tidal dissipation of kinetic energy (Kochanek 1994; Guillochon et al. 2014) and relativistic precession effects (Rees 1988; Evans & Kochanek 1989), the energy released during the accretion process will also be supercritical. It was during this supercritical phase that the event Swift J1644+57 was seen to have an associated jetted outflow (Bloom et al. 2011; Burrows et al. 2011; Cannizzo et al. 2011; Zauderer et al. 2011; the source Swift J2058+05 may provide another example of a jetted, super-Eddington TDE; Cenko et al. 2012). Although the jet launching mechanism for this event is uncertain, the magnetic field of the tidally disrupted star, assuming its flux is approximately conserved, is almost certainly insufficient to power the outflow. Therefore, unless one invokes the existence of a fossil magnetic field (Kelley et al. 2014; Tchekhovskoy et al. 2014), the radiation pressure associated with the accretion luminosity likely plays some role in powering the jet. At any rate, the radiation released during the supercritical accretion process affects the dynamics of the collimated outflow and contributes substantially to its overall energy and momentum.

Collapsar jets inject a significant amount of energy into the overlying stellar envelope as they punch their way into the circumstellar medium, creating a pressurized "cocoon" of shocked material with which the jet interacts (Lazzati et al. 2007; Morsony et al. 2007; López-Cámara et al. 2013), and the progenitors themselves—usually taken to be highly evolved, Wolf–Rayet stars (Matzner 2003; Woosley & Heger 2006)—are likely sustained primarily by radiation pressure. The initial stages of collapsar jet propagation are therefore characterized by the transfer of energy and momentum between two radiation-dominated fluids. Because of the high accretion rates and low specific angular momentum, the fallback disks generated during the super-Eddington phase of TDEs are likely optically and geometrically thick and radiation pressure-supported (Rees 1988; Loeb & Ulmer 1997; Coughlin & Begelman 2014a). In the zero-Bernoulli accretion model of Coughlin & Begelman (2014a), for example, the accreting material is inflated into a quasi-spherical envelope that surrounds the black hole. The dynamics of the jets of supercritical TDEs are therefore also modulated by the presence of a radiation pressure-supported, external medium.

Previous authors used the supersonic propagation of the outflow to model the interaction of the collapsar jet with the overlying envelope as an oblique shock-boundary layer structure (Bromberg & Levinson 2007; Kohler et al. 2012; Kohler & Begelman 2012, 2015). In these models, the outflow is assumed to consist of a perfect, single fluid with a relativistic equation of state, and some have included the presence of magnetic fields (e.g., Kohler & Begelman 2012). A more realistic picture, however, is obtained by considering the jet as a composite of the massive scatterers present in the outflow and the radiation that accompanies it. Indeed, this approach constitutes the underlying framework of radiation hydrodynamics. Furthermore, non-ideal, i.e., viscous, effects will tend to "smear" the discontinuities otherwise present in the system, resulting in a more gradual transition of the fluid quantities between the jet and the environment.

In the limit that the mean free path of a photon is small, radiation acts like an effective viscosity, with a coefficient of dynamic viscosity that depends both on the radiation pressure and the density of scatterers (see Section 2 of this paper), and transfers momentum and energy between neighboring fluid elements (Weinberg 1971; Loeb & Laor 1992). In the radiation-dominated interaction between collapsar jets and TDE jets and their respective ambient media, the effects of radiation viscosity should be quite large. Arav & Begelman (1992) considered the effects of radiation viscosity on the evolution of boundary layers in the two-stream approximation (see Section 3 of this paper); their treatment, however, was of a non-relativistic nature, meaning that their results have limited applicability to GRB and TDE jets.

In this paper, the first of two, we consider the effects of radiation viscosity on the propagation of relativistic jets in radiation-rich environments. In Section 2 we present the equations of radiation hydrodynamics in the viscous limit. Section 3 applies those equations to the two-stream problem, wherein the jet and ambient medium are considered to be two separate, interacting fluids, and we compare our results to the non-relativistic treatment of Arav & Begelman (1992). In Section 4 we discuss the results of the analysis and comment upon the application of our models to the jets produced by supercritical TDEs such as Swift J1644+57, GRBs, and other astrophysical systems. In a second paper (Coughlin & Begelman 2015, hereafter Paper II) we present an alternate model, the free-streaming jet solution, and compare it to the two-stream solution presented here.

As mentioned in the Introduction, radiative forces behave like an effective viscosity in the presence of shear, when the change in fluid quantities across the mean free path of a photon is small. The precise means by which this viscous coupling manifests itself can be determined by analyzing the Boltzmann equation.

Coughlin & Begelman (2014b) used the general relativistic Boltzmann equation for Thomson scattering to discern the effects of radiation viscosity in the presence of both relativistic velocities and gravitational fields (i.e., in accelerating reference frames). Instead of reproducing their work here, we will simply quote the equations of radiation hydrodynamics for a cold gas (gas pressure much less than both the gas rest mass density and the radiation pressure) that result from their analysis (see their Equation (49)):

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left[\left\{{\rho }^{\prime }+\displaystyle \frac{4}{3}e^{\prime} \left(1-\displaystyle \frac{10}{9}\displaystyle \frac{1}{{\rho }^{\prime }\kappa }{{\rm{\nabla }}}_{\alpha }{U}^{\alpha }\right)\right\}{U}^{\mu }{U}^{\nu }\right]\\ &&-\displaystyle \frac{8}{27}{{\rm{\nabla }}}_{\mu }\left[\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }{{\rm{\Pi }}}^{\mu \sigma }{{\rm{\Pi }}}^{\nu \beta }\left({{\rm{\nabla }}}_{\sigma }{U}_{\beta }+{{\rm{\nabla }}}_{\beta }{U}_{\sigma }+{g}_{\beta \sigma }{{\rm{\nabla }}}_{\alpha }{U}^{\alpha }\right)\right]\\ &&+\displaystyle \frac{1}{3}{g}^{\mu \nu }{\partial }_{\mu }e^{\prime} -\displaystyle \frac{1}{3}{{\rm{\nabla }}}_{\mu }\left[\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }\left({{\rm{\Pi }}}^{\mu \sigma }{U}^{\nu }+{{\rm{\Pi }}}^{\nu \sigma }{U}^{\mu }\right)\right.\\ &&\left.\times \left(4{U}^{\beta }{{\rm{\nabla }}}_{\beta }{U}_{\sigma }+{\partial }_{\sigma }\mathrm{ln}e^{\prime} \right)\right]=0.\end{eqnarray} \tag{ 1 }$$

Here the speed of light has been set to one, Greek indices range from 0 to 3, $\rho ^{\prime}$ is the fluid rest-frame density of scatterers, e' is the fluid rest-frame radiation energy density, κ is the scattering opacity (in units of cm² g⁻¹), ${g}_{\mu \nu }$ is the metric of the spacetime, ${{\rm{\nabla }}}_{\mu }$ is the covariant derivative, ${U}^{\mu }$ is the four-velocity of the flow, and ${{\rm{\Pi }}}^{\mu \nu }={U}^{\mu }{U}^{\nu }+{g}^{\mu \nu }$ is the projection tensor. The Einstein summation convention has been adopted here, meaning that repeated upper and lower indices imply summation. This equation shows that the coefficient of dynamic viscosity, η, for an optically thick, radiation-dominated gas is

$$\begin{eqnarray}&&\eta =\displaystyle \frac{8}{27}\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa },\end{eqnarray} \tag{ 2 }$$

which agrees with previous findings (Blandford et al. 1985; Loeb & Laor 1992); note that this specific coefficient is only for the case when the gas and radiation interact through Thomson scattering. The gas energy equation, which will also be useful for us, can be obtained by contracting Equation (1) with ${U}_{\nu }$, which we can show becomes (see Equation (50) of Coughlin & Begelman 2014b)

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }(e^{\prime} {U}^{\mu })+\displaystyle \frac{1}{3}e^{\prime} {{\rm{\nabla }}}_{\mu }{U}^{\mu }=\displaystyle \frac{4}{3}\displaystyle \frac{10}{9}{{\rm{\nabla }}}_{\mu }\left[\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }{U}^{\mu }{{\rm{\nabla }}}_{\alpha }{U}^{\alpha }\right]\\ &&+\displaystyle \frac{8}{27}\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }\left({{\rm{\nabla }}}_{\sigma }{U}_{\beta }+{{\rm{\nabla }}}_{\beta }{U}_{\sigma }+{g}_{\sigma \beta }{{\rm{\nabla }}}_{\alpha }{U}^{\alpha }\right){{\rm{\Pi }}}^{\mu \sigma }{{\rm{\nabla }}}_{\mu }{U}^{\beta }\\ &&+\displaystyle \frac{1}{3}{{\rm{\Pi }}}^{\mu \sigma }{{\rm{\nabla }}}_{\mu }\left[\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }\left(4{U}^{\beta }{{\rm{\nabla }}}_{\beta }{U}_{\sigma }+{\partial }_{\sigma }\mathrm{ln}e^{\prime} \right)\right]+\displaystyle \frac{1}{3}\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }\\ &&\times \left(4{U}^{\beta }{{\rm{\nabla }}}_{\beta }{U}_{\sigma }+{\partial }_{\sigma }\mathrm{ln}e^{\prime} \right)\left(2{U}^{\mu }{{\rm{\nabla }}}_{\mu }{U}^{\sigma }+{U}^{\sigma }{{\rm{\nabla }}}_{\mu }{U}^{\mu }\right).\end{eqnarray} \tag{ 3 }$$

To close the system, we require that the normalization of the four-velocity be upheld and that particle flux be conserved:

$$\begin{eqnarray}&&{U}_{\mu }{U}^{\mu }=-1,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left[\rho ^{\prime} {U}^{\mu }\right]=0.\end{eqnarray} \tag{ 5 }$$

Equations (1) and (3)–(5) constitute six linearly independent equations for the six unknowns ${U}^{\mu }$, e', and $\rho ^{\prime}$.

In addition to the energy density of the radiation, e', one can also calculate the number density of photons by requiring that the number flux, ${F}^{\mu }$, be conserved. One can show (Coughlin & Begelman 2014b) that the equation ${{\rm{\nabla }}}_{\mu }{F}^{\mu }=0$ becomes, in the viscous limit

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left[N^{\prime} {U}^{\mu }\right]={{\rm{\nabla }}}_{\mu }\left[\displaystyle \frac{1}{{\rho }^{\prime }\kappa }\left(\displaystyle \frac{10}{9}N^{\prime} {U}^{\mu }{{\rm{\nabla }}}_{\sigma }{U}^{\sigma }\right.\right.\\ &&\left.\left.+N^{\prime} {U}^{\alpha }{{\rm{\nabla }}}_{\alpha }{U}^{\mu }+\displaystyle \frac{1}{3}{{\rm{\Pi }}}^{\mu \sigma }{{\rm{\nabla }}}_{\sigma }N^{\prime} \right)\right],\end{eqnarray} \tag{ 6 }$$

where N' is the rest-frame number density of photons. Once we solve the equations of radiation hydrodynamics for the four-velocity of the fluid and the mass density of scatterers, we can solve Equation (6) for the number flux of photons.

The goal of the next two sections is to apply Equations (1) and (3)–(6) to the boundary layers established between fast-moving jets and their environments. For a more thorough discussion of the nature of the equations of radiation hydrodynamics in the viscous limit, we refer the reader to Coughlin & Begelman (2014b).

Arav & Begelman (1992) considered the Blasius (1908) boundary layer problem, wherein one analyzes the dynamics of viscous flow over a semi-infinite, rigid plate, with the viscous effects attributed to radiation. They showed, however, that the requirement that both velocity components vanish on the lower plate, the no-slip condition, results in a divergent boundary layer thickness. The authors then examined the more physical scenario of two interacting fluids, one moving at some asymptotic velocity and the other asymptotically stationary, known as the two-stream approximation. In this case the no-slip condition no longer applies, and they were able to show that the boundary layer thickness remains finite.

The treatment of Arav & Begelman (1992) was non-relativistic, meaning that their results have limited utility when one considers the boundary layers established between GRBs and super-Eddington TDEs and their ambient media. The enthalpy of the radiation, which is not ignorable in radiation-dominated flows, was also not included in their momentum equation. Here we extend their analysis to the relativistic regime and we include the radiation enthalpy.

3.1. Basic Setup

3.2. Boundary Layer Equations

The complexity of Equations (1) and (3)–(5) can be reduced by introducing the boundary layer thickness $\delta y$ such that $\delta \sim \delta y/\delta z$ is a small parameter when $\delta z$ is chosen to be a typical length scale in z. By keeping terms only to lowest order in δ, we will recover a set of reduced boundary layer equations.

To determine how δ depends on asymptotic fluid quantities (it is the reciprocal of the square root of the Reynolds number in the classical Blasius boundary layer), we compare the lowest-order terms in δ in the gas energy equation to the inviscid terms. By equating these terms we are requiring that the viscous heating of the radiation contribute a non-negligible increase in the entropy of the fluid, but because the inviscid terms are proportional to the divergence of the four-velocity (see Equation (10)), this equality can only be true when the gas is compressible, i.e., when the flow velocity is supersonic. When the flow becomes very subsonic, the energy equation can be ignored as the fluid is essentially incompressible. Making the substitutions $\partial /\partial y\sim 1/\delta y$, $\partial /\partial z\sim 1/z$, $v\sim {v}_{j}$ and ${\rho }^{\prime }\sim {\rho }_{0}^{\prime }$ in Equation (3), where ${\rho }_{0}^{\prime }={\rho }_{j}^{\prime }$ if we are in the jet ($y\gt 0$) or ${\rho }_{a}^{\prime }$ if we are in the ambient medium ($y\lt 0$), we find that the boundary layer thickness scales as

$$\begin{eqnarray}&&{\delta }^{2}\sim \displaystyle \frac{1}{{\rho }_{0}^{\prime }\;\kappa \;z\;{{\rm{\Gamma }}}_{j}\;{v}_{j}},\end{eqnarray} \tag{ 7 }$$

where ${{\rm{\Gamma }}}_{j}={\left(1-{v}_{j}^{2}\right)}^{-1/2}$ is the Lorentz factor of the jet. The boundary layer thickness therefore scales roughly as $1/\sqrt{\tau }$, where $\tau \sim {\rho }_{0}^{\prime }\kappa z$ is the fluid-frame optical depth along z.

We can now use our expression for δ (7) to keep only lowest-order terms in Equations (1) and (3)–(5). The resulting $\nu =y$, $\nu =z$, gas-energy and continuity boundary layer equations are, respectively,

$$\begin{eqnarray}&&\displaystyle \frac{\partial e^{\prime} }{\partial y}=0,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {x}^{\mu }}\left[\left({\rho }^{\prime }+\displaystyle \frac{4}{3}e^{\prime} \right){U}^{\mu }\;{\rm{\Gamma }}{v}_{z}\right]+\displaystyle \frac{1}{3}\displaystyle \frac{{de}^{\prime} }{{dz}}\\ &&-\displaystyle \frac{8}{27}\displaystyle \frac{\partial }{\partial y}\left[\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }\displaystyle \frac{\partial }{\partial y}\left[{\rm{\Gamma }}{v}_{z}\right]\right]=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {x}^{\mu }}\left[e^{\prime} {U}^{\mu }\right]+\displaystyle \frac{1}{3}e^{\prime} \displaystyle \frac{\partial {U}^{\mu }}{\partial {x}^{\mu }}=\displaystyle \frac{8}{27}\displaystyle \frac{e^{\prime} }{{\rho }^{\prime }\kappa }\displaystyle \frac{\partial {U}_{\mu }}{\partial y}\displaystyle \frac{\partial {U}^{\mu }}{\partial y},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial z}\left[{\rho }^{\prime }{\rm{\Gamma }}{v}_{z}\right]+\displaystyle \frac{\partial }{\partial y}\left[{\rho }^{\prime }{\rm{\Gamma }}{v}_{y}\right]=0.\end{eqnarray} \tag{ 11 }$$

The first of these demonstrates, as in the non-relativistic limit, that the pressure is constant across the boundary layer. For the remainder of this section we will assume that $e^{\prime} (z)={e}_{0}^{\prime }$, i.e., that the radiation energy density is independent of z.

3.3. Self-similar Approximation

The solution to the continuity Equation (11) can be obtained by introducing the stream function ψ through the definitions

$$\begin{eqnarray}&&{\rho }^{\prime }\kappa {\rm{\Gamma }}{v}_{z}=\displaystyle \frac{\partial \psi }{\partial y},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\rho }^{\prime }\kappa {\rm{\Gamma }}{v}_{y}=-\displaystyle \frac{\partial \psi }{\partial z},\end{eqnarray} \tag{ 13 }$$

which manifestly solves the continuity equation; we introduced a factor of κ, the opacity, to ensure that ψ is dimensionless. These relations also demonstrate that ${v}_{y}\sim \delta {v}_{z}$, which is what we expect: the velocity in the direction perpendicular to the majority of the motion is reduced by a factor of δ. Note, however, that the definition of the stream function now involves the density, which is not constant in this analysis.

As is done in the standard Blasius treatment, we assume that the stream function varies self-similarly as

$$\begin{eqnarray}&&\psi =\sqrt{\displaystyle \frac{8}{27}{\rho }_{0}^{\prime }\;\kappa \;{{\rm{\Gamma }}}_{j}\;{v}_{j}}\;{z}^{1/2}f(\alpha ),\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&\alpha =y/\delta y=\displaystyle \frac{y}{\sqrt{z}}\sqrt{\displaystyle \frac{27}{8}{\rho }_{0}^{\prime }\;\kappa \;{{\rm{\Gamma }}}_{j}\;{v}_{j}}\end{eqnarray} \tag{ 15 }$$

is our self-similar variable (we used Equation (7) for the boundary layer thickness) and f is a function to be determined from Equations (9) and (10). We will also assume that the density varies self-similarly as

$$\begin{eqnarray}&&{\rho }^{\prime }={\rho }_{0}^{\prime }\;g(\alpha ),\end{eqnarray} \tag{ 16 }$$

where g is a second function.

One can use these definitions in Equations (9) and (10) to derive a set of coupled, nonlinear, ordinary differential equations for f and g. It is mathematically convenient, however, to define a new self-similar variable ξ by

$$\begin{eqnarray}&&\xi ={\displaystyle \int }_{0}^{\alpha }g(\tilde{\alpha })d\tilde{\alpha },\end{eqnarray} \tag{ 17 }$$

where $\tilde{\alpha }$ is a dummy variable of integration, and write the functions f and g in terms of this variable. This approach is similar to the one taken by Arav & Begelman (1992).

With this parametrization, the velocities are

$$\begin{eqnarray}&&{\rm{\Gamma }}{v}_{z}={{\rm{\Gamma }}}_{j}{v}_{j}{f}_{\xi },\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\rho }^{\prime }{\rm{\Gamma }}{v}_{y}=\displaystyle \frac{1}{2\sqrt{z}}\sqrt{\displaystyle \frac{8}{27}\displaystyle \frac{{\rho }_{0}^{\prime }\;{{\rm{\Gamma }}}_{j}\;{v}_{j}}{\kappa }}\left(\alpha {f}_{\xi }g-f\right),\end{eqnarray} \tag{ 19 }$$

where a function with a subscript ξ denotes the derivative of that function with respect to ξ, i.e., ${f}_{\xi }={df}/d\xi$, ${f}_{\xi \xi }\;=\;{d}^{2}f/d{\xi }^{2}$, etc. Substituting these relations and the self-similar scaling for $\rho ^{\prime}$, Equation (16), into Equations (9) and (10), we find the following self-similar equations for f and g:

$$\begin{eqnarray}&&-\displaystyle \frac{1}{2}\left(g+\displaystyle \frac{4}{3}\mu \right)f\;{f}_{\xi \xi }+\mu \;{{\rm{\Gamma }}}_{j}^{2}{v}_{j}^{2}\displaystyle \frac{g\;{f}_{\xi }{({f}_{\xi \xi })}^{2}}{1+{v}_{j}^{2}{{\rm{\Gamma }}}_{j}^{2}{({f}_{\xi })}^{2}}=\mu \;g\;{f}_{\xi \xi \xi },\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{g}_{\xi }f=\displaystyle \frac{3}{2}{{\rm{\Gamma }}}_{j}^{2}{v}_{j}^{2}\displaystyle \frac{{g}^{2}{({f}_{\xi \xi })}^{2}}{1+{{\rm{\Gamma }}}_{j}^{2}{v}_{j}^{2}{({f}_{\xi })}^{2}},\end{eqnarray} \tag{ 21 }$$

where we have defined μ as the ratio ${e}_{0}^{\prime }/{\rho }_{0}^{\prime }$. The term $-2\;\mu \;f\;{f}_{\xi \xi }/3$ and the last term on the left-hand side of Equation (20) were absent in the treatment of Arav & Begelman (1992) because they did not include the enthalpy of the radiation. In the non-relativistic, $\mu \ll 1$ limit, Equation (20) reduces to the standard Blasius equation (by rescaling the self-similar variable) and g and f decouple from one another, as was found by Arav & Begelman (1992).

Before proceeding further, recall that the two-stream problem separates the jet and ambient material into two distinct media, the dividing line for our problem chosen to be y = 0, and that the asymptotic densities attained in these two media are ${\rho }_{j}^{\prime }$ and ${\rho }_{a}^{\prime }$ for the jet and ambient material, respectively. Therefore, in Equations (14)–(16), the parameter ${\rho }_{0}^{\prime }$ refers to either ${\rho }_{j}^{\prime }$ or ${\rho }_{a}^{\prime }$ depending on the sign of α and, hence, ξ. The functions f and g are thus piecewise defined about y = 0, as are the self-similar variables α and ξ, with solutions for $\xi \gt 0$ corresponding to jet quantities and those for $\xi \lt 0$ corresponding to the ambient medium. Equations (20) and (21) should therefore be interpreted as a total of four equations for four functions: f defined in the jet, f defined in the ambient medium, g defined in the jet, and g defined in the ambient medium. Note that μ can also differ between the two media, depending on the asymptotic density (but ${e}_{0}^{\prime }$ must be continuous across y = 0 because of Equation (8)).

Keeping in mind this subtlety of Equations (20) and (21), we must additionally determine the boundary conditions satisfied by f and g. The first two conditions satisfied by f can be determined by recalling that $v\to {v}_{j}$ as $y\to \infty$ and $v\to 0$ as $y\to -\infty$. Investigating Equation (18) and noting that ξ scales with y, these requirements translate to

$$\begin{eqnarray}&&{f}_{\xi }(\infty )=1,\quad {f}_{\xi }(-\infty )=0.\end{eqnarray} \tag{ 22 }$$

We also require that the density approach its asymptotic values in the limits of $y\to \pm \infty$. From Equation (16), this gives

$$\begin{eqnarray}&&g(\pm \infty )=1.\end{eqnarray} \tag{ 23 }$$

Now, note that if the jet and ambient materials are to retain their respective identities, then the flow along the surface of contact at y = 0 must remain parallel to that surface. In other words, there must not be any mass flow across the boundary, i.e., this surface is a contact discontinuity, which means that ${v}_{y}(0)=0$. From Equation (19), this shows that f must satisfy

$$\begin{eqnarray}&&f(0)=0.\end{eqnarray} \tag{ 24 }$$

The other two boundary conditions can be determined by requiring that the normal and transverse components of the energy–momentum tensor, the divergence of which we took to obtain Equation (1), be continuous across the point of contact y = 0. We can show that these restrictions demand that ${f}_{\xi }$ and ${f}_{\xi \xi }$ be continuous across y = 0, which closes the system.

We can also calculate the comoving-frame number density of photons N' by solving Equation (6). Keeping only terms to lowest order in the boundary layer thickness, the equation of conservation of photon number becomes

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{F}^{\mu }={{\rm{\nabla }}}_{\mu }\left(N^{\prime} {U}^{\mu }\right)-\displaystyle \frac{1}{3}\displaystyle \frac{\partial }{\partial y}\left(\displaystyle \frac{1}{{\rho }^{\prime }\kappa }\displaystyle \frac{\partial N^{\prime} }{\partial y}\right)=0.\end{eqnarray} \tag{ 25 }$$

As we did for the number density of scatterers, we will assume that the photon number density varies self-similarly across the boundary layer as

$$\begin{eqnarray}&&N^{\prime} ={N}_{0}^{\prime }\;h(\xi ),\end{eqnarray} \tag{ 26 }$$

where ${N}_{0}^{\prime }$ is the asymptotic number density of the jet or the ambient medium. With this form for the number density, we find that the equation for h becomes

$$\begin{eqnarray}&&\displaystyle \frac{3}{2}{{\rm{\Gamma }}}_{j}^{2}{v}_{j}^{2}\displaystyle \frac{g{({f}_{\xi \xi })}^{2}}{1+{{\rm{\Gamma }}}_{j}^{2}{v}_{j}^{2}{({f}_{\xi })}^{2}}h-f\;{h}_{\xi }=\displaystyle \frac{9}{4}g\;{h}_{\xi \xi }.\end{eqnarray} \tag{ 27 }$$

As for the functions f and g, h is really piecewise defined across the boundary y = 0, and so Equation (27) is really two equations—one for the photon number density in the jet, and another for that in the ambient medium. The requirement that the number density of photons asymptotically approach its jet and ambient values gives

$$\begin{eqnarray}&&h(\pm \infty )=1.\end{eqnarray} \tag{ 28 }$$

The y-component of the flux must also be continuous across the boundary, which, when written out, shows that the derivative of h must be continuous across $\xi =0$, which yields the final two boundary conditions.

Finally, the solutions to Equations (20), (21), and (27) will be in terms of the parameter ξ, which is itself a function of g (Equation (17)). We would like the solutions to be in terms of the parameter α, which is directly related to the physical coordinates y and z (Equation (15)). The transformation can be achieved by differentiating Equation (17), rearranging, and integrating to yield

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\xi }\displaystyle \frac{d\tilde{\xi }}{g(\tilde{\xi })}=\alpha ,\end{eqnarray} \tag{ 29 }$$

where $\tilde{\xi }$ is an integration variable and we set the integration constant to zero because we demand $\xi (\alpha =0)=0$. After solving Equations (20) and (21) for $g(\xi )$, we can numerically integrate and solve Equation (29) for ξ in terms of α.

3.4. Solutions

Equations (20), (21) and (27), together with the boundary conditions on f, g, and h, govern the behavior of the velocity, number density of scatterers, and number density of photons throughout the two-stream boundary layer. Here we present and analyze the solutions to those equations as we vary the asymptotic jet Lorentz factor ${{\rm{\Gamma }}}_{j}$ and the quantity $\mu ={e}_{0}^{\prime }/{\rho }_{0}^{\prime }$.

One caveat with the definition of α is that it depends on ${{\rm{\Gamma }}}_{j}$ via Equation (15). Therefore, if we plot solutions with different ${{\rm{\Gamma }}}_{j}$, we must be careful to incorporate this dependence so that the range of physical space we consider for each solution is the same. Because of this fact, in the following figures we will plot our solutions as functions of the variable

$$\begin{eqnarray}&&\tilde{\alpha }=\displaystyle \frac{\alpha }{\sqrt{{{\rm{\Gamma }}}_{j}{v}_{j}}}=\displaystyle \frac{y}{\sqrt{z}}\sqrt{\displaystyle \frac{27}{8}{\rho }_{0}^{\prime }\kappa }.\end{eqnarray} \tag{ 30 }$$

Figure 1 shows the solution for the normalized z-component of the four-velocity (${f}_{\xi }$) as a function of the self-similar variable $\tilde{\alpha }\propto y/\sqrt{z}$ for a number of jet Lorentz factors and $\mu =1$. We see that the width of the boundary layer is nearly unchanged as we vary the Lorentz factor of the asymptotic jet. The value of ${f}_{\xi }$ at the contact discontinuity is lower for larger ${{\rm{\Gamma }}}_{j}$, resulting in a greater shear ($\sim {f}_{\xi \xi }$) as one proceeds into the jet. The flattening of the velocity around $\tilde{\alpha }=0$ arises from the behavior of the density around this region and the function $\xi (\alpha )$ determined therefrom.

Zoom In Zoom Out Reset image size

Figure 2 illustrates the manner in which the density varies over the boundary layer for the same parameters as those chosen for Figure 1. Consistent with the findings of Arav & Begelman (1992), we find that the transition from the jet to the ambient medium carves out a region of low density material. This behavior can be understood by noting that the shear in the flow causes viscous heating of the fluid, which results in an increase in the specific entropy s'. Since the specific entropy scales as $s^{\prime} \propto e^{\prime} /{({\rho }^{\prime })}^{4/3}$ for a radiation-dominated gas, an increase in the entropy at constant pressure corresponds to an decrease in the density of scatterers. As was found by Arav & Begelman (1992), the density equals zero at the origin, which can be gleaned from the asymptotic behavior of Equation (21): for $\xi \ll 1$, we can let ${f}_{\xi }\sim {f}_{\xi \xi }\sim \xi$, and solving the resultant approximate differential equation shows that $g\propto -1/\mathrm{ln}(\xi )$.

Zoom In Zoom Out Reset image size

In Figure 3 we plot the solution to Equation (27), the normalized, rest-frame number density of photons, for the same set of Lorentz factors and $\mu =1$. By comparing this with Figure 2, we see that the photon number density roughly tracks that of the scatterers. However, the photon number density remains below its asymptotic value significantly farther into the ambient medium than does the particle density. Because the energy density of the radiation is constant across the boundary layer, the decrease in photon density corresponds to a higher average energy per photon increases within the boundary layer—a clear manifestation of viscous heating, as noted by Arav & Begelman (1992).

Zoom In Zoom Out Reset image size

Figures 4–6 show the z-component of the four-velocity, the number density of scatterers, and the photon number density, respectively, for ${{\rm{\Gamma }}}_{j}=25$ and $\mu =0.1$, 1, and 10. The overall qualitative behavior of the solutions is similar to that depicted in Figures 1–3. More specifically, however, we find that values of μ greater than unity cause the boundary layer thickness to decrease, but not appreciably. On the contrary, a value of μ only marginally less than one results in a significant widening of the boundary layer. This dependence is ultimately related to the compressibility of the fluid and the relation between that compressibility and the sound speed of a radiation-dominated gas (see discussion below).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The plots of the previous subsection demonstrate how radiation viscosity affects the boundary between a fast-moving flow and its surroundings under the two-stream approximation. Our analysis generalizes the treatment of Arav & Begelman (1992) by permitting the jet velocity to be relativistic and by incorporating the enthalpy of the radiation in the momentum equation. Our results are similar to those found by Arav & Begelman (1992; compare, e.g., their Figure 3 to our Figures 1 and 2); there are, however, a few notable differences.

For one, our boundary layer thickness, given by Equation (7), differs from that of Arav & Begelman (1992), who found ${\delta }^{2}\sim e^{\prime} /({\rho }_{0}^{\prime }{}^{2}\kappa \;{v}_{j}z)$. In addition to the Lorentz factor contained in ours, their boundary layer thickness has an additional factor of $e^{\prime} /{\rho }_{0}^{\prime }$, meaning that, in the non-relativistic limit, our results do not agree. This discrepancy arises from the fact that, while we compared the lowest-order terms in δ to the inviscid terms in the gas energy equation to obtain our boundary layer thickness, they compared the viscous term to the inertial term—the one proportional to $\rho ^{\prime}$—in the z-component of the momentum equation to obtain theirs. Because they ignored the enthalpy of the radiation, the inertial term was the only inviscid term present in the momentum equation, making it the only one available to balance the viscous part. However, if one does not ignore the radiative contribution to the momentum equation, one can now equate the viscous terms to either the term proportional to the mass density or the one proportional to the enthalpy. Because the density equals zero at $\alpha =0$, there is always some location at which the enthalpy term exceeds the inertial term, making it more reasonable to equate the former to the viscous contribution than the latter. If one follows this route, one recovers our ordering for the boundary layer thickness.

The second difference is that their solutions depend on the square of the Mach number, ${M}^{2}\equiv {v}^{2}/{c}_{{\rm{s}}}^{2}$, where ${c}_{{\rm{s}}}^{2}\;=4e^{\prime} /(9{\rho }_{0}^{\prime })$ is the non-relativistic sound speed. Our solutions, on the other hand, depend on both the jet velocity v_j and the ratio $\mu =e^{\prime} /{\rho }_{0}^{\prime }$, which is proportional to the non-relativistic sound speed. One reason for this difference arises from the discrepancy between our definitions of the boundary layer thickness. Another is due to the fact that our solutions included the enthalpy of the radiation; had Arav & Begelman (1992) included this term, factors of the sound speed would have arisen in their z-momentum equation. Finally, our solutions also depend on the value of the jet velocity—not just the ratio of the jet speed to the sound speed—because we included all relativistic effects, meaning that the speed of light now plays a role in determining the evolution of the system.

Our solutions show that, for fixed μ, the thickness of the boundary layer is approximately independent of the jet Lorentz factor, which is due to the competition between the scaling of the fiducial boundary layer thickness $\delta \propto 1/\sqrt{{{\rm{\Gamma }}}_{j}}$ (see Equation (7)) and the viscous heating. Specifically, a larger ${{\rm{\Gamma }}}_{j}$ results in a smaller δ and a greater shear; this shear increases the specific entropy $s^{\prime} \propto e^{\prime} /{({\rho }^{\prime })}^{4/3}$ and, since the pressure is a constant, this increase in the entropy implies a decrease in the density of scatterers which widens the boundary layer.

Figures 4–6 show how the solutions for the velocity, density of scatterers, and density of photons vary for a fixed jet Lorentz factor (we chose ${{\rm{\Gamma }}}_{j}=25$) but for a variable μ. Increasing μ relative to $\mu =1$ tends to slightly decrease the boundary layer thickness, while decreasing the value of μ relative to $\mu =1$ dramatically increases the thickness. This behavior arises from the fact that the viscous heating, which decreases the density of scatterers and widens the boundary layer, is effective when the gas is compressible. The compressibility of the fluid, however, is only important when the flow velocity is supersonic, and we can show that the sound speed of the gas is (Coughlin & Begelman 2014b)

$$\begin{eqnarray}&&{c}_{{\rm{s}}}=\displaystyle \frac{2}{3}\sqrt{\displaystyle \frac{\mu }{1+4\mu /3}},\end{eqnarray} \tag{ 31 }$$

which, noting that $e^{\prime} =3p^{\prime}$, where p' is the radiation pressure, correctly reduces to ${c}_{{\rm{s}}}\simeq \sqrt{4p^{\prime} /(3{\rho }^{\prime })}$ in the limit $\mu \ll 1$ and ${c}_{{\rm{s}}}\simeq 1/\sqrt{3}$ in the limit $\mu \gg 1$. Thus, when $\mu \ll 1$, the location of the sonic point extends farther into the ambient medium, resulting in a widening of the boundary layer. On the other hand, when $\mu \gg 1$, the sonic point approaches the jet, but only slightly due to the fact that the sound speed approaches a constant as $\mu \to \infty$. In fact, based on this reasoning, we expect our solutions to be independent of μ in the large μ limit, which is indeed reflected in Equation (20).

The densities of scatterers and photons both decrease within the boundary layer. This behavior has two interesting consequences. First, the lower number density of scatterers means that the optical depth is lower in the boundary layer. We are therefore able to see farther into the medium along lines of sight that probe regions of high shear. Second, because the energy density of the radiation stays unchanged as we traverse the media, the average energy per photon increases, resulting in a harder spectrum within the boundary layer.

Equation (8) shows that the radiation energy density, and hence the pressure, is constant across the boundary layer, which is ultimately a statement of the causal connectedness of the jet. This means, equivalently, that the boundary layer thickness $\delta y$ can be traversed by a sound wave in less time than it takes the jet material to cross the distance $\delta z$. Since the transverse sound speed is related to the isotropic sound speed by ${c}_{s\perp }={c}_{{\rm{s}}}/{{\rm{\Gamma }}}_{j}$, we find that the boundary layer thickness δ must satisfy $\delta \lesssim {c}_{{\rm{s}}}/{{\rm{\Gamma }}}_{j}$. When this inequality is no longer satisfied, Equation (8) no longer holds, and we must include more terms in all of the boundary layer equations that account for changes in the pressure.

The solution for g equals zero at $\alpha =0$ which, as we mentioned, can be determined by inspecting the asymptotic ($\xi \to 0$) limit of Equation (21). This feature was also found by Arav & Begelman (1992), and can be understood physically by noting that, when the jet and the ambient medium initially interacted, the surface separating them served as a discontinuity in velocity, resulting in an infinite shear and entropy generation. Although the solutions presented here illustrate the time-steady state of the system after radiation viscosity has smoothed out the discontinuity, the infinite entropy along the contact discontinuity persists and drives the density to zero. Mathematically, this behavior is ultimately due to the fact that $f(0)=0$, which itself came from the requirement that there be no mixing at the surface separating the jet and the external environment. This boundary condition is necessary to ensure that the two media retain their respective identities and underlies the two-stream assumption, and it allows us to prescribe different asymptotic properties of those media.

We applied the relativistic equations of radiation hydrodynamics in the viscous limit to the two-stream boundary layer, expanding on the past work of Arav & Begelman (1992). These equations are applicable as long as changes in the fluid quantities are small over the mean free path of a photon.

An interesting feature of the solutions presented here is the depression in the number density of scatterers within the boundary layer separating the jet and the ambient medium, which is consistent with the findings of Arav & Begelman (1992). We also showed that the number density of photons N' roughly tracks the density of scatterers, reaching a minimum toward, but not at, the contact discontinuity. Therefore, observers viewing a relativistic outflow with lines of sight that probe regions of high shear see farther into the outflow and they also see a higher energy per photon, and hence a harder spectrum.

Our solutions show that the physical boundary layer thickness does not depend strongly on the jet Lorentz factor ${{\rm{\Gamma }}}_{j}$, which results from a competition between viscous heating and the scaling of the fiducial boundary layer thickness $\delta \sim 1/\sqrt{{{\rm{\Gamma }}}_{j}}$. The dependence of our solutions on μ arises from the fact that the change in entropy of the flow is related to its compressibility, which is most important where the outflow velocity is supersonic. Since the sound speed scales as ${c}_{{\rm{s}}}\sim \sqrt{\mu }$ when $\mu \lesssim 1$, the point at which the outflow becomes subsonic extends farther into the ambient medium when μ is small, resulting in a widening of the boundary layer. When μ becomes larger than one, however, the sound speed does not increase much, asymptotically approaching $1/\sqrt{3}$, meaning that the sonic point only slightly approaches the jet, yielding a marginal decrease in the boundary layer thickness.

A number of assumptions about the nature of the jet and its surroundings are built in to the two-stream solutions presented here. Specifically, we adopted a two-dimensional, plane-parallel geometry for the flow and its surroundings, and we demanded that there be no pressure gradient ($e^{\prime} (z)={e}_{0}^{\prime }$) in the ambient medium. While these assumptions greatly enhanced the tractability of the problem, they somewhat hinder the astrophysical relevance of the solutions, as no systems likely conform precisely to these restrictions. However, in local regions of an outflow, where the radius of curvature is large and the pressure can be considered relatively constant, these solutions may be actualized. Therefore, while the physics of an entire jet-disk system may be poorly modeled by the two-stream scenario, local patches of the outflow, where the geometry can be treated as locally flat and the pressure gradient is small, are likely well-described by the solutions presented here.

The two-stream solutions could be applied to relativistic, radiation dominated jets, such as those that appear during super-Eddington TDEs, the event Swift J1644+57 being one such case. The event Swift J2058+05, observed shortly after Swift J1644, is another source that has been interpreted as a jetted TDE (Cenko et al. 2012). A comparison between the models presented here and these sources could provide valuable information concerning their progenitors and the properties of the jets themselves. The application of these models to long GRBs could likewise prove fruitful, potentially yielding, for example, information concerning the Lorentz factor of the jet, the stellar progenitor, and the direction of the line of sight between the observer and the source.

These models may also be applicable to jetted X-ray binary systems, or "microquasars" (Mirabel & Rodríguez 1999; Fender et al. 2004, 2009). For example, Arav & Begelman (1993) applied their non-relativistic, radiation-viscous solution (Arav & Begelman 1992) to the source SS 433 (Fabrika 2004; Begelman et al. 2006a). Since the jets of SS 433 are mildly relativistic, their speeds being ${v}_{j}\simeq 0.26\;c$ (Margon & Anderson 1989), including the relativistic terms arising from the treatment presented here may place new constraints on the properties of those jets and the surrounding envelope.

Finally, quasi-stars—protogalactic gas clouds supported by an accreting black hole—may also contain jets (Begelman et al. 2006b, 2008; Czerny et al. 2012). Since the power radiated by the black hole at the center of a quasi-star supports the overlying gaseous envelope, the mass of the envelope greatly exceeding that of the hole, the black hole accretes supercritically by several orders of magnitude. The jets launched from these systems are therefore radiation-dominated, and as they propagate through the quasi-star envelope, radiation-viscous effects likely dominate the interaction between the two media. Applying the solutions presented here to these systems could then yield information about the properties of these jets and the role they may have played during the epoch of galaxy formation.

One drawback to these models and their physical application, however, is that they require that the jet and ambient medium be separated by a contact discontinuity, which results in the non-physical vanishing of the density of scatterers at the interface. Furthermore, this prevents the jet from entraining ambient material; while this is not particularly problematic for the two-stream problem, in which the jet is considered infinite in extent, realistic jets have a finite width and total momentum, meaning that the entrainment of ambient material will cause a decrease in the outflow velocity with z that cannot be captured with the two-stream treatment.

In Paper II, we investigate a different type of viscous boundary layer—the free-streaming jet model. This model treats the entire system, jet and ambient medium, as a single fluid, which removes the need for a contact discontinuity and allows the density to remain non-zero throughout the boundary layer. We also show that the entrainment of ambient material causes an overall slowing of the outflow.

This work was supported in part by NASA Astrophysics Theory Program grant NNX14AB37G, NSF grant AST-1411879, and NASA's Fermi Guest Investigator Program.