THE RADIATION HYDRODYNAMICS OF RELATIVISTIC SHEAR FLOWS

Radiation interacts dynamically with matter in many astrophysical systems. Often, the density of the plasma in these systems is high enough that the photon mean free path is very small compared to the scales characterizing the fluid. In these optically thick situations, the radiation field is very nearly isotropic in the rest frame of the fluid, and the dynamical coupling is determined by local quantities. On the other hand, there are situations in which the mean free path of a photon greatly exceeds the physical length scale of the gas that scatters those photons. For these optically thin scenarios, the properties of the radiation field, and its dynamical effects, are largely determined nonlocally.

In these optically thick and thin limits, the manner in which radiation influences the gas is well-understood. For the former, the radiation field generates an energy density and pressure (in addition to the gas pressure) that are characterized by isotropy in the local comoving frame of the matter and an equation of state that reflects the relativistic nature of a photon gas (so that the pressure is one third the energy density; Mihalas & Mihalas 1984). For the latter, the radiation field does not respond locally to the properties of the gas, but can be considered to be imposed externally. The angular distribution of radiation coupled with the state of motion of the gas then generates a "radiation drag" that can cause particles in the solar system to gradually spiral into the Sun (the Poynting–Robertson effect; Robertson 1937) or an acceleration that can propel material away from active galactic nuclei (the Compton rocket effect; O'Dell 1981; Phinney 1982).

An interesting question then arises when one considers how radiation couples to a gas that is between the optically thick and thin limits. Some progress has been made on this question from the optically thick standpoint by assuming that the radiation field is approximately isotropic in the instantaneous rest frame of the fluid. The anisotropic contribution can then be determined from the Boltzmann equation, yielding the general relativistic equations of radiation hydrodynamics in the viscous limit (Coughlin & Begelman 2014b). These equations are "viscous" in the sense that the shear of the fluid—the change in the fluid velocity over the mean free path of a photon—serves to transfer energy and momentum between the radiation and the scatterers. However, it is not at all apparent how this viscous coupling in the optically thick limit transitions to the optically thin counterpart of radiation drag, where the relevant length scale changes from the (assumed-small) mean free path of the photon to the size of the entire fluid (or longer).

This question is also of practical, not just theoretical, importance, as some astrophysical systems are characterized by regions of only marginal optical depth. Such a region naturally arises when an optically thin jet or wind propagates alongside an optically thick disk or envelope, which occurs in ultraluminous X-ray sources (e.g., Arav & Begelman 1992; Begelman et al. 2006), the collapsar model of long gamma-ray bursts (Woosley 1993; MacFadyen & Woosley 1999), and super-Eddington tidal disruption events (Coughlin & Begelman 2014a), to name a few. The transition between the jet and the envelope is then of marginal optical depth, and one must self-consistently account for both the evolution of the fluid and the radiation field within this transition.

Numerically, most authors have resorted to some form of "closure" to capture the important physical effects that take place in these regions of marginal optical depth (but see Jiang et al. 2014 and Ryan et al. 2015), popular choices being flux-limited diffusion (FLD; Levermore & Pomraning 1981; Ohsuga et al. 2009) and M1 (Levermore 1984; McKinney et al. 2014). However, FLD does not capture any relativistic or velocity-dependent effects, which can be quite large in a region of high velocity and intense shear such as the transition between a fast-moving outflow and a hydrostatic envelope. An advantage of M1 is that it is easily extended to incorporate all of the effects of general and special relativity (Sa̧dowski et al. 2013). However, M1 breaks down when the radiation field is inherently very anisotropic (for example, along the axis of the jet as it propagates alongside the envelope) and it does not reduce correctly to the viscous limit (one must apply an additional, viscous term to the stress tensor in the radiation-rest frame; Coughlin & Begelman 2014b; Sa̧dowski et al. 2014).

In this paper we pursue a simple, but physically motivated approach to analyzing the manner by which scatterers in a fast-moving fluid interact with radiation when the medium is optically thick, thin, or in-between, with the specific application of the intense shear layers generated between relativistic jets and their surroundings in mind. In particular, we assume that when such a shear layer develops, the fluid rapidly assumes a self-similar form so that the fluid looks identical at every comoving point. In Section 2 we show that this assumption of self-similarity yields a very specific, one-parameter family of velocity profiles that is described solely by the amount of shear present in the flow.

We further assume that the photons interacting with a given fluid element were all scattered on a prescribed surface located at an optical depth $\tau \simeq 1$ relative to that fluid element, and that the scattering is isotropic in the rest frame of the scatterer and elastic (i.e., the scattering takes place in the Thomson limit). We calculate the shape of the $\tau =1$ surface in Section 3, and we show how it depends on the amount of shear within the flow.

Section 4 combines the assumption of the self-similarity of the flow and the $\tau =1$ scattering specification, demonstrating that the number density of photons is manifestly conserved, both in time and space, throughout the shear layer (as must be true, since we are only considering photon scattering). We also show that, in order for the radiation energy density to be simultaneously conserved in a scattering event and uniform across the shear layer, both necessary for the consistency of this model, the photon energy density must be an increasing function of time and thereby heating up in response to the shear present within the flow. In Section 5 we solve for this heating rate in the optically thin and thick limits, and we provide an approximate, interpolated solution for the heating rate when the optical depth is marginal. Section 6 presents the radiation energy–momentum tensor, and we derive the functional forms of the shear stress and the pressures of the radiation field, showing that they have relatively simple, approximate, analytic forms that match the viscous and streaming limits. We summarize and discuss the implications of our findings in Section 7.

Consider a two-dimensional, planar fluid where the fluid motion is purely along z and the variation in that motion is purely along y. The three-velocity of this fluid is characterized by ${\boldsymbol{v}}=v(y)\hat{z}$, and the four-velocity is, correspondingly,

$$\begin{eqnarray}&&{U}^{\alpha }=\left(\begin{array}{c}{\rm{\Gamma }}\\ 0\\ 0\\ {\rm{\Gamma }}v\end{array}\right),\end{eqnarray} \tag{ 1 }$$

where ${\rm{\Gamma }}={(1-{v}^{2})}^{-1/2}$ is the Lorentz factor.

The fundamental assumption we will make here is that the fluid appears identical to every comoving observer. This is essentially a statement of the self-similarity of the flow, and we suggest that it may apply in regions of shear that have "lost memory" of their boundary or initial conditions. This assumption then implies that the comoving density of matter satisfies ${\rho }^{\prime }(y)={\rho }^{\prime }$, with ${\rho }^{\prime }$ a constant, and that the velocity field varies as

$$\begin{eqnarray}&&{\rm{\Gamma }}=\mathrm{cosh}(\mu {\tau }_{y}^{\prime }),\quad {\rm{\Gamma }}v=\mathrm{sinh}(\mu {\tau }_{y}^{\prime }),\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{\tau }_{y}^{\prime }={\int }_{0}^{y}{\rho }^{\prime }\;\kappa \;d\tilde{y}={\rho }^{\prime }\kappa \;y\end{eqnarray} \tag{ 3 }$$

is the comoving optical depth along the y-direction (κ is the opacity) and μ is a constant that describes the amount of shear present in the flow. We can see that this expression for v has the required properties by considering a fluid element within the flow that moves with four-velocity ${U}^{\alpha }$ with respect to some frame α; relative to another fluid parcel within the shear layer, the frame of which we will denote by β, that same fluid element is observed to be moving at a velocity ${U}^{\beta }$ that is related to ${U}^{\alpha }$ by

$$\begin{eqnarray}&&{U}^{\beta }={{\rm{\Lambda }}}_{\;\alpha }^{\beta }{U}^{\alpha },\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Lambda }}}_{\;\alpha }^{\beta }$ is the local Lorentz transformation between the α and β frames. We thus have

$$\begin{eqnarray}{{\rm{\Lambda }}}_{\;\alpha }^{\beta }=\left(\begin{array}{cccc}{{\rm{\Gamma }}}_{1} & 0 & 0 & -{{\rm{\Gamma }}}_{1}{v}_{1}\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ -{{\rm{\Gamma }}}_{1}{v}_{1} & 0 & 0 & {{\rm{\Gamma }}}_{1}\end{array}\right),\end{eqnarray} \tag{ 5 }$$

where v₁ is the three-velocity of the β frame with respect to the α frame, and

$$\begin{eqnarray}&&{U}^{\alpha }=\left(\begin{array}{c}{{\rm{\Gamma }}}_{2}\\ 0\\ 0\\ {{\rm{\Gamma }}}_{2}{v}_{2}\end{array}\right),\end{eqnarray} \tag{ 6 }$$

where v₂ is the three-velocity of the fluid parcel as measured in the α frame. If we now use Equation (2) for the velocities, then the velocity of the fluid parcel as measured in the β frame is

