THE DYNAMICS, APPEARANCE, AND DEMOGRAPHICS OF RELATIVISTIC JETS TRIGGERED BY TIDAL DISRUPTION OF STARS IN QUIESCENT SUPERMASSIVE BLACK HOLES

To calculate the rate of mass return to the black hole, $\dot{M}_{\rm h}$, after tidal disruption, we use the same formalism as described in Guillochon et al. (2011). To minimize numerical artifacts in the solution (Guillochon et al. 2011), the simulations are performed in the rest frame of the star. The relative positions of the star and the black hole (taken as a point mass) are explicitly tracked over the course of the simulation, with the star having initially a null velocity (at infinity). We take the star to be a 1 M_☉ star described by a polytropic index n = 3 and an adiabatic equation of state (P∝ρ^γ) with the adiabatic index γ = 5/3. Our disruption simulations are performed using the FLASH hydrodynamics framework (Fryxell et al. 2000), which includes an adaptive mesh refinement scheme that permits the wide range of scales necessary to resolve the star and the debris streams simultaneously. During the disruption, the debris tails are adaptively refined based on their density relative to the maximum density in the simulation, ρ_max, at each time step, with a cutoff density for the lowest refinement level of 10⁻¹⁹ times the maximum density. All matter with ρ > 10⁻³ρ_max is refined to the highest level. The disruption is performed in a box that is several orders of magnitude larger than the initial star's radius of r_* = R_☉, chosen to facilitate tracking of the tidal debris over long enough timescales for hydrodynamical effects to no longer play a major role in determining the distribution of mass as a function of binding energy.

To calculate the rate of mass return to the black hole $\dot{M}_{\rm h}$ as a function of the black hole mass M_h, we performed a single simulation of the disruption of a 1 M_☉ star by a 10⁶ M_☉ black hole with impact parameter β ≡ r_t/r_p = 2, where r_t and r_p are the tidal and the periastron radius, respectively (see Figure 1). This encounter is deep enough to leave no surviving core (Guillochon & Ramirez-Ruiz 2012). To accurately model this disruption, we perform a simulation using a cubical volume with a basic grid of 8³ cells (along each Cartesian coordinate axis), with a width of 4 × 10¹⁴ cm, and 16 levels of refinement, resulting in the initial diameter (=2 R_☉) of the star being resolved by 90 grid cells. The shape of $\dot{M}_{\rm h}$ in the Newtonian approximation for full disruptions is only affected by the degree of symmetry in the tides (Guillochon et al. 2011) and the star's initial density profile (Lodato et al. 2009; Ramirez-Ruiz & Rosswog 2009). As shown in Guillochon et al. (2011), tides are symmetric for black holes with M_h ≳ 10⁶ M_☉, with the maximum force difference between the near- and far-side of the star being 6% for a β = 2 encounter with a M_h = 10⁶ M_☉, and even less for black holes of larger mass. Assuming that the shape of $\dot{M}_{\rm h}$ does not vary for black holes of mass ⩾10⁶ M_☉, $\dot{M}_{\rm h}(M_{\rm h})$ can be related to our benchmark simulation via

$$\begin{equation} \dot{M}_{\rm h}(M_{\rm h}, t) = \dot{M}_{\rm h} (10^{6} \,M_\odot , t^{\prime })\left(\frac{M_{\rm h}}{10^{6}\,M_\odot }\right)^{-1/2} \end{equation} \tag{ 1 }$$

$$\begin{equation} \hspace{-1.5pc}t^{\prime } = t\left(\frac{M_{\rm h}}{10^{6}\,M_\odot }\right)^{-1/2}, \end{equation} \tag{ 2 }$$

where t and t' are the times since disruption for a black hole of mass M_h and our benchmark simulation, respectively (Evans & Kochanek 1989). Stars of different masses and ages can have different degrees of central concentrations (Tout et al. 1996). Lodato et al. (2009) showed that, for full disruptions, the star's initial density profile primarily affects $\dot{M}_{\rm h}$ at early times, but that all reasonable initial stellar profiles eventually lead to a power-law decay where $\dot{M}_{\rm h} \propto t^{-5/3}$ at late times. Therefore, we assume that a typical disruption of a 1 M_☉ by a black hole of mass M_h ⩾ 10⁶ M_☉ can be reasonably modeled by a single simulation assuming a star with a similar structure to the Sun using the above scalings, which yields the $\dot{M}_{\rm h}$ shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The returning gas on bound orbits from the tidally disrupted star does not immediately produce a flare of activity from the black hole. First material must enter quasi-circular orbits and form an accretion torus (Ramirez-Ruiz & Rosswog 2009). Only then will viscous effects release enough binding energy to power a flare. Once the torus is formed, it will evolve under the influence of viscosity and radiative cooling, although the viscosity would have to be implausibly low (i.e., the usual viscosity dissipation time t_ν for a thick disk would be inversely proportional to the viscosity parameter α) for the bulk of the mass to be stored for longer than t_peak in a reservoir at r ≈ r_t, where t_peak is the time corresponding to the peak in the accretion rate and r_t is the tidal radius. After t_ν ≪ t_peak the mass accretion rate would continue to fade as $\dot{M}_{\rm h}$.

Little is known about the relation between jet production in supermassive accreting black holes and the state of their constituent accretion disks. However, an association between hot thick accretion flow with the strongest jets, as observed in binary black holes, is expected (e.g., Krolik & Piran 2012). It is the poloidal field protruding from the disk that is thought to drive the jet and because its strength increases with disk thickness, jets are expected to be stronger in thick disks than in thin ones (Meier et al. 2001). In GX 339-4, for example, a jet is produced when the X-ray source is in the low/hard state. In this state, the disk temperature is T ∼ 10⁹ K, suggestive of a thick disk. On the other hand, when the source enters the high/soft state (with a thin 10⁷ K disk) the radio jet emission drops to a level at least tens of times weaker (Fender et al. 1999). In this paper, we make the assumption that jets will be produced only when the accretion disk is geometrically thick and hot.