$$\begin{eqnarray}{U}^{\beta } & = & \left(\begin{array}{c}\mathrm{cosh}(\mu {\tau }_{1}^{\prime })\mathrm{cosh}(\mu {\tau }_{2}^{\prime })-\mathrm{sinh}(\mu {\tau }_{1}^{\prime })\mathrm{sinh}(\mu {\tau }_{2}^{\prime })\\ 0\\ 0\\ -\mathrm{sinh}(\mu {\tau }_{1}^{\prime })\mathrm{cosh}(\mu {\tau }_{2}^{\prime })+\mathrm{cosh}(\mu {\tau }_{1}^{\prime })\mathrm{sinh}(\mu {\tau }_{2}^{\prime })\end{array}\right)\\ & = & \left(\begin{array}{c}\mathrm{cosh}[\mu ({\tau }_{2}^{\prime }-{\tau }_{1}^{\prime })]\\ 0\\ 0\\ \mathrm{sinh}[\mu ({\tau }_{2}^{\prime }-{\tau }_{1}^{\prime })]\end{array}\right),\end{eqnarray} \tag{ 7 }$$

where ${\tau }_{1}^{\prime }$ is the optical depth of the β frame and ${\tau }_{2}^{\prime }$ is the optical depth of the moving gas parcel, both with respect to the α frame. However, it is apparent that

$$\begin{eqnarray}{\tau }_{2}^{\prime }-{\tau }_{1}^{\prime } & = & {\displaystyle \int }_{0}^{{y}_{2}}{\rho }^{\prime }\;\kappa \;d\tilde{y}-{\displaystyle \int }_{0}^{{y}_{1}}{\rho }^{\prime }\;\kappa \;d\tilde{y}\\ & = & {\displaystyle \int }_{{y}_{1}}^{{y}_{2}}{\rho }^{\prime }\;\kappa \;d\tilde{y},\end{eqnarray} \tag{ 8 }$$

where y₁ and y₂ are the positions of the β frame and the gas parcel, respectively, as measured by the α frame. This expression is just the optical depth to the gas parcel as measured by the β frame, and we thus see that the velocity field as measured in the β frame is identical to the velocity field measured in the α frame—precisely the attribute we require for the flow.

Interestingly, this form for the velocity field, Equation (2), also describes the ultrarelativistic limit of optically thick jet propagation mediated by radiation viscosity. We refer the reader to the Appendix for a demonstration of this result.

According to a given fluid element, which we will consider to be the origin, neighboring gas parcels are all observed to have Doppler-shifted volumes of

$$\begin{eqnarray}&&{V}_{o}=\displaystyle \frac{{V}^{\prime }}{{\rm{\Gamma }}(1+v\mathrm{cos}\theta )}.\end{eqnarray} \tag{ 9 }$$

In this equation, V_o is the observed volume of the gas parcel, ${V}^{\prime }$ is its comoving volume, v is its velocity, and θ is the angle made between the z-axis of the origin and the velocity vector of the moving gas parcel and is, therefore, identical to the normal definition of θ in spherical-polar coordinates.

The perceived optical depth measured to a distance r from the origin is given by

$$\begin{eqnarray}&&\tau ={\int }_{0}^{r}{\rho }_{o}\;\kappa \;d\tilde{r},\end{eqnarray} \tag{ 10 }$$

where $\tilde{r}$ is a dummy variable of integration and ${\rho }_{o}$ is the observed mass density within the fluid; note that this is different from the comoving optical depth, defined by Equation (3), as here we are taking into account light travel time and Lorentz contraction effects (which are encapsulated in the Doppler factor). Since the observed density is related to the observed volume via ${\rho }_{o}\propto 1/{V}_{o}$, Equation (10) becomes

$$\begin{eqnarray}&&\tau ={\int }_{0}^{r}{\rm{\Gamma }}(1+v\mathrm{cos}\theta ){\rho }^{\prime }\;\kappa \;d\tilde{r}.\end{eqnarray} \tag{ 11 }$$

With Equation (2) for the velocity, Equation (11) is

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1}{{\lambda }^{\prime }}{\int }_{0}^{r}\left[\mathrm{cosh}\left(\displaystyle \frac{\mu y}{{\lambda }^{\prime }}\right)+\mathrm{sinh}\left(\displaystyle \frac{\mu y}{{\lambda }^{\prime }}\right)\mathrm{cos}\theta \right]d\tilde{r},\end{eqnarray} \tag{ 12 }$$

where ${\lambda }^{\prime }=1/({\rho }^{\prime }\kappa )$ is the comoving mean free path of the radiation. If we further note that our definition of θ is just that of spherical-polar coordinates, so we can write $y=r\mathrm{sin}\theta \mathrm{sin}\phi$, then this integral can be evaluated and inverted to yield r in terms of τ, μ, and the polar angles. Doing so gives

$$\begin{eqnarray}&&r=\displaystyle \frac{{\lambda }^{\prime }}{\mu \mathrm{sin}\theta \mathrm{sin}\phi }\mathrm{ln}(h),\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}h & = & {(1+\mathrm{cos}\theta )}^{-1}(\mathrm{cos}\theta +\mu \tau \mathrm{sin}\theta \mathrm{sin}\phi \\ & & +\sqrt{1+2\mu \tau \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi +{\mu }^{2}{\tau }^{2}{\mathrm{sin}}^{2}\theta {\mathrm{sin}}^{2}\phi })\end{eqnarray} \tag{ 14 }$$

is a function of μ, τ and the polar angles. Equation (13) gives the perceived distance to neighboring fluid elements in terms of their optical depth τ, which is a two-dimensional surface in θ and ϕ for fixed τ.

We expect that the radiation field at the origin is mainly determined by photons scattered at $\tau \simeq 1$. Photons scattered at $\tau \gtrsim 1$ will, on average, suffer additional scattering before reaching the origin, while relatively few photons will be scattered at $\tau \lesssim 1$.

Figure 1 shows the $\tau =1$ surface in the y−z plane for a number of different μ, which we recall parameterizes the change in the velocity in the y-direction relative to the optical depth, with larger μ corresponding to greater shear (see Equation (2)). The blue circle is the solution for $\mu =0$, the most elongated surface has $\mu =5$, and each curve in between differs from the previous by 1. The arrows in the figure show the direction of the velocity, with the length of the arrow scaling directly with the magnitude of the velocity.

Zoom In Zoom Out Reset image size

This figure demonstrates how the optical depth of the fluid responds to changes in the velocity field: when $\mu =0$, the fluid is stationary everywhere, and correspondingly the $\tau =1$ surface is just a circle where $r={\lambda }^{\prime }$. As the shear starts to increase, the fluid elements that are directly along y at z = 0 are Lorentz contracted because their motion is perpendicular to the line of sight, which increases the density of those fluid elements and brings the $\tau =1$ surface closer to the origin. As we look along y at non-zero locations on the z-axis, the Doppler shift competes with the Lorentz contraction to give a more complicated surface. In particular, for $y\lt 0$ and $z\lt 0$, the fluid is moving away from the origin, and the Doppler shift works with the Lorentz contraction to give an overall smaller fluid volume, increasing the density and moving the $\tau =1$ surface closer to the origin. On the other hand, for $y\gt 0$ and $z\lt 0$, the fluid is moving toward the origin, and the Doppler shift serves to lengthen the fluid element, decreasing the optical depth and extending the $\tau =1$ surface to greater distances from the origin. This behavior is inverted when we consider positions within the fluid characterized by $z\gt 0$.

Using Equations (2) and (13), we find that the fluid four-velocity varies along the $\tau =1$ surface as

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{1}{2}\left(h+\displaystyle \frac{1}{h}\right),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}v=\displaystyle \frac{1}{2}\left(h-\displaystyle \frac{1}{h}\right).\end{eqnarray} \tag{ 16 }$$

where h depends on θ and ϕ via Equation (14). Figure 2 shows the Lorentz factor—Equation (15)—along the $\tau =1$ surface with $\phi =\pi /2$ for the same set of μ used in Figure 1. This figure shows that the Lorentz factor, and correspondingly the velocity, is maximized approximately at the location where the $\tau =1$ surface has receded to the largest distance from the origin, which is consistent with the fact that the Doppler shift is modifying the observed mean free path of the radiation. We will use these expressions in the following section.

Zoom In Zoom Out Reset image size

4.1. Photon Number

Because we are only considering scattering processes, the number of photons must be conserved and, to preserve the similarity of the flow, must also be the same at every point within the fluid. On the $\tau \simeq 1$ surface, we have

$$\begin{eqnarray}&&{n}^{\prime }=\int f\cdot {({k}^{\prime })}^{2}{{dk}}^{\prime }\;d{{\rm{\Omega }}}^{\prime },\end{eqnarray} \tag{ 17 }$$

where f is the photon distribution function (units of number per momentum cubed per length cubed; see, e.g., Mihalas & Mihalas 1984), ${k}^{\prime }$ is the magnitude of the photon three-momentum, and ${{\rm{\Omega }}}^{\prime }$ is the comoving solid angle. A prime on a quantity denotes that it is measured in the comoving frame of the $\tau \simeq 1$ surface. We will adopt the assumption that photons are scattered isotropically in the rest frame of the fluid, and so the distribution function should be independent of angle. It is thus tempting to write $f=f({k}^{\prime });$ however, this assertion is problematic as the radiation field should be heating up if there is a non-zero shear ($\mu \ne 0$), and we therefore expect the photon spectrum to be an increasing function of time. We will therefore adopt the expression

$$\begin{eqnarray}&&f=g(t)q\left(\displaystyle \frac{{k}^{\prime }}{j(t)}\right),\end{eqnarray} \tag{ 18 }$$

where g and j are functions of time and q is an unspecified function. Inserting this expression into Equation (17), we find

$$\begin{eqnarray}&&{n}^{\prime }=4\pi g(t)j{(t)}^{3}\int q(x){x}^{2}{dx},\end{eqnarray} \tag{ 19 }$$

where we have changed variables to $x\equiv {k}^{\prime }/j(t)$. If we now enforce the fact that the photon number density should be time-independent (scattering does not alter the number of photons), then we immediately see that

$$\begin{eqnarray}&&g(t)=j{(t)}^{-3}.\end{eqnarray} \tag{ 20 }$$

The number density of photons at the origin is

$$\begin{eqnarray}&&n=\int {f}_{o}\;{k}^{2}{dk}\;d{\rm{\Omega }},\end{eqnarray} \tag{ 21 }$$

where f_o is the distribution function at the origin. Since the distribution function is a general relativistic invariant (Debbasch & van Leeuwen 2009), we have that the distribution function observed at the origin is the same as that on the scattering surface:

$$\begin{eqnarray}&&{f}_{o}=j{({t}_{r})}^{-3}q\left(\displaystyle \frac{{k}^{\prime }}{j({t}_{r})}\right),\end{eqnarray} \tag{ 22 }$$

where ${t}_{r}=t-r(\theta ,\phi )$ is the retarded time; r is the distance to the $\tau =1$ surface that is a function of angle, specifically given by Equation (13). Inserting this expression into Equation (21) gives

$$\begin{eqnarray}&&n=\int q\left(\displaystyle \frac{{k}^{\prime }}{j({t}_{r})}\right)\displaystyle \frac{{k}^{2}}{j({t}_{r})}d\left(\displaystyle \frac{k}{j({t}_{r})}\right)d{\rm{\Omega }}.\end{eqnarray} \tag{ 23 }$$

Because the photon four-momentum transforms as a vector, we have²

$$\begin{eqnarray}&&{k}^{\prime }={\rm{\Gamma }}k(1+v\mathrm{cos}\theta ),\end{eqnarray} \tag{ 24 }$$

and so Equation (23) becomes

$$\begin{eqnarray}&&n=\displaystyle \frac{{n}^{\prime }}{4\pi }\int {{\mathcal{D}}}^{3}d{\rm{\Omega }},\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{\mathcal{D}}=\displaystyle \frac{1}{{\rm{\Gamma }}(1+v\mathrm{cos}\theta )}\end{eqnarray} \tag{ 26 }$$

is the relativistic Doppler factor—the same one that appears in Equation (9). By using Equations (14)–(16), we can show that the Doppler factor is given by

$$\begin{eqnarray}&&{\mathcal{D}}=\displaystyle \frac{1}{\sqrt{1+2\mu \tau \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi +{\mu }^{2}{\tau }^{2}{\mathrm{sin}}^{2}\theta {\mathrm{sin}}^{2}\phi }},\end{eqnarray} \tag{ 27 }$$

and so we have

$$\begin{eqnarray}n & = & \displaystyle \frac{{n}^{\prime }}{4\pi }\\ & & \times \displaystyle \int \displaystyle \frac{\mathrm{sin}\theta \;d\theta \;d\phi }{{(1+2\mu \tau \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi +{\mu }^{2}{\tau }^{2}{\mathrm{sin}}^{2}\theta {\mathrm{sin}}^{2}\phi )}^{3/2}}.\end{eqnarray} \tag{ 28 }$$

Now, to preserve the self-similarity of the flow, the origin can equally well be considered to be located on the $\tau =1$ surface appropriate to some other region of the fluid, and hence the number of emitted photons must be equal to the number of observed photons. We therefore require that $n={n}^{\prime }$, and the consistency of this approach then demands that

$$\begin{eqnarray}&&\int \displaystyle \frac{\mathrm{sin}\theta \;d\theta \;d\phi }{{(1+2\mu \tau \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi +{\mu }^{2}{\tau }^{2}{\mathrm{sin}}^{2}\theta {\mathrm{sin}}^{2}\phi )}^{3/2}}=4\pi .\end{eqnarray} \tag{ 29 }$$

Remarkably, we can show that the above condition is satisfied for any product $\mu \tau$, and it is therefore an identity. This model, in which the velocity varies according to Equation (2), thus exactly conserves particle number (as a function of time) and yields precisely the same particle number at every comoving point within the flow.

4.2. Photon Energy

The energy density of the radiation scattered by a given location on the $\tau =1$ surface can be written as

$$\begin{eqnarray}&&{e}^{\prime }=\int f\cdot {({k}^{\prime })}^{3}{{dk}}^{\prime }d{{\rm{\Omega }}}^{\prime }=4\pi j(t)\int q(x){x}^{3}{dx},\end{eqnarray} \tag{ 30 }$$

where here we have simply made the same transformation that turned Equation (17) into (19). As was done for the number density of photons, we can also construct the energy density observed at the origin, and thereby scattered by the $\tau =1$ surface, by

$$\begin{eqnarray}&&e=\int f\;{k}^{3}{dk}\;d{\rm{\Omega }}.\end{eqnarray} \tag{ 31 }$$

By now using the same procedure that allowed us to arrive at Equation (21), we can show that this expression becomes

$$\begin{eqnarray}&&e=\displaystyle \frac{{e}^{\prime }}{4\pi }\int \displaystyle \frac{j({t}_{r})}{j(t)}{{\mathcal{D}}}^{4}d{\rm{\Omega }}.\end{eqnarray} \tag{ 32 }$$

Because this is a scattering process that is assumed to be elastic in the rest frame of the scatterer, we require $e={e}^{\prime }$. This requirement then imposes the restriction

$$\begin{eqnarray}&&\int \displaystyle \frac{j({t}_{r})}{j(t)}{{\mathcal{D}}}^{4}d{\rm{\Omega }}=4\pi .\end{eqnarray} \tag{ 33 }$$

Since this equation must be true for all time, we see that j must have the form of an exponential:

$$\begin{eqnarray}&&j=\mathrm{exp}\left(\displaystyle \frac{\nu t}{{\lambda }^{\prime }\tau }\right),\end{eqnarray} \tag{ 34 }$$

where ν is the dimensionless heating rate. We see that, in order for our treatment to be consistent, the heating rate must uniquely yield

$$\begin{eqnarray}&&\int {{\mathcal{D}}}^{4}{e}^{-\nu r/({\lambda }^{\prime }\tau )}d{\rm{\Omega }}=4\pi .\end{eqnarray} \tag{ 35 }$$

Using Equation (27) for ${\mathcal{D}}$ and (13) for r, this relation becomes

$$\begin{eqnarray}&&I(\nu ,\mu \tau )=4\pi ,\end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray}&&I\equiv \int \displaystyle \frac{{h}^{-\tfrac{\nu }{\mu \tau \mathrm{sin}\theta \mathrm{sin}\phi }}\mathrm{sin}\theta \;d\theta \;d\phi }{{(1+2\mu \tau \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi +{\mu }^{2}{\tau }^{2}{\mathrm{sin}}^{2}\theta {\mathrm{sin}}^{2}\phi )}^{2}}\end{eqnarray} \tag{ 37 }$$

and h is given by Equation (14).

From Equation (30), the energy density of the photon field evolves as ${e}^{\prime }\sim \mathrm{exp}(\nu t/({\lambda }^{\prime }\tau ))$, and ν therefore represents the rate at which the shear present in the flow provides energy to the radiation field. This heating rate is not arbitrary, however, as Equation (36) must be an identity for all $\mu \tau$ to preserve the self-similarity of the flow. Ideally we would like to solve Equation (36) for $\nu (\mu ,\tau )$. However, the complexity of the integral I makes this a formidable task, and we must resort to approximate solutions.

5.1. Viscous Limit

When the product $\mu \tau$ is small, we can Taylor-expand Equation (36) in powers of $\mu \tau$ about zero. Doing so to third order and evaluating the integrals, we can show that the heating rate must be

$$\begin{eqnarray}&&\nu =\displaystyle \frac{2}{15}{\mu }^{2}{\tau }^{2}+{ \mathcal O }({\mu }^{4}{\tau }^{4})\end{eqnarray} \tag{ 38 }$$

to satisfy Equation (36) identically (i.e., for all $\mu \tau$). This expression shows that the rate at which the radiation field heats viscously is proportional to the shear squared to lowest order, which is reasonable from a physical standpoint: if there is no shear ($\mu =0$), we do not expect the radiation field to heat up. Likewise, the heating rate should not depend on the sign of μ, and so we would predict that the lowest-order correction to the heating rate is proportional to ${\mu }^{2}$ (and all odd powers of μ should drop out of the expression).

5.2. Streaming Limit

When the change in the Lorentz factor across the photon mean free path is much greater than one, so $\mu \tau \gg 1$, simply Taylor-expanding ν about $\mu \tau =0$ is no longer a valid approach. In this case we must find an alternative method of approximating the integral I.