It is clear from Figure 3 that most of the debris would be fed to the black hole much more rapidly than it could be accepted if the radiative efficiency was high. As the disk material advects onto black hole, we thus assume a relativistic jet with $L_{\rm j} (t)\propto \dot{M}_{\rm h} (t)$ can be powered for as long as $\dot{M}_{\rm h} (t)\gtrsim \dot{M}_{\rm Edd}$.

Zoom In Zoom Out Reset image size

The structure of the optically thick, high-entropy envelope formed as a result of the tidal disruption of a star by a massive black hole has been described by Loeb & Ulmer (1997). Their results, used as initial conditions in Section 4, are briefly described here. During the tidal disruption process, about half of the stellar material escapes on hyperbolic orbits with speeds ∼9500(M_h/10⁶ M_☉)^1/6 km s⁻¹, while the rest falls back onto the black hole. The bound gas, after pericenter passage, is on orbits which collide with the infalling stream near the original orbital plane at apocenter, giving rise to a shock which redistributes angular momentum (Rosswog et al. 2009; Ramirez-Ruiz & Rosswog 2009). The debris raining down would, after little more than its free-fall time, settle into a disk surrounded by a radiation-dominated envelope, whose inner radius is

$$\begin{equation} r_{\rm t} = r_{\ast } \left(\frac{M_{\rm h}}{M_{\ast }} \right)^{1/3} \approx 2.15 \times 10^{13} \left(\frac{M_{\rm h}}{10^6 \,M_\odot } \right)^{1/3}\,\mathrm{cm} , \end{equation} \tag{ 3 }$$

where r_* ≈ R_☉ and M_* ≈ M_☉.

The structure of such an envelope is simplified by the fact that a fully ionized gas dominated by Thomson opacity and radiation pressure tends to approach a uniform entropy state. Assuming hydrostatic equilibrium and an equation of state dominated by radiation pressure (p∝ρ^4/3), the density stratification in the optically thick envelope can be written as

$$\begin{equation} \rho _{\tau } = \frac{f M_\ast }{4 \pi \ln (r_{\tau }/r_{\rm t}) r^3}, \end{equation} \tag{ 4 }$$

where f is the fraction of stellar material in the envelope and r_τ is the radius at which the envelope becomes optically thin, which for electron scattering opacity is given by

$$\begin{equation} r_{\tau } \approx 1.7 \times 10^{15} \left(\frac{f M_{\ast }}{0.5 \,M_\odot }\right)^{1/2}\,\mathrm{cm.} \end{equation} \tag{ 5 }$$

The density distribution derived in Equation (4) is probably valid only at late times (t ∼ t_peak), as the returning gas does not immediately produce a thick disk of radiation-dominated gas. First material must enter quasi-circular orbits. The bound orbits are initially very eccentric and the range of orbital periods is large. The orbital semimajor axis of the most tightly bound debris is

$$\begin{equation} a \approx 10^3 \left({M_{\rm h}\over 10^6 \,M_\odot }\right)^{-1/3}\left({r_\ast \over R_\odot } \right)\left({M_\ast \over M_\odot }\right)^{-2/3} r_{\rm g}, \end{equation} \tag{ 6 }$$

and the period is

$$\begin{equation} t_{\rm a} \approx 7.3 \left({a \over 10^3 r_{\rm g}}\right)^{3/2} \left({M_{\rm h}\over 10^6 \,M_\odot }\right)^{-1/2}\,{ \rm days}, \end{equation} \tag{ 7 }$$

where r_g ≈ 1.5 × 10¹¹(M_h/10⁶ M_☉) cm is the gravitational radius. As a result of internal dissipation due to high viscosity or shocks, the debris raining down would, after one or two orbital periods (t_a ≪ t_peak), form a highly elliptical disk surrounded by an extended atmosphere with a wide spread in apocentric distances, whose progressively increasing vertical scale is likely to be smaller than the one predicted by Equation (5). For a 10⁷ M_☉ black hole, t_a at a ≈ r_τ is approximately 161 days while t_peak ≈ 89 days.

Aside from our own Galactic center and M31, we have no direct measurement of the gas content at sub-parsec scales near massive black holes. The stellar density, however, is more well known. After all, if the stars were not closely packed near the center of the galaxy, we would not have evidence for central black holes within quiescent galaxies. These same stars should provide mass winds, whose strength depends on the concentration of stars enclosed in the black hole's sphere of influence as well as the rates and velocity of the mass injection of those winds.

To determine the gas structures surrounding quiescent supermassive black holes, we follow the formalism developed by Quataert (2004), who modeled the distribution of hot gas around the central parsec of the galactic center under the assumption of spherical symmetry and an adiabatic equation of state. The one-dimensional hydrodynamical equations are also solved using FLASH (Fryxell et al. 2000), following the method described in Naiman et al. (2012). The winds from the closely packed stellar members are assumed to shock and thermalize such that density and energy contributions can be treated as source terms in the hydrodynamical equations. In spherical symmetry, these equations can be written as (Holzer & Axford 1970)

$$\begin{equation} \frac{1}{r^2} \frac{d}{dr} (\rho _{\rm a} v r^2 ) = q(r) \end{equation} \tag{ 8 }$$

$$\begin{equation} \rho _{\rm a} v \frac{dv}{dr} = - \frac{d p}{dr} - \frac{G M_{\rm h} \rho _{\rm a}}{r^2} - q(r) v \end{equation} \tag{ 9 }$$

$$\begin{equation} \frac{1}{r^2} \frac{d}{dr} \left[ \rho _{\rm a} v r^2 \left(\frac{v^2}{2} + \frac{c_{\rm s}^2}{\gamma -1} \right)\right] + \frac{ \rho _{\rm a} v G M_{\rm h}}{r^2} = \frac{q(r) v_{\rm w}^2}{\gamma (\gamma - 1)} , \end{equation} \tag{ 10 }$$

where the velocity, v, and density of the ambient medium, ρ_a, depend solely on r, and c_s is the sound speed of the gas.

Within the central parsec of the galactic nuclei, gas is assumed to be supplied by winds which originate from massive stars. In the Milky Way's galactic center, these stars include blue super giants each with mass loss rates $\dot{M}_w \approx 10^{-4} \,M_\odot \,{\rm yr}^{-1}$ and wind speeds of v_w ≈ 600–1000 km s⁻¹ (e.g., Najarro et al. 1997). Here we use the term q(r) to quantify the total rate of mass injection from the stellar winds:

$$\begin{equation} \Sigma \dot{M}_{\rm w} = \int 4 \pi r^2 q(r) dr. \end{equation} \tag{ 11 }$$

Following Quataert (2004) we set q(r)∝r^−η with η = 0, 2, 3 for r  [r_m, r_M] and q(r) = 0 otherwise. Different values of η correspond to different mass injection distributions. A value of η = 0, for example, describes mass that is injected preferentially at large radii while η = 3 corresponds to equal mass injection for all radii within [r_m, r_M]. To uniquely specify q(r), we must determine $\Sigma \dot{M}_{\rm w}$, η, and [r_m, r_M]. For simplicity, we assume the central star cluster properties from Quataert can be scaled with the central black hole's mass, so that [r_m, r_M]∝M^1/3_h and $\Sigma \dot{M}_{\rm w} \propto \Sigma N_\ast \propto M_{\rm h}$.

Figure 5 shows the results of our simulations for different black hole masses. After many sound crossing times the flow of gas from stellar winds settles into a steady state. Far away from the cluster, as expected, the flow is driven out by the aggregate influence of the stellar winds while interior to the cluster the gas is captured and accreted onto the central black hole. The stagnation radius, the boundary where the flow is divided between inflowing and outflowing, increases with the mass of the central black hole. As can be seen from Figures 5 and 6, the density stratification is strongly dependent on the black hole mass but nearly independent on the power law assumed for the distribution of the mass injection (i.e., η).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The density distribution at small radii can be roughly described by ρ_a∝r^−k with k ≈ 1, 1.4, 1.5 for M_h = 10⁶, 10⁷, and 10⁸ M_☉, respectively. In the outer regions, the density is seen to rapidly converge to a k = 2 wind profile. Realistically, we would expect the ambient medium around the supermassive black hole to have a complex multiphase structure as inferred from models of the line-emitting region (Barai et al. 2012).

A simple analytical argument can be used to understand the dynamical evolution of jets. Let us consider a cylindrically symmetric, relativistic jet moving through a stratified medium. As long as the flow is continuously injected, the head of the jet (or working surface) will have a double shock structure, composed by a forward and a reverse shock where, in the frame of reference where the contact discontinuity is at rest, the ambient medium and the jet material respectively are decelerated and heated, transforming kinetic into thermal energy.

The speed v_h with which the head of the channel advances into the surrounding medium is obtained by balancing the momentum flux in the shocked jet material against that of the shocked surrounding medium, measured in the frame comoving with the advancing head (e.g., Begelman & Cioffi 1989; Ramirez-Ruiz et al. 2002; Matzner 2003; Bromberg et al. 2011; Lazzati et al. 2012; Metzger et al. 2012):

$$\begin{equation} \rho _{\rm j} h_{\rm j} \Gamma _{\rm j}^2\Gamma _{\rm h}^2 (v_{\rm j}-v_{\rm h})^2 + p_{\rm j} = \rho _{\rm a} h_{\rm a} \Gamma _{\rm h}^2v_{\rm h}^2 + p_{\rm a}{,} \end{equation} \tag{ 12 }$$

where ρ_{j, a}, h_{j, a}, Γ_{j, a}, and p_{j, a} are the mass density, specific enthalpy, Lorentz factor, and pressure of the jet and ambient medium, respectively, and v_{j, h} is the velocity of the jet and the advancing head (where u = Γβ = Γv/c is the proper fluid velocity and h_a ≃ 1). If v_h is supersonic with respect to the ambient gas, p_{a, j} can be neglected and Equation (12) gives

$$\begin{equation} v_{\rm h}={v_{\rm j}}\left[1+\left(\frac{\rho _{\rm a}}{\rho _{\rm j} h_{\rm j}\Gamma _{\rm j}^2}\right)^{1/2}\right]^{-1}. \end{equation} \tag{ 13 }$$

The resulting flow pattern will depend crucially upon the jet Mach number and the density ratio between the jet and the given ambient medium. When the jet density significantly exceeds the ambient density, v_h ≈ v_j.

A relativistic jet emerging from the central region of the galaxy will pass through a region of steadily decreasing ambient density for up to several tens of parsecs (Figure 4). Even if the jet is confined, this decrement in density will result in an increase in cross section, and the degree of collimation θ = s/r will either decrease or increase, depending on whether the size of the evacuated channel s increases more or less rapidly than r. This indicates that the dynamics of the expanding jet is expected to be modified by changes in collimation and cannot be properly captured by spherically symmetric solutions. These difficulties motivate consideration of the hydrodynamical confinement in axisymmetric numerical calculations, to which we now turn our attention.

The propagation of relativistic jets triggered by the tidal disruption of stars in quiescent supermassive black holes is studied here using the adaptive mesh refinement code Mezcal. The code solves the special relativistic hydrodynamics (SRHD) equations in two-dimensional spherical (polar) coordinates. The Mezcal code integrates the SRHD equations by using a second-order upwind scheme. The equation of state, relating enthalpy to pressure and density, assumes an adiabatic index γ = 4/3. A detailed description of the SRHD version of the Mezcal code is presented in De Colle et al. (2012a) together with a series of standard tests.⁴

Common to all calculations is the initiation of a conical jet with an initial opening angle θ₀ and uniform velocity u₀ = Γ₀β₀, with matter injected along the symmetry axis at the inner boundary of the computational domain r = r₀. The luminosity of the jet is assumed here to follow the mass accretion rate, $L_{\rm j} (t) =\epsilon _{\rm j} c^2 \dot{M}_{\rm h} (t)$. This relation together with mass conservation determines the density of the jet, ρ_j(r₀), at the injection boundary. Given the large range in scales, the propagation of the jet in the optically thick envelope (r ⩽ r_τ) and its subsequent expansion through the medium shaped by the interaction of massive stellar winds (r ⩽ r_M ≈ 10³r_τ) are studied separately.