To this end, note that when $\mu \tau \gg 1$, the denominator in Equation (37) is very large unless $\mu \tau \mathrm{sin}\theta \mathrm{sin}\phi \simeq 1$. Thus, it must be the case that ${\mathcal{D}}$ has a relative maximum near this location, and it is therefore this region of angular space that contributes predominantly to the integral I. By differentiating ${\mathcal{D}}$ with respect to θ and ϕ, we find that the relative maxima occur at the points

$$\begin{eqnarray}&&{\phi }_{m}=\displaystyle \frac{\pi }{2},\displaystyle \frac{3\pi }{2},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&\mathrm{cos}{\theta }_{m}=\pm \sqrt{\displaystyle \frac{1}{2}\left(1+\displaystyle \frac{m}{\sqrt{4+{m}^{2}}}\right)}\simeq \pm \;\left(1-\displaystyle \frac{1}{2{m}^{2}}\right),\end{eqnarray} \tag{ 40 }$$

where for ease of notation we have set $m\equiv \mu \tau ;$ in this expression the negative solution corresponds to $\phi =\pi /2$ and the positive solution to $\phi =3\pi /2$. We can show that the integrand is symmetric about $\phi \to \pi +\phi$ and $\theta \to \pi -\theta$, and for this reason we focus only on the maximum ${\phi }_{m}=3\pi /2$. The last line in Equation (40) follows from the fact that we are interested in the large-m limit of this solution, and this also shows that

$$\begin{eqnarray}&&{\theta }_{m}\simeq \displaystyle \frac{1}{m}+{ \mathcal O }\left(\displaystyle \frac{1}{{m}^{3}}\right).\end{eqnarray} \tag{ 41 }$$

With this value of ${\theta }_{m}$, we can approximate the function ${\mathcal{D}}$ by its second-order Taylor series about the point ${\theta }_{m}$. Calculating the second derivatives of ${\mathcal{D}}$, we can then show that this Taylor series is

$$\begin{eqnarray}&&{\mathcal{D}}\simeq m-\displaystyle \frac{{m}^{5}}{2}{\left(\theta -\displaystyle \frac{1}{m}\right)}^{2}-\displaystyle \frac{m}{2}{\left(\phi -\displaystyle \frac{3\pi }{2}\right)}^{2}.\end{eqnarray} \tag{ 42 }$$

This expression shows that the angular width in the θ-direction subtended by ${\mathcal{D}}$ about its maximum value is

$$\begin{eqnarray}&&{\rm{\Delta }}{\theta }_{m}\simeq \displaystyle \frac{1}{{m}^{2}},\end{eqnarray} \tag{ 43 }$$

while ${\rm{\Delta }}{\phi }_{m}$—the angular width in the ϕ-direction—is of order unity. From this expression, it then follows that the angular integral of ${\mathcal{D}}$ is approximately

$$\begin{eqnarray}&&\int {\mathcal{D}}d\theta d\phi \simeq m{\rm{\Delta }}{\theta }_{m}{\rm{\Delta }}{\phi }_{m}\simeq \displaystyle \frac{1}{m}.\end{eqnarray} \tag{ 44 }$$

We would like to follow a similar procedure for approximating the integral in Equation (36), but the integrand ${\mathcal{I}}\equiv {{\mathcal{D}}}^{4}\mathrm{exp}[-\nu r/(\lambda \tau )]\mathrm{sin}\theta$ is clearly much more complicated than just ${\mathcal{D}}$. The maximum value of ${\mathcal{I}}$ will therefore not coincide exactly with the point ${\theta }_{m}=1/m$. However, it is possible to show that the function

$$\begin{eqnarray}&&f\equiv {{\mathcal{D}}}^{4}\mathrm{sin}\theta \end{eqnarray} \tag{ 45 }$$

has its maximum value at the point

$$\begin{eqnarray}&&{\theta }_{f}\simeq 1/m+{ \mathcal O }(1/{m}^{3}),\end{eqnarray} \tag{ 46 }$$

while we find numerically that the function

$$\begin{eqnarray}&&g\equiv {e}^{-\tfrac{\nu r}{\lambda \tau }}\end{eqnarray} \tag{ 47 }$$

has a relative minimum near the point $\theta \simeq 1/m$. It is also straightforward to show that the extrema of both of these functions occur at $\phi =\pi /2,3\pi /2$.

Since ${\mathcal{I}}$ is the product of the functions f and g, the former possessing a relative maximum at the point ${\theta }_{m}\simeq 1/m$, the latter a relative minimum at that point, it is not immediately obvious that the relative extremum of ${\mathcal{I}}$ (which does occur at ${\theta }_{m}\simeq 1/m$) will be a maximum. However, we find that the absolute value of the second derivatives of f at ${\theta }_{m}$ are much larger than the second derivatives of g. Therefore, the second derivatives of ${\mathcal{I}}$, which are just the sum of the second derivatives of f and g, are dominated by the function f, meaning that the extremum of I is indeed a maximum. We therefore have

$$\begin{eqnarray}&&{\mathcal{I}}\simeq {{\mathcal{I}}}_{m}+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}{{\mathcal{I}}}_{m}}{\partial {\theta }^{2}}{\left(\theta -\displaystyle \frac{1}{m}\right)}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}{{\mathcal{I}}}_{m}}{\partial {\phi }^{2}}{\left(\phi -\displaystyle \frac{3\pi }{2}\right)}^{2},\end{eqnarray} \tag{ 48 }$$

where subscripts m denote that we are evaluating the function at the maximum. We can show that the derivatives, to lowest order in $1/m$, are

$$\begin{eqnarray}&&{{\mathcal{I}}}_{m}\simeq {m}^{3-\nu },\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{{\mathcal{I}}}_{m}}{\partial {\theta }^{2}}\simeq -{m}^{7-\nu },\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{{\mathcal{I}}}_{m}}{\partial {\phi }^{2}}\simeq -{m}^{3-\nu },\end{eqnarray} \tag{ 51 }$$

which demonstrates, as was true for ${\mathcal{D}}$ itself, that the width of the maximum of ${\mathcal{I}}$ is ${\rm{\Delta }}{\theta }_{m}\simeq 1/{m}^{2}$ and ${\rm{\Delta }}{\phi }_{m}\simeq 1$. It then follows that

$$\begin{eqnarray}&&I\simeq \int {\mathcal{I}}d\theta d\phi \simeq {{\mathcal{I}}}_{m}{\rm{\Delta }}{\theta }_{m}{\rm{\Delta }}{\phi }_{m}\simeq {m}^{1-\nu }.\end{eqnarray} \tag{ 52 }$$

From Equation (36), the integral I must be equal to $4\pi$, so $\nu (m)$ must satisfy

$$\begin{eqnarray}&&{m}^{1-\nu }\simeq 4\pi .\end{eqnarray} \tag{ 53 }$$

This shows that, in the large-m limit, ν is given by

$$\begin{eqnarray}&&\nu =1-\displaystyle \frac{C}{\mathrm{ln}(\mu \tau )},\end{eqnarray} \tag{ 54 }$$

where C is a numerical constant. Since we only took the leading-order (in $1/m$) expressions for ${{\mathcal{I}}}_{m}$ and its derivatives, the precise value of C cannot be directly computed here. This approach does show, however, that the heating rate must be equal to exactly one in the limit that $\mu \tau \to \infty$, but the convergence is slow ($\propto 1/\mathrm{ln}(\mu \tau )$).

5.3. Interpolated Solution

In the above two subsections we found the following asymptotic limits for the heating rate:

$$\begin{eqnarray}\nu =\left\{\begin{array}{ll}\tfrac{2}{15}{\mu }^{2}{\tau }^{2} & \mathrm{for}\ \mu \tau \ll 1,\\ 1-\frac{C}{\mathrm{ln}(\mu \tau )} & \mathrm{for}\ \mu \tau \gg 1.\end{array}\right.\end{eqnarray} \tag{ 55 }$$

To determine the μ dependence of ν in between these limits, we will write

$$\begin{eqnarray}&&\nu =1-\displaystyle \frac{1}{1+\tfrac{1}{2C}\mathrm{ln}\left(\tfrac{P(\mu \tau )}{Q(\mu \tau )}\right)}.\end{eqnarray} \tag{ 56 }$$

where P and Q are polynomials in the quantity $\mu \tau$. We can then determine these polynomials by requiring that Equation (56) reduce correctly to the viscous limit to a predetermined order. For example, if we only want to match the lowest-order viscous approximation to the heating rate, so $\nu =2{\mu }^{2}{\tau }^{2}/15$, then we can show that

$$\begin{eqnarray}&&P(\mu \tau )=1+\displaystyle \frac{4C}{15}{\mu }^{2}{\tau }^{2},\end{eqnarray} \tag{ 57 }$$

and

$$\begin{eqnarray}&&Q(\mu \tau )=1.\end{eqnarray} \tag{ 58 }$$

The expression for the heating rate that matches both the $\mu \tau \gg 1$ and $\mu \tau \ll 1$ limits is then

$$\begin{eqnarray}&&\nu =1-\displaystyle \frac{1}{1+\tfrac{1}{2C}\mathrm{ln}\left(1+\tfrac{4C}{15}{\mu }^{2}{\tau }^{2}\right)}.\end{eqnarray} \tag{ 59 }$$

To determine C we have adopted a brute-force method of numerically integrating the left-hand side of Equation (36) with (56) for ν for a number of different C. We find that the value of C that solves Equation (36) when $\mu \tau \gg 1$ is

$$\begin{eqnarray}&&C\simeq 0.812.\end{eqnarray} \tag{ 60 }$$

Figure 3 shows the integral in Equation (36), I, normalized by $4\pi$ when ν is given by Equation (59) and C = 0.812. It is apparent that this expression for the heating rate almost exactly satisfies the integral constraint (36), with the maximum deviation from unity being 0.98 at $\mu \tau \simeq 10$.

Zoom In Zoom Out Reset image size

If we want to match higher-order viscous corrections, then it is apparent that P and Q will be of the general form

$$\begin{eqnarray}&&P(\mu \tau )=\displaystyle \sum _{n=0}^{j}{p}_{n}{(\mu \tau )}^{n},\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&Q(\mu \tau )=\displaystyle \sum _{n=0}^{{\ell }}{q}_{n}{(\mu \tau )}^{n}.\end{eqnarray} \tag{ 62 }$$

Because the heating rate should only depend on even powers of $\mu \tau$, it follows that the odd coefficients in these expansions are zero. Likewise, since the asymptotic limit should be $\sim 1-C/\mathrm{ln}(\mu \tau )$, we find $j={\ell }+2$. Thus the next order approximation to the heating rate is

$$\begin{eqnarray}&&\nu =1-\displaystyle \frac{1}{1+\tfrac{1}{2C}\mathrm{ln}\left(\tfrac{1+{p}_{2}{\mu }^{2}{\tau }^{2}+{p}_{4}{\mu }^{4}{\tau }^{4}}{1+{q}_{2}{\mu }^{2}{\tau }^{2}}\right)},\end{eqnarray} \tag{ 63 }$$

where the coefficients can be determined by equating the Taylor series of this function to the viscous approximation of the heating rate. Notice that, because we have three unknowns here, we must expand the viscous limit to sixth order. Similarly, the next highest order will have five unknowns, meaning that we need to expand the viscous limit to tenth order, etc.

The values of the coefficients p_n and q_n will depend on where we truncate the viscous approximation to the heating rate, and thus our higher-order interpolated solutions for ν will differ from those at lower orders. However, we find that the overall solution appears to converge very rapidly to a unique solution, as depicted by Figure 4. This figure shows the lowest-order interpolated solution (i.e., Equation (59); blue curve) and the solution accurate to sixth order (yellow curve), both with C = 0.812. We see by eye that there are only very small differences between these heating rates. Therefore, we can, to a very high degree of accuracy, neglect the higher-order viscous corrections and take Equation (59) to be the heating rate that correctly reproduces both the viscous and streaming limits of radiation propagation and interpolates well between these extremes. Using the fact that $C\simeq 0.812$, the heating rate that preserves the self-similarity of the flow is

$$\begin{eqnarray}&&\nu =1-\displaystyle \frac{1}{1+0.616\mathrm{ln}(1+0.217{\mu }^{2}{\tau }^{2})}.\end{eqnarray} \tag{ 64 }$$

Zoom In Zoom Out Reset image size

6.1. Comoving Frame

The energy–momentum tensor of the radiation field is (Mihalas & Mihalas 1984)

$$\begin{eqnarray}&&{R}^{\mu \nu }=\int f\;{k}^{\mu }{k}^{\nu }\displaystyle \frac{{d}^{3}k}{{k}^{0}}.\end{eqnarray} \tag{ 65 }$$

We can evaluate this tensor at the origin by recalling that the distribution function is given by Equation (22), using the fact that ${k}^{x}=k\mathrm{cos}\theta$ and ${k}^{y}=k\mathrm{sin}\theta \mathrm{sin}\phi$, and using Equation (24) to write k in terms of ${k}^{\prime }$ and θ. Doing so then gives

$$\begin{eqnarray}&&{R}^{\bar{\mu }\bar{\nu }}=\displaystyle \frac{{e}^{\prime }}{4\pi }\int {S}^{\bar{\mu }\bar{\nu }}{{\mathcal{D}}}^{4}{h}^{-\tfrac{\nu }{m\mathrm{sin}\theta \mathrm{sin}\phi }}d{\rm{\Omega }},\end{eqnarray} \tag{ 66 }$$

where

$$\begin{eqnarray}{S}^{\bar{\mu }\bar{\nu }}\;= & \left(\begin{array}{cccc}1 & \mathrm{sin}\theta \mathrm{cos}\phi & \mathrm{sin}\theta \mathrm{sin}\phi & \mathrm{cos}\theta \\ \mathrm{sin}\theta \mathrm{cos}\phi & {\mathrm{sin}}^{2}\theta {\mathrm{cos}}^{2}\phi & {\mathrm{sin}}^{2}\theta \mathrm{sin}\phi \mathrm{cos}\phi & \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{cos}\phi \\ \mathrm{sin}\theta \mathrm{sin}\phi & {\mathrm{sin}}^{2}\theta \mathrm{sin}\phi \mathrm{cos}\phi & {\mathrm{sin}}^{2}\theta {\mathrm{sin}}^{2}\phi & \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi \\ \mathrm{cos}\theta & \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{cos}\phi & \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{sin}\phi & {\mathrm{cos}}^{2}\theta \end{array}\right).\end{eqnarray} \tag{ 67 }$$

We have placed bars on these tensors because they are evaluated in the comoving frame of the gas parcel, and Equation (66), therefore, represents the comoving radiation energy–momentum tensor.

The heating rate ν is constructed such that ${R}^{\bar{0}\bar{0}}={e}^{\prime }$, and the symmetry of the integrand means that the energy fluxes ${R}^{\bar{0}\bar{z}}$, ${R}^{\bar{0}\bar{y}}$, and ${R}^{\bar{0}\bar{x}}$, and the stresses ${R}^{\bar{y}\bar{x}}$ and ${R}^{\bar{z}\bar{x}};$ are zero. It is also straightforward to show that the x-component of the pressure is equal to ${R}^{\bar{x}\bar{x}}={R}^{\bar{0}\bar{0}}-{R}^{\bar{z}\bar{z}}-{R}^{\bar{y}\bar{y}}$. The only non-zero and non-trivial components of the stress tensor are therefore the pressures, ${R}^{\bar{z}\bar{z}}$ and ${R}^{\bar{y}\bar{y}}$, and the shear stress ${R}^{\bar{y}\bar{z}}$.

Figure 5 shows these three quantities—the absolute value of the shear stress (solid, blue curve), the z-momentum flux (the z-component of the pressure; dashed, red curve) and the y-momentum flux (the y-component of the pressure; dotted–dashed, purple curve) all normalized by the energy density ${e}^{\prime }$—as functions of $\mu \tau$. When the shear across the photon mean free path is small, so $\mu \ll 1$, we expect the radiation field to reduce to the optically thick limit. We see from Figure 5 that this consistency check is met: the shear stress approaches zero, and the momentum fluxes are ${R}^{\bar{z}\bar{z}}={R}^{\bar{y}\bar{y}}={e}^{\prime }/3$—the value associated with an isotropic photon gas. As μ increases, we see that the shear stress initially rises in a linear fashion, reaches a peak value of ${R}^{\bar{y}\bar{z}}\simeq 0.25$ at $\mu \tau \simeq 1.8$, and decays asymptotically as ${R}^{\bar{y}\bar{z}}\propto 1/(\mu \tau );$ the y-component of the pressure falls off rather steeply with the shear, being well approximated by ${R}^{\bar{y}\bar{y}}\simeq 1/{(\mu \tau )}^{2}$ for $\mu \gtrsim 1;$ and the z-component of the pressure, while it rises steeply initially, eventually levels off to ${R}^{\bar{z}\bar{z}}={e}^{\prime }$ for $\mu \tau \gg 1$.

Zoom In Zoom Out Reset image size

The solutions for these quantities were obtained by numerically integrating Equation (66) (with Equation (64) for the heating rate), and therefore have a complicated dependence on the quantity $\mu \tau$. However, we find that each can be fit well by fairly simple functions. In particular, we find that

$$\begin{eqnarray}&&{R}^{\bar{y}\bar{z}}=-\displaystyle \frac{4}{15}\tfrac{\mu \tau \;{e}^{\prime }}{1+\tfrac{3}{10}{\mu }^{2}{\tau }^{2}},\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&{R}^{\bar{z}\bar{z}}=\left(1-\displaystyle \frac{2}{3}\tfrac{\tfrac{10}{9}\mathrm{tanh}\left(\tfrac{9}{10}\mu \tau \right)}{\mu \tau }\right){e}^{\prime },\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&{R}^{\bar{y}\bar{y}}=\displaystyle \frac{1}{3}\tfrac{{e}^{\prime }}{1+\tfrac{1}{3}{\mu }^{2}{\tau }^{2}}.\end{eqnarray} \tag{ 70 }$$

All of the properties of the radiation field that we have investigated thus far have been functions of the combination $\mu \tau$, thereby rendering the exact value of τ relatively unimportant. Equation (68), however, gives us one means to determine its value: in the viscous limit, this equation gives

$$\begin{eqnarray}&&{R}^{\bar{y}\bar{z}}=-\displaystyle \frac{4}{15}\mu \tau +{ \mathcal O }({\mu }^{3}{\tau }^{3}),\end{eqnarray} \tag{ 71 }$$

which is, self-consistently, the relation we determine by Taylor-expanding Equation (66) to first order in $\mu \tau$. On the other hand, recent investigations of the equations of radiation hydrodynamics in the viscous limit have shown that (Coughlin & Begelman 2014b; see also Blandford et al. 1985)

$$\begin{eqnarray}&&{R}^{\bar{y}\bar{z}}=-\displaystyle \frac{8}{27}\mu .\end{eqnarray} \tag{ 72 }$$

Comparing this expression with Equation (71), we see that the value of τ that ensures that our self-similar approach is consistent with the equations of radiation hydrodynamics in the viscous limit is

$$\begin{eqnarray}&&{\tau }_{}=\displaystyle \frac{10}{9}.\end{eqnarray} \tag{ 73 }$$

However, it should be noted that Coughlin & Begelman (2014b) used a dipole scattering kernel to treat the interactions between the photons and the electrons in the gas. If one uses a spherical kernel, which is more appropriate for our analysis here (as we assumed that the distribution function was isotropic in the rest frame of the $\tau \simeq 1$ surface; see Section 4), then one finds (Weinberg 1971)

$$\begin{eqnarray}&&{R}^{\bar{y}\bar{z}}=-\displaystyle \frac{4}{15}\mu .\end{eqnarray} \tag{ 74 }$$

If we adopt this expression for the viscous limit of the radiation energy–momentum tensor, then we see that the appropriate optical depth is

$$\begin{eqnarray}&&\tau =1.\end{eqnarray} \tag{ 75 }$$

6.2. Covariant Formulation

Equation (66) gives the radiation energy–momentum tensor in the comoving frame of the fluid. However, the equations of radiation hydrodynamics, which govern the interaction between the scatterers and the radiation field, are given by

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\alpha }({T}^{\alpha \beta }+{R}^{\alpha \beta })=0,\end{eqnarray} \tag{ 76 }$$