The existence of a steady, spherical, optically thick envelope around the black hole should be regarded as an extreme assumption because the extension and structure will be modulated at a variable rate by fallback. The envelope starts to form when the most tightly bound debris falls back. Our stellar disruption simulation (Section 2) shows that the first material returns at a time t_a, with an infall rate of about $\dot{M}_{\rm a}\approx 10^{-4} \,M_\odot \,{\rm yr}^{-1}\le \dot{M}_{\rm Edd}$ (Figure 2). The infall rates are expected to increase, relatively steadily, for at least tens of orbital periods, before the majority of the bound material rains down at t ≈ t_peak. Once the envelope is formed, it will evolve under the influence of viscosity, radiative cooling, and time-dependent mass injection from both the jetted outflow and the angular momentum redistribution region.

To study the evolution of the relativistic jet in the optically thick region, for simplicity we assume the envelope to be fully formed (t ∼ t_peak) and accurately described by the steady, spherical hydrostatic solutions presented in Section 3.1. In addition, since the jet head breaks out of the envelope after a time r_τ/c ≪ t_peak, the jet power is assumed to be steady during this phase with $L_{\rm j} \approx \epsilon _{\rm j} \dot{M}_{\rm peak} c^2$. While in principle the fall-back matter may interrupt the jet propagation (if the ram pressure of the infalling material is larger than the jet kinetic energy), this is unlikely to happen, as the accretion on the central black hole happens mainly on the equatorial plane.

Snapshots of detailed hydrodynamic simulations of the evolution of the jet through the optically thick envelope are presented in Figure 7, where the density maps of the expanding jet at various times after t_peak in its hydrodynamical evolution are plotted. Initially, the low-density jet is unable to move the envelope material at a speed comparable to its own and thus is abruptly decelerated. Most of the energy output during that period is deposited into a cocoon surrounding the jet, in which the energy supplied by the jet exceeds that imparted to the swept-up envelope material by a factor β⁻¹_h. In the cocoon region, the jet and ambient medium remain separated by a contact discontinuity where shearing instabilities are prominent. The cocoon region exhibits two important dynamical effects (e.g., Bromberg et al. 2011). First, it forms a weak shock that moves laterally at the relativistic internal sound speed, and second, it acts to partially collimate the jet itself. As the jet expands further into the envelope, the drastic density drop permits the jet head to accelerate to velocities close to the speed of light (β_h ≈ 1). Thus, if it is a general property that the jet becomes relativistic near the boundary of the envelope, the outer edge of the envelope is reached in a crossing time ≈r_τ/(cβ_hΓ²_h) as measured by an observer along the line of sight.

Zoom In Zoom Out Reset image size

After the jet has broken free from the envelope region, the fraction of relativistic plasma injected into the cocoon will be much reduced as the jet accelerates and β_h ≈ 1. The amount of energy that accumulated in the cocoon while the jet was advancing sub-relativistically is

$$\begin{equation} E_{\rm c}\approx 10^{50}\bar{\beta }^{-1}_{\rm h} \left({r_{\tau } \over 10^{15}\,{\rm cm}}\right)\left({L_{\rm j}\over 10^{46}\,{\rm erg\,s^{-1}}}\right)\,{\rm erg}, \end{equation} \tag{ 14 }$$

where $r_{\tau }/\bar{\beta }_{\rm h}$ is the envelope traversal time and $\bar{\beta }_{\rm h}$ is the average speed of the jet head. The energy accumulated in this phase is thus larger than the binding energy of the envelope.⁵ At the radius ≈r_τ where the head of the jet starts to advance relativistically, the volume of the region incorporated into the cocoon is related to both the jet and cocoon expansion velocities by V_c ≈ (π/3)r³_τ(β_c/β_h)². Unlike the jet, this cocoon material does not have a relativistic outward motion, although it has a relativistic internal sound speed. At first an elongated bubble (since pressure balance may never be reached within a radiation-dominated isentropic atmosphere) will be inflated, which can expand most rapidly along the rotation axis and will eventually unbind the envelope (Figure 8). Even under the extreme assumption that the envelope is fully formed at t ≈ t_peak, the energy deposited by the jet is expected to eject most of the envelope material in the optically thick region. Thus, passage through this region cannot significantly alter the jet launching conditions or enhance collimation.

Zoom In Zoom Out Reset image size

4.3.1. Internal Shocks

Instead of assuming the jet to be a steady outflow, here we suppose that it is irregular on timescales much shorter than t_peak. For instance, if the Lorentz factor in the outflowing collimated ejecta varied by a factor of more than two, then the shocks that developed when fast material overtook slower material would be internally relativistic (Rees & Meszaros 1994). Dissipation would then take place whenever internal shocks developed in the ejecta, which can then be reconverted into energetic particles and radiation (the jet may also lose energy as it propagates through the photon field of the accretion disk). Internal shocks generated either as a consequence of fluctuations at the source (e.g., Kobayashi et al. 1997; Ramirez-Ruiz & Lloyd-Ronning 2002) or arising from the development of large-amplitude instabilities provide an attractive explanation for the large-scale variability seen in Swift 1644+57 (Bloom et al. 2011; Burrows et al. 2011).

To illustrate the basic idea, suppose that two blobs of equal mass, but with different Lorentz factors Γ_i and Γ_j (with Γ_i > Γ_j ≫ 1) are ejected at times t₁ and t₂, where t₂ − t₁ = δt. In the case of highly relativistic ejecta, the shock develops after a distance of order

$$\begin{equation} r_{\rm \iota } \approx c \delta t \frac{2 \Gamma _i^2\Gamma _j^2}{\Gamma _i^2-\Gamma _j^2} \approx 6\times 10^{15} \left({\delta t \over 10^3\,{\rm s}} \right) \left({\Gamma _j\over 10} \right)^2\,\rm cm. \end{equation} \tag{ 15 }$$

The reconversion of bulk energy can be very efficient: when the two blobs share their momentum, they move with $\Gamma _{ij}=\sqrt{\Gamma _i\Gamma _j}$, so the fraction of the energy dissipated is

$$\begin{equation} \epsilon _{\iota }={\Gamma _i +\Gamma _j -2\sqrt{\Gamma _i\Gamma _j} \over \Gamma _i +\Gamma _j}. \end{equation} \tag{ 16 }$$

High efficiency does not, therefore, require an impact on matter at rest; all that is needed is that the relative motions in the comoving frame be relativistic.

The evolution of a non-steady relativistic jet with randomly varying Lorentz factor Γ (between 2 and 20) over a range of timescales 10² s ≲ δt ≲ 10³ s is shown in Figure 9. The average parameters of the jet are chosen so that they are similar to those displayed in Figure 7. Compared to Figure 7, significant structural differences along the evacuated channel appear when faster material catches up with slower material and a strong shock forms. An internal shock in the relativistic jet will move with a Lorentz factor of up to Γ_ij ≈ 10, and the emitting material behind the shock will be subject to a large Doppler boost. To be most effective, dissipation must occur when the envelope is optically thin, i.e., r_ι ≳ r_τ (Figure 9).

Zoom In Zoom Out Reset image size

The deceleration of the working surface allows slower ejecta to catch up with the head of the jet, replenishing and re-energizing the reverse shock and boosting the momentum in the working surface. Since the efficiency of converting bulk motion to radiation in an internal shock depends on the difference between Γ's, the fact that the working surface is decelerating (β_h ≲ 1) implies that the efficiency will be much higher than in the standard GRB-like case, where Γ is assumed to fluctuate by a factor of a few within the jet and among blobs on a typical timescale δt ≳ r_g/c ≈ 50(M_h/10⁷ M_☉) s.

We now turn our attention to the evolution of the jet as it expands through the large-scale environment of the galactic nuclei where the density stratification is thought to be regulated by the interaction of winds from the surrounding massive stars. The deceleration of the jet within this region gives rise to the non-thermal, long-lived radio emission, such as that observed in Sw 1644+57 (Zauderer et al. 2011; Berger et al. 2012). The density profile within this region is taken from the simulations of the interacting stellar winds described in Section 3.2 and shown in Figures 5 and 6. The power of the jet is assumed to follow the mass feeding rate, $L_{\rm j} (t) =\epsilon _{\rm j} c^2 \dot{M}_{\rm h} (t)$, with $\dot{M}_{\rm h} (t)$ taken from the results of the tidal disruption simulations presented in Section 2 and Figure 2.

Little is known about the relation between jet production in the more massive accreting black holes and the state of their constituent accretion disks. However, a similar association of hot thick accretion flows with the strongest jets as seen in binary black holes is expected. Motivated by this, we assume that strong relativistic jets preceding tidal disruption are only triggered when $\dot{M}_{\rm h} (t)\gtrsim \dot{M}_{\rm Edd}$. The properties of these jets are thus severely constrained by the mass feeding rate history. The emitted radiation is an observable diagnostic which provides constraints on the processes occurring at their point of origin.

The propagation of a jet injected at the onset of the super-Eddington accretion phase over successive decades in radius ranging from 10⁻² pc to 1 pc is shown in Figure 10 for M_h = 10⁷ M_☉. Two illustrative cases are depicted. In case a, the jet power increases with time to achieve a maximum at t = t_peak ≈ 0.27 yr, it then subsequently decreases as ∼t^−5/3 until it is finally halted at the time when $\dot{M}_{\rm h} (t) < \dot{M}_{\rm Edd} \approx 0.71$ yr. In case b, the injection history is the same until the jet power is prematurely halted at t = t_peak. Figure 11 shows the corresponding schematic world lines of the bow shock evolution for the two illustrative cases depicted in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As seen in Figure 10, the jet evacuates a channel out to some location where it impinges on the surrounding medium. A continuous flow of relativistic fluid emanating from the nucleus supplies this region with mass, momentum, and energy. The jet, unable to move the surrounding material at a speed comparable to its own, is decelerated and, as a result, most of the energy output during that period is deposited into a cocoon surrounding the jet. As a first approximation, evolution of the cocoon's shape is governed both by the advance of the head and by the cocoon's own pressure-driven sideways expansion into the surrounding medium. Heavy jets therefore propagate almost ballistically and are naturally surrounded by small cocoons. Light jets have large cocoons and, as seen in Figure 10, can be confined by the pressure of the shocked material.

The cocoon will effectively collimate the jet (Bromberg et al. 2011) as long as

$$\begin{equation} \rho _{\rm j} h_{\rm j} \Gamma _{\rm j}^2/\rho _{\rm a} < \theta _0^{-4/3}{,} \end{equation} \tag{ 17 }$$

where ρ_jh_jΓ²_j/ρ_a gives the ratio between the jet's energy density and the rest-mass energy density of the surrounding medium at the location of the head. Collimation is then seen to increase with decreasing k for k ⩽ 2. When k = 2 the jet's head velocity is constant and consequently most of the energy flowing into the cocoon cannot effectively counterbalance the jet's expansion. The density profile ρ_a∝r^−k in the stellar wind mass injection region (region II in Figure 4) varies from 1.4 ≲ k ≲ 2 with increasing r for a black hole with M_h = 10⁷ M_☉ and, as a result, the cocoon's pressure is expected to collimate the jet and decrease its opening angle. This is clearly illustrated in Figure 11 in which the results of the simulations are directly compared with the evolution of the working surface as predicted by Equation (13). In a jetted source like Sw 1644+57, the expansion of the working surface within the external environment on scales r ⩽ 1 pc is clearly incompatible with that predicted on the assumption of a constant jet opening angle.