where ${{\rm{\nabla }}}_{\alpha }$ is the covariant derivative and ${T}^{\alpha \beta }$ is the energy–momentum tensor of the scatterers. It is therefore necessary to differentiate the energy–momentum tensor of the radiation field, which is tantamount to evaluating this tensor at different locations within the fluid, and hence just knowing its form at the origin is insufficient. However, from the self-similar nature of the shear flow, we know that at any comoving point the radiation field must have the form given by (66). Since the comoving frame can be obtained by making a local Lorentz transformation, we therefore have

$$\begin{eqnarray}&&{R}^{\alpha \beta }={{\rm{\Lambda }}}_{\;\bar{\alpha }}^{\alpha }{{\rm{\Lambda }}}_{\;\bar{\beta }}^{\beta }{R}^{\bar{\alpha }\bar{\beta }}.\end{eqnarray} \tag{ 77 }$$

where

$$\begin{eqnarray}{{\rm{\Lambda }}}_{\;\mu }^{\bar{\mu }}=\left(\begin{array}{cccc}{\rm{\Gamma }} & 0 & 0 & {\rm{\Gamma }}v\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ {\rm{\Gamma }}v & 0 & 0 & {\rm{\Gamma }}\end{array}\right).\end{eqnarray} \tag{ 78 }$$

Carrying out the multiplication, this gives

$$\begin{eqnarray}{R}^{\alpha \beta }=\left(\begin{array}{cccc}{{\rm{\Gamma }}}^{2}({e}^{\prime }+{v}^{2}{R}^{\bar{z}\bar{z}}) & 0 & {\rm{\Gamma }}{{vR}}^{\bar{y}\bar{z}} & {{\rm{\Gamma }}}^{2}v({e}^{\prime }+{R}^{\bar{z}\bar{z}})\\ 0 & {R}^{\bar{x}\bar{x}} & 0 & 0\\ {\rm{\Gamma }}{{vR}}^{\bar{y}\bar{z}} & 0 & {R}^{\bar{y}\bar{y}} & {\rm{\Gamma }}{R}^{\bar{y}\bar{z}}\\ {{\rm{\Gamma }}}^{2}v({e}^{\prime }+{R}^{\bar{z}\bar{z}}) & 0 & {\rm{\Gamma }}{R}^{\bar{y}\bar{z}} & {{\rm{\Gamma }}}^{2}({v}^{2}{e}^{\prime }+{R}^{\bar{z}\bar{z}})\end{array}\right).\end{eqnarray} \tag{ 79 }$$

We can show that this matrix can be written

$$\begin{eqnarray}{R}^{\alpha \beta } & = & e^{\prime} {U}^{\alpha }{U}^{\beta }+{R}^{\bar{z}\bar{z}}{{\rm{\Pi }}}^{\alpha \beta }\\ & & +\displaystyle \frac{\lambda ^{\prime} }{\mu }{R}^{\bar{y}\bar{z}}{{\rm{\Pi }}}^{\alpha \sigma }{{\rm{\Pi }}}^{\beta \rho }({{\rm{\nabla }}}_{\sigma }{U}_{\rho }+{{\rm{\nabla }}}_{\rho }{U}_{\sigma })+{\rm{\Delta }}{P}^{\alpha \beta },\end{eqnarray} \tag{ 80 }$$

where

$$\begin{eqnarray}&&{{\rm{\Pi }}}^{\alpha \beta }={U}^{\alpha }{U}^{\beta }+{g}^{\alpha \beta }\end{eqnarray} \tag{ 81 }$$

is the projection tensor and we defined

$$\begin{eqnarray}{\rm{\Delta }}{P}^{\alpha \beta }\equiv \left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & {R}^{\bar{x}\bar{x}}-{R}^{\bar{z}\bar{z}} & 0 & 0\\ 0 & 0 & {R}^{\bar{y}\bar{y}}-{R}^{\bar{z}\bar{z}} & 0\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ 82 }$$

which measures the degree of pressure anisotropy exhibited by the radiation field.

When the amount of shear present in the flow is vanishingly small, the components of the comoving radiation energy–momentum tensor reduce to ${R}^{\bar{y}\bar{z}}\simeq 0$, ${R}^{\bar{z}\bar{z}}\simeq {R}^{\bar{y}\bar{y}}\simeq {R}^{\bar{x}\bar{x}}\simeq {e}^{\prime }/3$. Using these expressions in Equation (80), we see that we self-consistently recover the result for the energy–momentum tensor of an isotropic, relativistic fluid:

$$\begin{eqnarray}&&{R}^{\alpha \beta }={e}^{\prime }{U}^{\alpha }{U}^{\beta }+\displaystyle \frac{1}{3}{e}^{\prime }{{\rm{\Pi }}}^{\alpha \beta }.\end{eqnarray} \tag{ 83 }$$

If we keep first-order shear corrections to the comoving energy–momentum tensor, so we maintain pressure anisotropy but have ${R}^{\bar{y}\bar{z}}=-4\mu \tau {e}^{\prime }/15$, then we find

$$\begin{eqnarray}&&{R}^{\alpha \beta }={e}^{\prime }{U}^{\alpha }{U}^{\beta }+\displaystyle \frac{1}{3}{e}^{\prime }{{\rm{\Pi }}}^{\alpha \beta }-\displaystyle \frac{4{e}^{\prime }\tau }{15{\rho }^{\prime }\kappa }{{\rm{\Pi }}}^{\alpha \sigma }{{\rm{\Pi }}}^{\beta \rho }({{\rm{\nabla }}}_{\sigma }{U}_{\rho }+{{\rm{\nabla }}}_{\rho }{U}_{\sigma }).\end{eqnarray} \tag{ 84 }$$

Comparing this expression to Equation (37) of Coughlin & Begelman (2014b), we see that this approach to analyzing relativistic shear is in agreement with the viscous equations of radiation hydrodynamics for divergenceless flow (${{\rm{\nabla }}}_{\mu }{U}^{\mu }=0$) if we adopt $\tau =10/9$.

On the other hand, when the amount of shear becomes very large across the mean free path of the photon ($\mu \gg 1$), we have ${R}^{\bar{x}\bar{x}}\simeq {R}^{\bar{y}\bar{y}}\simeq {R}^{\bar{y}\bar{z}}\simeq 0$ and ${R}^{\bar{z}\bar{z}}\simeq {e}^{\prime }$. Furthermore, since the flow is highly relativistic in this case even in the immediate vicinity of any comoving gas parcel, we have ${{\rm{\Gamma }}}^{2}\simeq {{\rm{\Gamma }}}^{2}{v}^{2}$ and it therefore follows from Equation (80) that

$$\begin{eqnarray}&&{R}^{\alpha \beta }=2{e}^{\prime }{U}^{\alpha }{U}^{\beta }.\end{eqnarray} \tag{ 85 }$$

When the shear is between these two limits, relatively simple covariant expressions for the energy–momentum tensor of the radiation field are not able to be obtained. However, we note that Equation (80) is valid for arbitrary μ, but the full μ-dependence of the comoving stress tensor must be incorporated. Furthermore, by using Equation (2), we can show that the shear parameter μ can be expressed as a covariant scalar via:

$$\begin{eqnarray}&&{\mu }^{2}={({\lambda }^{\prime })}^{2}{{\rm{\Pi }}}^{\mu \sigma }({{\rm{\nabla }}}_{\mu }{U}_{\nu })({{\rm{\nabla }}}_{\sigma }{U}^{\nu }).\end{eqnarray} \tag{ 86 }$$

The right-hand side of this equation can be interpreted as the relativistic "square" of the shear over the mean free path of the photon.

In this paper we analyzed how a radiation field responds to regions of intense, relativistic shear, which likely arise in extreme astrophysical environments such as collapsars (MacFadyen & Woosley 1999) and tidal disruption events (Coughlin & Begelman 2014a). We considered a two-dimensional, planar shear flow, with the motion along the z-direction and the variation in that motion along the y-direction, in which the fluid appears identical at every comoving point. We demonstrated that this self-similar assumption requires the velocity profile of the fluid to have the form ${\rm{\Gamma }}=\mathrm{cosh}(\mu {\tau }_{y}^{\prime })$ with μ a constant and ${\tau }_{y}^{\prime }={\rho }^{\prime }\kappa \;y$ (and the comoving density ${\rho }^{\prime }$ is a constant to preserve the self-similarity of the flow). Using this velocity field, we determined the $\tau \simeq 1$ surface—the location within the fluid where the integrated optical depth along the line of sight equals roughly one—which is a complicated function of viewing angle and shear owing to relativistic Doppler beaming (see Figure 1 and Equation (13)).

Using the structure of the $\tau \simeq 1$ surface and the assumption that photons are scattered isotropically in the rest frame of the scatterer, we showed that this type of shear flow exactly conserves photon number if the distribution function of the radiation field is given by Equation (18). We also demonstrated that, if the energy density is to be uniform throughout the shear layer, which must be true if the scattering is elastic and the self-similarity of the fluid is upheld, then the radiation field must be heating up exponentially with a heating rate that is given by the solution to an integro-algebraic equation (see Equations (36) and (37)). The solutions to this equation were determined in the limits of small and large shear, and an approximate heating rate that interpolates between these limits and almost exactly satisfies the integro-algebraic equation (see Figure 3) was found (Equation (64)).

Finally, we constructed the comoving energy–momentum tensor of the radiation field, the only non-zero and non-trivial components of which were the shear, ${R}^{\bar{y}\bar{z}}$, and the z- and y-components of the pressure, ${R}^{\bar{z}\bar{z}}$ and ${R}^{\bar{y}\bar{y}}$. These quantities were shown to smoothly transition from their viscous (small shear) to their streaming (large shear) limits, the former characterized by a shear stress that varies linearly with the shear (i.e., Newtonian in nature) and isotropic pressure, the latter portraying vanishing shear stress and highly anisotropic pressure (${R}^{\bar{y}\bar{y}}\simeq {R}^{\bar{x}\bar{x}}\simeq 0$, ${R}^{\bar{z}\bar{z}}\simeq {e}^{\prime };$ see Figure 5). We showed that ${R}^{\bar{y}\bar{z}}$, ${R}^{\bar{y}\bar{y}}$, and ${R}^{\bar{z}\bar{z}}$ were very well-fit by approximate, analytic functions (see Equations (68)–(70)), and found that the value of the optical depth must be equal to $\tau =10/9$ if the viscous limits of our equations match those pursued by other authors (though a value of $\tau =1$ matches the results if an isotropic scattering kernel is used in the Boltzmann equation). By using the self-similarity of the fluid, we were able to construct the form of the energy–momentum tensor at any point within the flow, showing that the result agreed with the isotropic and viscous limits.

Figure 5, together with Equations (68)–(70), is perhaps the most important result of this investigation: this figure shows how the radiation field in a relativistically moving plasma transitions from the viscous to the streaming limit. Interestingly, the shear stress, ${R}^{\bar{y}\bar{z}}$, does not grow to arbitrarily large values as the shear increases (as one might predict from its linear growth at small μ), but reaches a maximum when $\mu \simeq 1.8$ and thereafter decays approximately as $1/\mu$. The origin of this behavior can be understood as follows: the shear stress ${R}^{\bar{y}\bar{z}}$ gives the amount of y-momentum transferred in the z-direction. As the shear starts to increase, the structure of the $\tau =1$ surface deviates from a sphere at $r={\lambda }^{\prime }$, allowing photons with negative y-momentum to be transferred in the positive z-direction (and vice versa; see Figure 1) and generating the stress. When the shear becomes very large, however, the surface becomes increasingly aligned with the z-axis as a consequence of relativistic Doppler beaming, meaning that only photons possessing a small amount of transverse momentum are perceived in the comoving frame. The radiation field thus becomes highly beamed in the direction of motion of the fluid, thus inhibiting the transfer of transverse momentum to the scatterers. This behavior also shows that the effects of radiation drag can actually be quenched in an optically thin, relativistic shear layer, ultimately due to the fact that the perceived radiation field adapts to the presence of the shear itself.

Figure 5 also demonstrates that the z-component of the pressure approaches the energy density as the shear becomes very large, as one would suspect for a highly beamed source. In fact, the red dashed curve in that figure represents the Eddington factor $f(\mu \tau )$, which relates the z-component of the pressure to the energy density via ${R}^{\bar{z}\bar{z}}=f(\mu \tau ){e}^{\prime }$. Investigating Equation (69), we thus see that shear in a scattering medium generates an effective Eddington factor that can be well approximated by

$$\begin{eqnarray}&&f(\mu )=1-\displaystyle \frac{2}{3}\displaystyle \frac{\mathrm{tanh}(\mu )}{\mu },\end{eqnarray} \tag{ 87 }$$

where we have set $\tau =10/9$ (which, as we demonstrated above, is the value we expect if this method is to reduce to the viscous limit when the scattering is modeled by a dipole kernel).

Our approach adopted a very specific form for the velocity of the fluid within the shear layer. As we argued above, this velocity has the property that the fluid looks the same at every comoving point, meaning that such a shear layer may develop naturally in regions where the fluid has a "lost memory" of the boundary conditions or initial conditions. In further support of this notion, we can show that this exact velocity profile develops in the treatment of radiation–viscous boundary layers when the flow becomes ultrarelativistic, and we refer the reader to the Appendix for a demonstration of this fact.

The investigation we have undertaken here is related to the radiation drag limit of the interaction between photons and matter, where the radiation field is considered a constant background that is unaltered by the scattering processes that take place. On the other hand, here we have accounted for the evolution of the radiation field—it is heated, exponentially so, by the shear present in the flow. However, we have ignored the back reaction that this heating (and the viscous stress) must have on the flow itself; similar to the fictitious entity that maintains the isotropy of the radiation field in the radiation drag limit, there must be some external force in our model that maintains the shear profile of the self-similar flow. It is the work done by this force that ultimately heats the radiation field.

We also reiterate the point made in the Introduction and in Section 2: owing to its self-similar nature, certain regions of astrophysical plasmas may naturally conform to the shear profile given by Equation (2). In such cases, the amount of shear present in the flow is maintained by the non-self-similar aspects of the problem, being, e.g., the boundary conditions or gradients along the direction of motion of the fluid. We verified this notion by showing, in the Appendix, that the velocity profile of the two-stream, radiation–viscous boundary layer manifestly yields the self-similar form given by Equation (2). In this case it is the boundary conditions—that the fluid match the speed of the "jet" in one limit and approach the static envelope in another limit—that provide and maintain the shear of the flow. Nevertheless, far away from these boundaries the specifics of those constraints are lost and the velocity transitions to the self-similar flow field of Equation (2).

In our treatment, the comoving density, ${\rho }^{\prime }$, was considered independent of position and time, which is fundamental to the assumption of self-similarity within the boundary layer. However, it is possible to generalize this aspect of the problem by simply keeping the integral expression for the variable ${\tau }_{y}^{\prime }$, which appears in Equation (3). Doing so, we can construct the $\tau \simeq 1$ surface in an identical manner to what was done in Section 3, integrate the expression exactly, and solve for the radial position of the surface as a function of angle, which gives

$$\begin{eqnarray}&&{\int }_{0}^{r\mathrm{sin}\theta \mathrm{sin}\phi }{\rho }^{\prime }\kappa \;d\tilde{y}={\tau }_{y}^{\prime }=\displaystyle \frac{1}{\mu \mathrm{sin}\theta \mathrm{sin}\phi }\mathrm{ln}(h),\end{eqnarray} \tag{ 88 }$$

where h is still given by Equation (14). What this finding demonstrates is that, since the velocity is only a function of $\mu {\tau }_{y}^{\prime }$ (Equation (2)), the same vector transformations to obtain the comoving properties of the radiation on the $\tau =1$ surface hold, and thus the particle number is still exactly conserved if we make the same ansatz for the distribution function, Equation (18). However, the heating rate will not, in general, be the same, as this property of the radiation field depends on the retarded time between the comoving frame and the $\tau =1$ surface. We therefore must assume a specific form for ${\rho }^{\prime }$ to calculate ν.