If the jet were suddenly to turn off (in case a this occurs at t = t_Edd while in case b at t = t_peak) it would be preceded by the collapse of the jet channel which, suddenly evacuated, would be filled in on the transverse sound crossing time by whatever material happened to surround it. This ambient or cocoon material could be driven into the jet channel at the internal sound speed (Figure 10). After the evacuated channel collapses, a cylindrical rarefaction wave will propagate along the jet channel, reaching the head at t = t_r ≫ t_peak when the jet has traveled a much larger distance (r_r) than when it was initially turned off. This is followed by a series of compression and rarefaction waves of decreasing amplitude. In case a, the inclusion of the ∝t^−5/3 injection phase produces only a small change in t_r. This is clearly illustrated in Figure 12, where we plot the total internal energy available at the working surface. After the rarefaction wave reaches the jet head, which is observed to take place at comparable times in both of the cases depicted in Figure 10, the available internal energy drops although it does so less abruptly in case a where the jet continues to be powered until t_Edd. The intensity of the radiation from the radio-emitting electrons, which have recently emerged from the hot spot rapidly decreases after a time t_r ≫ t_peak, which corresponds to a time ≈r_r/(cβ_hΓ²_h) as measured by an observer along the line of sight.

Zoom In Zoom Out Reset image size

When the jet becomes impulsive at t ≳ t_r, the necessary condition for the jet to remain in pressure balance with its surroundings is no longer satisfied and strong shocks are driven into the jet. As a result, the contact discontinuity and the forward shock are abruptly decelerated. This deceleration allows ejecta to catch up and pass through a reverse shock just inside the contact discontinuity, increasing the dissipated internal energy at the shock front (Figure 13) and in principle giving rise to a longer-lived afterglow than that predicted by the communication delay between the working surface and the base of the jet. As the jet become non-relativistic, the contribution to the total luminosity of the jet coming from the non-relativistic material (both at the jet head and the cocoon) will become prominent. The late radio observations can, if coupled with a detailed modeling, therefore be used to infer the deceleration of the system. The study of the detailed radiation properties of the system is left to a future study.

Zoom In Zoom Out Reset image size

Quiescent black holes appear to be capable of producing directed flows of relativistic matter. Given the twin requirements of super-Eddington luminosity and short flaring timescales, the currently favored models involve the disruption of a star with the subsequently energy dissipation from accreted debris providing the required luminosity. Here we examine, within this model, the consequences of assuming a relation between jet production and the state of their constituent accretion disks. In those sources that display jets, production of the jet should only occur when the disk is in a hard accretion state and not when the disk emission is soft. Testing this prediction in quiescent galactic nuclei, as argued in Section 5.3, requires statistical studies because the typical jet production duty cycle time is expected to be thousands of years. Confirmation of this effect will provide additional evidence that accretion states similar to those in binary black holes also should exist in supermassive black holes. As we show in Section 5.2, many of the observed properties of Swift 1644+57 can be understood as resulting from a hot accretion flow, driven into a previously quiescent supermassive black hole, which is then subsequently collimated into a pair of anti-parallel jets. If we were to venture a general classification scheme for tidal disruption jets, we would obviously expect the black hole mass, the rate at which the gas is supplied to the black hole, the angular momentum of the black hole, the flow velocity, and the orientation relative to our line of sight to be the essential parameters.

Most of the radiation we receive is reprocessed by matter quite distant from the black hole. The fact that jets are detectable at all indicates that some re-randomization of kinetic energy and re-acceleration of particles must be occurring along the jet length to counteract radiation and adiabatic losses. Under the assumption that these prime movers are able to form collimated, relativistic outflows, the observable effects of the two major radiating regions, the working surface and internal shocks, are studied. This leads to a unified picture where internal shocks (and/or external Compton at the base of the jet) can provide most of the variable, high photon energy luminosity, and where the region of greatest radio emissivity is located at the head of the jet.

Hot radio spots are naturally interpreted as the working surface at the end of a jet. The jet is decelerated at a strong collisionless shock where particle acceleration and field amplification can occur. The relativistic electrons and magnetic field pressure are balanced by the ram pressure of the ambient medium. Because the source moves relativistically, aberration of light must be accounted for when calculating the observed radio light curve and the local direction of polarization. And as the radio hot spot is unlikely to be resolved observationally, we can only measure the average polarization over the whole image. As a result, the breaking of symmetry of the emitting region around our line of sight is required in order to produce a net polarization. Statistical studies over a sample of tidal disruption jets, or time-resolved polarimetry of different emission episodes within a single event, may teach us about the dominant emission mechanism, the jet structure, or the magnetic field configuration within the hot spot.

It is obvious from the discussions in this paper that the dynamics of tidal disruption jets are complex, especially because of rich interactions between the jet and the surrounding medium. Confirmation was provided of the important notion that the visibility of jets, in particular the emission emanating from the head of the jet, is determined largely by the luminosity and velocity injection history of the transient jet as well as the properties of the surrounding medium. Axisymmetric hydrodynamical calculations of these jets show that the resulting dynamics are different from those predicted by one-dimensional models; interaction with the surrounding medium cannot only enhance collimation but also mediate the communication delay between the base and the head of the jet. Even in the simplest case of a jet whose thrust directly scales with $\dot{M_{\rm h}}$, complex behavior with multiple possible transitions in the observable part of the radio afterglow lifetime may be seen. The eventual resulting afterglow light curve depends fairly strongly on the properties of $\dot{M_{\rm h}}$, especially the impact parameter β and mass of the disrupted star and the black hole. There are advantages and disadvantages due to this dependence on $\dot{M_{\rm h}}$, β, M_☉, and M_h. On the negative side, these dependencies imply that one cannot be too specific about the times at which we expect to see transitions in the observed emission. On the positive side, if and when we do see these transitions, they can be fairly constraining on the properties of the system.