Fukue (2008) followed a similar procedure to what we outlined here for calculating the properties of the radiation field: he constructed a surface in the comoving frame of the fluid that satisfied $\tau =1$ (his "one-tau photo-oval"), and he calculated the properties of the radiation field in the comoving frame based on the appearance of the field on the $\tau =1$ surface. However, his treatment assumed that the fluid satisfied ${\boldsymbol{v}}=v(z)\hat{z}$, i.e., a one-dimensional flow in which there is no transverse shear. He also let the comoving velocity field be $v(z)={v}_{0}+({dv}/{dz})(z-{z}_{0})$ and similarly for other fluid quantities, meaning that his results are only valid when the shear over the mean free path is small. The heating of the radiation field due to relativistic time delays was also ignored in his model, which was an essential aspect of our formulation.

Most radiation hydrodynamics codes employ a closure scheme—either some variant of M1 or FLD—to calculate some of the moments (namely the shear stresses and the pressures) of the radiation energy–momentum tensor. However, there are a few authors who have opted to directly solve the radiative transfer equation (Jiang et al. 2014; Ryan et al. 2015) alongside the equations of radiation hydrodynamics (or MHD), thereby directly computing the moments of the radiation field and coupling them to the equations of motion (and vice versa). While the methods we have outlined here were not directly based on the relativistic transfer equation, many of the properties of the fluid and the radiation field—the self-similar appearance of the flow, the consistent transition between the optically thick and thin limits, the conservation of energy and particle number—should also result from an analysis of the Boltzmann equation. We therefore feel that Equations (68)–(70) (or something close to them) should arise from an investigation of the relativistic transfer equation, and hence the comoving components of the radiation stress tensor obtained in this paper can be considered as a test for relativistic radiation–MHD solvers.

Our results concerning the properties of the radiation field were rigorously obtained only for flow in which the velocity varies as ${\rm{\Gamma }}=\mathrm{cosh}(\mu {\tau }_{y}^{\prime })$. However, investigating Equation (80), we see that the form for the stress tensor that we derived only depends on this assumption through the inclusion of the parameter μ (and its assumed-constant nature). Furthermore, if we recall that μ can be characterized as a covariant scalar via Equation (86) that also does not depend on the explicit form for the velocity profile, then we can plausibly interpret Equation (80) as the radiation stress tensor that is valid for more general, but divergenceless, flows when the shear becomes large. If we additionally want to include flows that contain non-zero divergence, then comparison between our Equations (80) and (37) of Coughlin & Begelman (2014b) suggests that the generalization of the viscous stress tensor to arbitrary shear is

$$\begin{eqnarray}&&{R}^{\alpha \beta }=-\displaystyle \frac{8}{27}\displaystyle \frac{{e}^{\prime }}{{\rho }^{\prime }\kappa }\displaystyle \frac{{{\rm{\Pi }}}^{\alpha \sigma }{{\rm{\Pi }}}^{\beta \rho }\left({{\rm{\nabla }}}_{\sigma }{U}_{\rho }+{{\rm{\nabla }}}_{\rho }{U}_{\sigma }-\tfrac{2}{3}{g}_{\sigma \rho }{{\rm{\nabla }}}_{\nu }{U}^{\nu }\right)}{1+\tfrac{10}{27}\tfrac{1}{{({\rho }^{\prime }\kappa )}^{2}}{{\rm{\Pi }}}^{\mu \sigma }({{\rm{\nabla }}}_{\mu }{U}_{\nu })({{\rm{\nabla }}}_{\sigma }{U}^{\nu })}.\end{eqnarray} \tag{ 89 }$$

We plan to investigate the validity of this relation and explore its uses in analyzing simple shear flows in a future paper.

This work was supported in part by NASA Astrophysics Theory Program grants NNX14AB37G, NSF grant AST-1411879, and NASA's Fermi Guest Investigator Program. We thank Charles Gammie for useful comments, particularly for suggesting the use of these solutions as a test problem for current radiation–MHD codes.

Coughlin & Begelman (2015a) used the viscous equations of radiation hydrodynamics to analyze the two-stream boundary layer, which treats the transition between a relativistic jet and a surrounding envelope as confined to a thin (on the order of the square root of the mean free path of the photon) layer and considers the jet and the ambient medium as two distinct fluids. In this problem, the velocity of the fluid is primarily along the z-direction and the variation in the properties of the fluid occurs predominantly along the y-direction, meaning that the geometry of the problem is identical to that established in the preceding sections. Following the procedure outlined in Coughlin & Begelman (2015a), we define the four-velocity in the z-direction as

$$\begin{eqnarray}&&{\rm{\Gamma }}{v}_{z}={{\rm{\Gamma }}}_{j}{v}_{j}{f}_{\xi }(\xi )\end{eqnarray} \tag{ 90 }$$

and the comoving density as

$$\begin{eqnarray}&&{\rho }^{\prime }={\rho }_{0}^{\prime }g(\xi ).\end{eqnarray} \tag{ 91 }$$

In these definitions, v_j is the asymptotic velocity of the jet (and ${{\rm{\Gamma }}}_{j}={(1-{v}_{j}^{2})}^{-1/2}$ is its Lorentz factor), ${\rho }_{0}^{\prime }$ is the density of the ambient medium, f and g are functions of the self-similar variable

$$\begin{eqnarray}&&\xi ={\int }_{0}^{y}{\rho }^{\prime }\kappa \;{dy},\end{eqnarray} \tag{ 92 }$$

which we note is identical to our definition of ${\tau }_{y}^{\prime }$ (see Equation (3)), and subscripts ξ denote differentiation with respect to ξ (i.e., ${f}_{\xi }={df}/d\xi$, ${f}_{\xi \xi }={d}^{2}f/d{\xi }^{2}$, etc.). Inserting Equations (90) and (91) into the z-component of the momentum equation and the gas energy equation (and keeping only lowest-order terms in the boundary layer thickness; see Coughlin & Begelman (2015a) for details of this procedure) then yields the following two self-similar equations:

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

$$\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{ 94 }$$

In these equations $\chi ={e}_{0}^{\prime }/{\rho }_{0}^{\prime }$ is the ratio of the comoving radiation energy density in the jet to the comoving density of scatterers (they defined this quantity by μ, which we avoided for obvious reasons). By direct substitution, we can show that the assumption

$$\begin{eqnarray}&&{\rm{\Gamma }}{v}_{z}=\mathrm{sinh}(\mu \;\xi ),\end{eqnarray} \tag{ 95 }$$

where μ is an unspecified constant, exactly cancels the second term on the left-hand side of Equation (93) with the right-hand side. Furthermore, the first term in this same equation is proportional to $1/({{\rm{\Gamma }}}_{j}^{2}{v}_{j}^{2})$ so that, in the relativistic limit, this solution satisfies Equation (93) to order ${ \mathcal O }(1/{{\rm{\Gamma }}}_{j}^{2})$.

Inserting this solution into Equation (94), the equation for the density becomes

$$\begin{eqnarray}&&\displaystyle \frac{{g}_{\xi }}{{g}^{2}}=\displaystyle \frac{3}{2}{\mu }^{3}{{\rm{\Gamma }}}_{j}{v}_{j}\displaystyle \frac{1}{\mathrm{cosh}(\mu \;\xi )}.\end{eqnarray} \tag{ 96 }$$

This equation can be integrated and, redefining the density ${\rho }_{0}^{\prime }$ appearing in the equation as the density at $\xi =0$, so that $g(0)=1$, we find

$$\begin{eqnarray}&&g=\displaystyle \frac{1}{1-3{\mu }^{2}{{\rm{\Gamma }}}_{j}{v}_{j}\mathrm{arctan}\left(\mathrm{tanh},\left(\tfrac{\mu \xi }{2}\right)\right)}.\end{eqnarray} \tag{ 97 }$$

Since $\mathrm{tanh}(\pm \infty )=\pm 1$, this expression demonstrates that the comoving density varies between ${(1\pm 3\pi {\mu }^{2}{{\rm{\Gamma }}}_{j}{v}_{j}/4)}^{-1}$. Thus, as long as the inequality $\mu \lesssim 1/\sqrt{{{\rm{\Gamma }}}_{j}{v}_{j}}$ is satisfied, the comoving density is nearly constant in regions where the velocity varies as ${\rm{\Gamma }}{v}_{z}\simeq \mathrm{sinh}(\mu \xi )$, which is consistent with what we assumed based on the requirement of self-similarity in Section 2 of this paper.

Equations (95) and (97) cannot provide the solution throughout the entire two-stream boundary layer because they do not satisfy the boundary conditions (i.e., the velocity must approach zero as we proceed into the ambient medium and it must equal the jet velocity as we go into the jet). However, what we have shown is that, at any comoving point within the flow that is sufficiently far from the boundaries, the velocity and density profile do approach the self-similar forms that we assumed in our treatment here. Since the viscous equations of radiation hydrodynamics in the boundary layer limit are the same regardless of the boundary conditions, we likewise expect this velocity profile to appear in other relativistic, radiation–viscous problems (e.g., the free-streaming jet boundary layer; Coughlin & Begelman 2015b). We thus expect that the self-similar velocity profile ${\rm{\Gamma }}v=\mathrm{sinh}(\mu {\tau }_{y}^{\prime })$ is an inherent feature of most relativistic shear flows.