The X-ray and radio flux from Swift 1644+57 has been interpreted as being emitted by the jet produced from the super-Eddington accretion of material resulting from the disruption of a star by a supermassive black hole. For the first ∼30 days since the flare triggered the Swift/BAT instrument, the object exhibited tremendous variability on ≲ 10⁴ s timescales (Burrows et al. 2011), with isotropic equivalent flare luminosities >10⁴⁸ erg s⁻¹ in the soft X-rays alone. Beyond this first month of violent activity, the degree of variability decreased, and appeared to follow the predicted t^−5/3 power law associated with the late-time evolution of fallback of tidal stream onto the black hole (Berger et al. 2012).

The initial month-long X-ray light curve plateau, punctuated by periods of extreme variability, has been presumed to be associated with the phase of the event in which the accretion rate exceeds the Eddington limit (Berger et al. 2012; Metzger et al. 2012). However, as the Eddington limit is ∼100 times smaller than the peak accretion rate for a complete disruption by a 10⁶ M_☉ black hole, only a disruption by a more massive black hole or a partial disruption, for which the peak accretion rates are comparable to the Eddington rate, are capable of producing such a short period of super-Eddington accretion (Figures 2 and 3). In other words, we would expect the X-ray plateau to last longer for a tidal disruption on a central black hole with M_h ∼ 10⁶. But while a partial disruption or a disruption by a more massive black hole can yield a plateau of the appropriate duration, these events peak at a much later time relative to the time of disruption, resulting in a decay slope that yields a poor match to the observed X-ray light curve.

A proposed solution to this problem is to allow the star to penetrate more deeply at pericenter (Cannizzo et al. 2011), for which analytical approximations in which the binding energy distribution across the star is presumed to be "frozen in" at pericenter. Evans & Kochanek (1989) and Lodato et al. (2009) can produce shorter flare times. These shorter flares could reproduce both the month-long plateau and the observed decay rate. However, numerical simulations by Guillochon & Ramirez-Ruiz (2012) have shown that an increase of β beyond the value for which the star is completely destroyed do not result in faster transients, with the time of peak actually tending to slightly larger values with increasing β for complete disruptions. For very deep encounters in which general relativistic effects are important, Laguna et al. (1993) show that while t_peak tends to smaller values, it does not scale as sharply as predicted by Evans & Kochanek (1989) and Lodato et al. (2009), and the time at which a β = 5 and a β = 10 encounter cross the Eddington limit is almost identical. Therefore, a deep encounter is incapable of reproducing the duration of the two phases exhibited by the X-ray light curve.

We suggest an alternative model in which the flare is a standard full disruption of a main-sequence star, but that the triggering of the BAT instrument and the month-long plateau are actually artifacts of variability driven at a timescale comparable to the orbital period at the tidal radius, $\tau _{\rm t} = 2\pi \sqrt{r_{\rm t}/{\it GM}_{\rm h}}$. In this model, $\dot{M}_{\rm h}$ exceeds Eddington even during the observed t^−5/3 decay phase. Variability on this timescale could be driven by the interaction of the material that returns to pericenter with previously accreted material within the accretion disk. Hydrodynamical studies of the fallback stream have shown that the material that returns to pericenter after a disruption passes through a small nozzle region where it is compressed violently (Ramirez-Ruiz & Rosswog 2009; Rosswog et al. 2009), resulting in a rapid heating of material that may affect the nozzle dynamics or the dynamics of the accretion disk. As both the production of the jet and the seed photons emitted by the disk are dependent on the instantaneous value of $\dot{M}_{\rm h}$, the variability introduced by the fallback can have nonlinear feedback effects which may be able to produce the observed variations in flux, as is often seen in three-dimensional magnetohydrodynamical simulations of thick accretion disks (e.g., Hawley 2009)

To model this variability, we use the method of Kelly et al. (2011), in which the stochastic luminosity variation is modeled as a single Ornstein–Uhlenbeck process with characteristic frequency ω = 2π/τ_t, which is then added to the mean luminosity $L_{\rm j}\propto \dot{M}_{\rm h}$. To match the basic character of the observed light curve, we adjust the amplitude of the driving noise depending on the degree to which L_j exceeds the Eddington ratio, with a minimum noise amplitude chosen to match the observed variability at late times. At late times, the driving noise we have selected is within a factor of a few of what is measured for steady-state active galactic nuclei (AGNs).

Figure 14 shows the results of this simple model in comparison to the Swift event for 10³ random realizations of a flare generated by the disruption of a 1 M_☉ star by a 10⁶ M_☉ black hole, along with the realization that happens to give the best overall fit. The observed power, assumed to follow $\dot{M}_{\rm h}(t)$, is amplified by two powers of the Doppler factor, δ, with δ² = 20 and with 10% of the luminosity assumed to be emitted in the 2–10 keV energy range. For each realization, the triggering time is set when the isotropic equivalent luminosity of 10^47.5 erg s⁻¹ is surpassed. In general, this occurs close to the peak luminosity for the choice of parameters we have selected, with the trigger itself being caused by the apex of one of the mini-flares associated with our variable component. This leads to an envelope shape that shows a subtle decline from the initial triggering event within a few τ_t, and then flattens out over tens of days until t > a few times t_peak. If the light curve near t_peak is characterized by a large degree of variability, as our models predicts, the triggering event occurs either before or after t_peak with almost equal probability, and thus identifying the actual time of peak within the light curve itself would be difficult, if not impossible.

Zoom In Zoom Out Reset image size

As argued in Section 4.4, even for the simplest assumption of $L_{\rm j}\propto \dot{M}_{\rm h}$, we find that the resulting structure and dynamics of the radio hot spots are very different from those predicted by the standard spherical (one-dimensional) solutions (Bloom et al. 2011; Berger et al. 2012; Metzger et al. 2012). This is mainly because at early stages when t ⩽ t_peak the jet evolution is governed both by the advance of the head and by the cocoon's pressure-driven sideways expansion into the ambient medium. The rapid decrease in luminosity at t_peak would be preceded by the collapse of the jet channel. The collapse of both the channel and the cocoon would proceed from the galactic nuclei outward. A model in which the hot spots are simply turned off is incompatible with hydrodynamical collapse of the source. The fact that information about the rapid decrease in luminosity at t ≳ t_peak is communicated to the jet head at t_r ≫ t_peak (Figure 14) provides an attractive explanation (without the need of late time energy ejection as advocated by Berger et al. 2012) for the delay between the jet power X-ray luminosity and the radio afterglow (thought to be produced by the synchrotron-emitting electrons at the jet's head) seen in Sw 1644+57 (Berger et al. 2012). The most difficult task at present is to predict t_r, the time at which the intensity of the radiation from the radio-emitting electrons which have recently emerged from the hot spot starts to rapidly decrease, which corresponds to a time ≈r_r/(cβ_hΓ²_h) as measured by an observer along the line of sight. This is mainly because the hydrodynamical collapse of the channel depends sensitively on the initial Lorentz factor and opening angle of the jet as well as on the poorly known state of the surrounding environment, from the initial density gradient created by stellar wind injection to the large-scale ambient structures transversed by the jet's head at r ≳ r_r.

The rate of transient events coming as the result of the tidal disruption of a main-sequence star is highly dependent on the population of black holes and their immediate environments. The steady-state rate of stars entering the loss cone has been estimated by many authors (Bahcall & Wolf 1976; Lightman & Shapiro 1977; Shapiro 1977; Syer & Ulmer 1999; Wang & Merritt 2004), but these calculated rates ignore mass segregation, the mass function of stars within the sphere of influence (Bahcall & Wolf 1977), the importance of black hole mergers (Chen et al. 2009; Stone & Loeb 2011), and interactions between the central cluster and the host galaxy (Magorrian & Tremaine 1999; Merritt 2009), all of which may affect the rate of stellar disruption. As implied by the number of AGN flares thought to result from tidal disruptions (Gezari et al. 2009), the average rate of disruption is known to be ∼10⁻⁵ events yr⁻¹ gal⁻¹, but how these events are distributed among galaxies remains an open question. Given these uncertainties, we assume that the rate of disruption is not dependent on either black hole or galaxy properties for the remainder of this section.

Jets are thought to form in optically thick accretion flows resulting from accretion rates that are in excess of Eddington (Quataert 2001; Narayan & McClintock 2008). The mass function of stars around black holes may be substantially different than the standard Kroupa initial mass function (Kroupa 2001), but the average stellar mass is on the order of 1 M_☉ for all black holes. As a result, the peak accretion rate simply scales with the orbital period at the tidal radius, which means low-mass black holes can more easily produce jets from the disruption of a star given their smaller Eddington limits. This biases the production of jets to low-mass black holes, of which our census is growing (Ramya et al. 2011; Xiao et al. 2011) but is still too small to definitively relate to galactic velocity dispersion as has been demonstrated convincingly with black holes of M_h ≳ 10⁶ M_☉ (Merritt & Ferrarese 2001; Gültekin et al. 2009; Beifiori et al. 2012). Emission resulting from disks, on the other hand, is thought to extend well below the Eddington rate, and thus the only requirement to produce a transient disk from the disruption of a star is that the pericenter of the orbit does not lie within the black hole's Schwarzschild radius. As a result, disks are not likely to be produced by disruptions for black holes of mass ≳ 10⁸ M_☉.

Figure 15 shows the fraction of the total bolometric flux received by an observer from both the jet and disk components resulting from the tidal disruptions of a 1 M_☉ stars. The majority of stars that are fed into the black hole originate near the black hole's sphere of influence from the pinhole regime of disruption (Lightman & Shapiro 1977), and thus the differential contribution of stars as a function of the impact parameter β approximately scales as β⁻². The left panel in Figure 15 shows the expected contribution of events to each emission type assuming that black holes of all masses are equally common, an assumption which is not unreasonable given the potentially flat black hole mass function at low black hole mass, as measured by Greene & Ho (2007). If the trend of M_h − σ continued to lower black hole masses, this plot would be further biased to black holes of lower mass given the large population of low-mass galaxies. If the trend does not continue and black holes of ≲ 10⁶ M_☉ are rare, both emission components are biased to black holes of larger mass, as shown in the right panel of Figure 15. The limits on the mass of the black hole associated with the transient event Sw 1644+57 coming from the variability timescale and the M_h − σ relation (Bloom et al. 2011; Burrows et al. 2011; Miller & Gültekin 2011) agree quite well with a black hole mass function that excludes low-mass black holes, but given the associated uncertainties in the mass determination even a completely flat distribution of black hole masses cannot be completely ruled out.

Zoom In Zoom Out Reset image size

An additional constraint on the progenitors of events similar to Sw 1644+57 is the time at which $\dot{M}_{\rm h}$ reaches its peak value, t_peak, and the time at which $\dot{M}_{\rm h}$ drops below M_Edd, t_Edd. Under the same assumptions used to generate Figure 15, we show in Figure 16 the probability P of observing a flare with a given value of t_peak and t_Edd. We find that the most probable t_peak is a few to tens of days, whereas the most probable t_Edd ∼ 1–10 years, depending on whether or not low-mass central black holes are common. Both distributions are consistent with the Swift event, which our model suggests was caught within days of the peak, and remained above the Eddington limit for ∼1 year.

Zoom In Zoom Out Reset image size

The discovery of flaring black hole candidates in nearby galaxies by Swift will continue to elucidate the demography of the quasar and AGN population, while also enabling us to distinguish between various theoretical models of tidal disruption. As the feeding rate evolves dramatically over only a period of weeks to months, the powering of these formerly dead quasars offers a unique perspective into how the production of jets from black holes depends on the rate at which they are fed and the state of the environment in which they reside. These events can provide valuable insight into the physical mechanisms that operate near massive black holes that are not continually active, a prospect that is otherwise only possible through the study of the nearest few massive black holes. As the number of events increases, the range of possible models that can simultaneously explain the luminosity, variability, color, and time evolution of these events will diminish, enabling a firm characterization of one of the dominant feeding mechanisms of massive black holes in the local universe.