Estimating the Jet Power of Mrk 231 during the 2017–2018 Flare

The first 43 GHz observations were not triggered but were obtained to see the benefit of higher resolution. They occurred in two epochs in 2006, separated by three months. In hindsight, the observations were in a relatively low state for the core, ∼20 mJy at 43 GHz (compare to Table 2). The results were analyzed in Reynolds et al. (2009). Three observation frequencies were chosen: 15.3, 22.2, and 43.1 GHz. Using archival data, it was determined that the time variability indicated that the line of sight (LOS) was restricted in order to avoid the inverse Compton catastrophe, $\lt 25\buildrel{\circ}\over{.} {6}_{-2\buildrel{\circ}\over{.} \,6}^{+3\buildrel{\circ}\over{.} \,2}$. The core spectrum changed dramatically between the two epochs, with the 22 GHz flux density more than doubling. We were able to make a crude estimate of the jet power based on this spectral change: Q ≈ 3 × 10⁴³ erg s⁻¹. It was also argued based on the lower frequency VLBA in Ulvestad et al. (1999) that the spectrum of K1 is free–free absorbed near 5.4 GHz and 1.6 GHz. Based on the emission measure, K1 was interpreted as a radio lobe that results from the jet emitted from the core being stopped by the entrainment of the BAL wind.

It was determined that the 2006 campaign lacked the frequency coverage to properly describe the core spectrum, i.e., the SSA turnover. Thus, we proposed 8.4 GHz observations as well as 15.3, 22.2, and 43.1 GHz observations. Furthermore, the observations were triggered by a strong flare at 17.6 GHz. Our primary goal was to detect a discrete ejection. We did find one clear and one marginal superluminal ejections at 43 GHz, but they were very faint, ∼1–2 mJy, and faded very quickly (Reynolds et al. 2017). Thus, they were not detectable in subsequent epochs two to four weeks later. Our three epochs of NuSTAR observations did not detect any X-ray evolution during the flare. The VLBA coverage stopped during the flare peak. There was not enough information to constrain the jet power. The spectra appeared to be SSA steep spectrum power laws (see Figure 4 for examples of SSA power laws in the current campaign), except in the last epoch. Curiously, in the last epoch, the core spectrum changed to a flat spectrum power law with no evidence of SSA at frequencies lower than 8.4 GHz.

The clear detection of the ejection of a weak superluminal component was used in conjunction with Doppler aberration arguments to restrict the LOS to <235, similar to the upper bound found in the 2006 campaign from the time-variability Doppler argument (Reynolds et al. 2009).

Even though a large amount of information was lost due to poorer than expected resolution, we did anticipate the additional information that can be found in the core spectrum with a multifrequency observation plan. The spectra of the unresolved core in each epoch are shown in Figure 4. In Reynolds et al. (2009), we used the detailed nature of the spectral shape to argue that the spectrum of the radio core was most likely a power-law synchrotron spectrum that was seen through an SSA screen as opposed to a synchrotron power-law spectrum that was seen through an opaque screen caused by free–free absorption. The plots in Figure 4 show SSA power-law fits to the data. This is a simple spectral model that has been used to describe previous observations (Reynolds et al. 2009, 2017). In Section 6, we will interpret these fits in terms of the physical models that emit these spectra.

Zoom In Zoom Out Reset image size

Figure 5 shows the other important piece of information that can be retrieved from Table 2 about the core dynamics during the flare. It has been shown that K1 is stationary (with respect to the phase reference calibrator) within the uncertainties of the VLBA over two decades (Ulvestad et al. 1999, 1999; Reynolds et al. 2009, 2017). Thus, as in Reynolds et al. (2017), we use K1 as the origin for an astrometric description of core motion during the flare. There is more scatter in the 2017 data than we see in the 2015 data, due to the loss of observing stations in the current analysis. We do not consider the poor epoch C observation in this analysis, due to the problems discussed earlier in this section. The apparent velocity is v_app ≈ 0.97c based on the best-fit line, directed toward K1. What we identify as the core appears to be a discrete ejection that dominates the total nuclear flux density at 43 GHz. The nuclear region and the ejected components are too small to be resolved based on the fits in Table 2. The displacements are about one-third of the beam width (Figure 2) in our two best observation epochs A and B, and less in later epochs. The next section will present models of the physical nature of the discrete high-brightness feature moving toward K1.

Zoom In Zoom Out Reset image size

It is important to distinguish between quantities measured in the plasmoid frame of reference and those measured in the Earth observer frame of reference. The strategy will be to evaluate the physics in the plasma rest frame then transform the results to the observer's frame for comparison with observation. The underlying power law for the flux density is defined as ${S}_{\nu }(\nu ={\nu }_{o})=S{\nu }_{o}^{-\alpha }$, where S is a constant. Observed quantities will be designated with a subscript, "o," in the following expressions. The observed frequency is related to the emitted frequency, ν, by ν_o = δν. The bulk flow Doppler factor of the plasmoid

$$\begin{eqnarray}&&\delta =\displaystyle \frac{{\gamma }^{-1}}{1-\beta \cos \theta },\ {\gamma }^{-1}=1-{\beta }^{2},\end{eqnarray} \tag{ 1 }$$

where β is the normalized three-velocity of bulk motion and θ is the angle of the motion to the LOS to the observer. The SSA attenuation coefficient in the plasma rest frame is given by (Ginzburg & Syrovatskii 1969)

$$\begin{eqnarray}\begin{array}{rcl}\mu (\nu ) & = & \overline{g(n)}\displaystyle \frac{{e}^{3}}{2\pi {m}_{e}}{N}_{{\rm{\Gamma }}}{\left({m}_{e}{c}^{2}\right)}^{2\alpha }\\ & & \times {\left(\displaystyle \frac{3e}{2\pi {m}_{e}^{3}{c}^{5}}\right)}^{\tfrac{1+2\alpha }{2}}{\left(B\right)}^{(1.5+\alpha )}{\left(\nu \right)}^{-(2.5+\alpha )},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\overline{g(n)} & = & \displaystyle \frac{\sqrt{3\pi }}{8}\\ & & \times \displaystyle \frac{\overline{{\rm{\Gamma }}}[(3n+22)/12]\overline{{\rm{\Gamma }}}[(3n+2)/12]\overline{{\rm{\Gamma }}}[(n+6)/4]}{\overline{{\rm{\Gamma }}}[(n+8)/4]},\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&N={\int }_{{{\rm{\Gamma }}}_{\min }}^{{{\rm{\Gamma }}}_{\max }}{N}_{{\rm{\Gamma }}}{{\rm{\Gamma }}}^{-n}\,d{\rm{\Gamma }},\ n=2\alpha +1,\end{eqnarray} \tag{ 4 }$$

where Γ is the ratio of lepton energy to rest-mass energy m_ec², and $\overline{{\rm{\Gamma }}}$ is the gamma function. B is the magnitude of the total magnetic field. The power-law spectral index for the flux density is α = (n − 1)/2. The low-energy cutoff, E_min = Γ_minm_ec², is constrained loosely by the data in Table 3 and Figure 4. The fact that the core is very luminous at ν_o = 8.4 GHz means that the lepton energy distribution is not cut off near this frequency. This is not a very stringent bound as it is just a consequence of the fact that 8.4 GHz is our lowest observing frequency. It would be a major coincidence if this were not to provide a loose upper bound on Γ_min. Note that the SSA opacity in the observer's frame, μ(ν_o), is obtained by direct substitution of $\nu ={\nu }_{o}/\delta$ into Equation (2).

Table 3.  Fits to the Core Spectra

Epoch	Date	Peak Frequency	Peak Luminosity	Spectral Index	Number Index
	(MJD)	(GHz)	(erg s−1)	α	n
		νpeak	$\nu {L}_{\nu }(\nu ={\nu }_{\mathrm{peak}})$
A	58063	13.9	8.08 × 1040	0.83	2.66
B	58077	14.1	1.12 × 1041	0.75	2.50
C	58102	13.2	1.05 × 1041	0.97	2.94
D	58132	11.0	8.15 × 1040	1.03	3.06

Download table as:  ASCIITypeset image

The uniform, homogeneous approximation yields a simplified solution to the radiative transfer equation (Ginzburg & Syrovatskii 1965; van der Laan 1966)

$$\begin{eqnarray}\begin{array}{rcl}{S}_{{\nu }_{\unicode{x000F8}}} & = & \displaystyle \frac{{S}_{o}{\nu }_{o}^{-\alpha }}{\tau ({\nu }_{o})}\times \left(1-{e}^{-\tau ({\nu }_{o})}\right),\\ \tau ({\nu }_{o}) & \equiv & \mu ({\nu }_{o})L,\ \tau ({\nu }_{o})=\overline{\tau }{\nu }_{o}^{(-2.5+\alpha )},\end{array}\end{eqnarray} \tag{ 5 }$$

where τ(ν) is the SSA opacity, L is the path length in the rest frame of the plasma, S_o is a normalization factor, and $\overline{\tau }$ is a constant. There are three unknowns in Equation (5): $\overline{\tau }$, α, and S_o. These are the three constraints on the following theoretical model that are estimated from the observational data. These three constraints are used to establish the uniqueness and existence of solutions in the following.

The theoretical spectrum is parameterized by Equations (2)–(5) and the synchrotron emissivity that is given in Tucker (1975) as

$$\begin{eqnarray}&&{j}_{\nu }=1.7\times {10}^{-21}[4\pi {N}_{{\rm{\Gamma }}}]a(n){B}^{(1+\alpha )}{\left(\displaystyle \frac{4\times {10}^{6}}{\nu }\right)}^{\alpha },\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&a(n)=\displaystyle \frac{\left({2}^{\tfrac{n-1}{2}}\sqrt{3}\right)\overline{{\rm{\Gamma }}}\left(\tfrac{3n-1}{12}\right)\overline{{\rm{\Gamma }}}\left(\tfrac{3n+19}{12}\right)\overline{{\rm{\Gamma }}}\left(\tfrac{n+5}{4}\right)}{8\sqrt{\pi }(n+1)\overline{{\rm{\Gamma }}}\left(\tfrac{n+7}{4}\right)}.\end{eqnarray} \tag{ 7 }$$

One can transform this to the observed flux density, S(ν_o), in the optically thin region of the spectrum using the relativistic transformation relations from Lind & Blandford (1985),

$$\begin{eqnarray}&&S({\nu }_{o})=\displaystyle \frac{{\delta }^{(3+\alpha )}}{4\pi {D}_{L}^{2}}\int {j}_{\nu }^{{\prime} }{dV}^{\prime} ,\end{eqnarray} \tag{ 8 }$$

where D_L is the luminosity distance, and in this expression, the primed frame is the rest frame of the plasma. These are the basic equations needed to fit the data in Figure 4.

The apparent velocity constrains the kinematics of the plasmoid (Rees 1966):

$$\begin{eqnarray}&&{v}_{\mathrm{app}}/c={\beta }_{\mathrm{app}}=\beta \sin \theta /(1-\beta \cos \theta )\approx 0.97,\end{eqnarray} \tag{ 9 }$$

where we used the result from Figure 5. The apparent motion and the LOS upper limit $25\buildrel{\circ}\over{.} {6}_{-2\buildrel{\circ}\over{.} \,6}^{+3\buildrel{\circ}\over{.} \,2}$ noted in Section 2.1 constrain the dynamics. The validity of these numbers is a major assumption of this analysis.

The models that produce the fits in Figure 4 are of interest if we can deduce the physical parameters in Equations (1)–(9) that are responsible for the spectra. Furthermore, we want to know the kinematic consequences of these parameters within a specific model. There are two basic classes of models. The first are discrete ejections of plasma, a ballistic ejection of a plasmoid. This is essentially a blob of magnetized gas. The gas is turbulent, and the magnetic field is randomized and tangled (Moffet 1975). In this model, the magnetic field strength behaves as an extra component of the internal gas pressure and acts as an expansive force on the plasmoid (Punsly 2008). Second, the discrete unresolved emission might be "knots" or regions of enhanced dissipation within a continuous jet. In the jet scenario, the magnetic field is predominantly an organized toroidal magnetic field, ${B}_{\phi }^{{\prime} }$, in the rest frame of the plasma in contrast to the turbulent magnetic field of the discrete ejections (Blandford & Königl 1979). The behavior of this magnetic field is much different than that of the turbulent magnetic field. It acts as a hoop stress that provides a confining force on the outflow (Punsly 2008).

5.2.1. The Poynting Flux in a Strong Knot in the Jet

We review the derivation of the dominance of the toroidal component of the magnetic field in order to derive an estimate of the poloidal Poynting flux. The dominance of the toroidal component of the magnetic field is a consequence of the perfect magnetohydrodynamic (MHD) assumption and approximate angular momentum conservation in the jet (Blandford & Königl 1979). The total angular momentum flux per unit magnetic flux in the observer's coordinate system is (Punsly 2008)

$$\begin{eqnarray}&&L\approx k\mu {r}_{\perp }\gamma {v}^{\phi }-\displaystyle \frac{c}{4\pi }{B}^{\phi }{r}_{\perp },\ k\equiv \displaystyle \frac{N\gamma {v}^{P}}{{B}^{P}},\end{eqnarray} \tag{ 10 }$$

where k is the perfect MHD conserved mass flux per unit poloidal magnetic flux, and B^P and v^P are the poloidal magnetic field and the poloidal velocity, respectively. The quantity N is the number density in the plasma rest frame, μ is the specific enthalpy, and r_⊥ is the cylindrical radius. The first term is the mechanical angular momentum flux per unit magnetic flux, and the second term is the electromagnetic angular momentum flux per unit magnetic flux. The notion of describing the angular momentum per unit poloidal magnetic flux is physically motivated by angular momentum conservation and perfect MHD. In perfect MHD, angular momentum is transported conserved along each poloidal magnetic field line in a wind or jet (Punsly 2008). Jets expand as they propagate (Blandford & Königl 1979). Thus, r_⊥ is much larger in the knot than at the jet base. Consider a jet that is highly magnetic, either Poynting flux–dominated or at a minimum the mechanical angular momentum flux per unit magnetic flux and the electromagnetic angular momentum flux per unit magnetic flux are comparable. Then, by Equation (10), if L is approximately constant, in the limit of large r_⊥, the following must be true for highly magnetic jets (assuming angular momentum outflow, i.e., L > 0, B^P > 0, B^ϕ < 0):

$$\begin{eqnarray}&&{v}^{\phi }\lt \displaystyle \frac{L}{k\mu {r}_{\perp }\gamma },\ 0.5\displaystyle \frac{4\pi }{{{cr}}_{\perp }}L\lt | {B}^{\phi }| \lt \displaystyle \frac{4\pi }{{{cr}}_{\perp }}L.\end{eqnarray} \tag{ 11 }$$

The lower bound in the second relationship merely expresses the fact that by assumption (magnetic jet), the electromagnetic angular momentum flux is at least equal to the mechanical angular momentum flux. By poloidal magnetic flux conservation B^P ∼ r_⊥⁻², so by Equation (11), $| {B}^{\phi }| \gg {B}^{P}$ at large r_⊥. To estimate the poloidal Poynting flux, S^P, in the plasmoid, first transform fields to the observer's frame:

$$\begin{eqnarray}&&{B}^{\phi }=\gamma {B}_{\phi }^{{\prime} }\quad {E}^{\perp }=\displaystyle \frac{{v}^{P}}{c}\gamma {B}_{\phi }^{{\prime} }-\displaystyle \frac{{v}^{\phi }}{c}\gamma {B}^{P}\approx \displaystyle \frac{{v}^{P}}{c}\gamma {B}_{\phi }^{{\prime} },\end{eqnarray} \tag{ 12 }$$

where ${E}^{\perp }$ is the poloidal electric field orthogonal to the magnetic field direction. The poloidal Poynting flux in the observer's frame, S^P, along the jet direction is (Punsly 2008)

$$\begin{eqnarray}&&{S}^{P}=\displaystyle \frac{c}{4\pi }{E}^{\perp }{B}^{\phi }\approx \displaystyle \frac{c}{4\pi }{\gamma }^{2}\beta {\left[{B}_{\phi }^{{\prime} }\right]}^{2},\ {B}_{\phi }^{{\prime} }\approx B.\end{eqnarray} \tag{ 13 }$$

5.2.2. Mechanical Contributions to the Energy Flux

The energy content is separated into two pieces. The first is the kinetic energy of the protons, ${ \mathcal K }(\mathrm{protonic})$,

$$\begin{eqnarray}&&{ \mathcal K }(\mathrm{protonic})=(\gamma -1){{Mc}}^{2},\end{eqnarray} \tag{ 14 }$$

where M is the mass of the plasmoid. The other piece is named the lepto-magnetic energy, $E(\mathrm{lm})$, and is composed of the volume integral of the leptonic internal energy density, U_e, and the magnetic field energy density, U_B. It is straightforward to compute the lepto-magnetic energy in a spherical volume,

$$\begin{eqnarray}\begin{array}{rcl}E(\mathrm{lm}) & = & \displaystyle \int ({U}_{B}+{U}_{e})\,{dV}=\displaystyle \frac{4}{3}\pi {R}^{3}\\ & & \times \left[\displaystyle \frac{{B}^{2}}{8\pi }+{\displaystyle \int }_{{{\rm{\Gamma }}}_{\min }}^{{{\rm{\Gamma }}}_{\max }}({m}_{e}{c}^{2})({N}_{{\rm{\Gamma }}}{E}^{-n+1})d\,E\right].\end{array}\end{eqnarray} \tag{ 15 }$$

The lepto-magnetic energy is often argued to be minimized in astrophysical sources (Hardcastle et al. 2004; Croston et al. 2005; Kataoka & Stawarz 2005). The relevance of this will be discussed in terms of the time evolution of various models of the flare in the next section. The leptons also have a kinetic energy analogous to Equation (14),

$$\begin{eqnarray}&&{ \mathcal K }(\mathrm{leptonic})=(\gamma -1){{ \mathcal N }}_{e}{m}_{e}{c}^{2},\end{eqnarray} \tag{ 16 }$$

where ${{ \mathcal N }}_{e}$ is the total number of leptons in the plasmoid.

The other quantities of interest are the protonic and leptonic field-aligned (jet-axis-aligned) poloidal energy fluxes in the frame of the observer. The protonic energy flux in the frame of the observer is approximately the kinetic energy flux density given by

$$\begin{eqnarray}&&{ \mathcal E }(\mathrm{proton})=N(\gamma -1)\gamma {v}^{P}{m}_{p}{c}^{2},\end{eqnarray} \tag{ 17 }$$

where m_p is the mass of the proton. The leptonic jet aligned poloidal energy flux density in the frame of the observer is

$$\begin{eqnarray}&&{ke}(\mathrm{leptonic})=N\gamma {v}^{P}\left[\gamma \mu {c}^{2}\right],\end{eqnarray} \tag{ 18 }$$

where k is the mass flux defined in Equation (10) and e is the specific energy of a lepton (Punsly 2008). From an energetics standpoint, it is useful to subtract off the rest-mass term, which merely reflects the conservation of particle number from the source to the jet and cannot be converted into different forms of energy flux. This has been called the free thermo-kinetic energy flux and was defined in McKinney et al. (2012) as

$$\begin{eqnarray}&&{ \mathcal E }(\mathrm{leptonic})=N\gamma {v}^{P}\left[\gamma \mu {c}^{2}-{m}_{e}{c}^{2}\right].\end{eqnarray} \tag{ 19 }$$

This would be appropriate for protonic outflows. However, for an electron–positron outflow, the mass is not conserved due to pair creation, and in fact in the models considered, the number of leptons in the plasmoid increases substantially over time. Thus, ${ke}(\mathrm{leptonic})$ is the more appropriate energy flux in this study. The specific enthalpy decomposes as

$$\begin{eqnarray}&&N\mu ={U}_{e}+P,\end{eqnarray} \tag{ 20 }$$

where the relativistic pressure $P\approx (1/3)({U}_{e}-{{Nm}}_{e}{c}^{2})$ (Willott et al. 1999).

Figure 7 shows the dependence of E(lm) in Equation (15) on R for a wide span of preassigned values of δ and E_min. The figure allows one to see how E(lm)(R) varies as one adjusts the preassigned δ and E_min. More importantly, it shows how the plasmoid changes, independent of the preassigned δ and E_min from epoch A to epoch D. The dashed line and the red dots in each panel show the time evolution of a possible physical solution. It is the solution found for the major flares in GRS 1915+105 (Punsly 2012). A blob of turbulent electron–positron plasma is ejected from the black hole accretion system. It obeys energy conservation, and the radiation losses are negligible in the 69 days of monitoring (∼10⁴⁷ erg). To the right of the minimum of each E(lm) curve, the stored plasmoid energy is more magnetic than mechanical and to the left of the minimum of E(lm) of each curve in Figure 7, the energy is more mechanical than magnetic. The dashed curve indicates a scenario in which the flare evolves from magnetically dominated to near equipartition at late times. This satisfies physical requirement (3) from Section 6. Magnetic energy is converted to mechanical form. This is consistent with physical requirement (2) of Section 6, magnetic forces can launch the jet, and this energy is converted to bulk motion and internal leptonic energy.

Zoom In Zoom Out Reset image size

The main concern is physical requirement (1) from Section 6. We can estimate the radio luminosity of K1 from the VLBA observations. For example, in our best observation, epoch A, we note that the K1 luminosity is down about 10% from its historic high in 2015. The data in Table 2 is well fit by a power law with α = 1.61. Most of the radio luminosity is at low frequency, but our data only go down to 8.4 GHz. A previous VLBA campaign in 1996 covered frequencies of 1.5, 2.3, 4.9, 8.4, and 15.4 GHz (Ulvestad et al. 1999). Even though the double is not resolved (even partially) below 8.4 GHz, the core was very weak at 8.4 and 17.6 GHz (∼15% of that in 2017). Thus, the low-frequency emission can all be considered to be from K1 as a good approximation. It was argued that the power-law behavior began to be attenuated by free–free absorption below 8.4 GHz (Ulvestad et al. 1999). The flux density peaked at 5 GHz. Thus, we estimate the radio luminosity from K1 in epoch A by integrating the power-law fit from 5 to 50 GHz, ≈8.4 × 10⁴⁰ erg s⁻¹. Thus, from 2015 to 2018, this equates to ∼2.5 × 10⁴⁸ erg yr⁻¹. The total energy supplied every year must be 2.5 × 10⁴⁸/_R erg, where _R is the radiative efficiency of the plasma from 5 GHz to 50 GHz in K1. For the case LOS = 145 and Γ_min = 1, in Figure 7, there needs to be $\approx [2.5\,\times \,{10}^{48}\,\mathrm{erg}\,{\mathrm{yr}}^{-1}/1.25\,\times \,{10}^{49}\,\mathrm{erg}]/{\epsilon }_{R}\approx 0.2/{\epsilon }_{R}$ ejections per year, as the light curve in Figure 1 shows at most two ejections per year. This implies a radiative efficiency of larger than 10% in K1 which would be a very extreme requirement for the radio lobe, but such a high efficiency might be consistent with a compact steep spectrum object (O'Dea 1998).

Another possibility is a ballistic ejection of a blob of protonic plasma. Figure 8 plots ${ \mathcal K }(\mathrm{protonic})$ as a function of E(lm) for the case where LOS = 145 and Γ_min = 1. This is useful for understanding the time evolution of the system. For a flare that decays, we posit that the final state is near minimum energy, minimum E(lm), as in the discrete ejections in GRS 1915+105. If this is true, then a solution associated with the dashed black line is indicated. The plasmoid is ejected with ${ \mathcal K }(\mathrm{protonic})=2.7\times {10}^{51}$ erg, a constant total energy from epoch A to epoch D. Epoch D is relaxed to the minimum energy configuration. Figure 8 shows that the time evolution is one in which plasma internal energy is converted to magnetic turbulent energy. Based on Figure 1, the longest time frame in which energy is pumped into the plasmoid is the time lapse, T, from epoch A to the start of the flare (∼120 days; see Section 8 for more details). The instantaneous jet power required to energize and eject the plasmoid is therefore ${Q}_{o}\approx { \mathcal K }(\mathrm{protonic})/T\gt 2.7\times {10}^{44}$ erg s⁻¹. This is a large jet power by astrophysical standards. Theoretical efforts to explain the initiation of such powerful jets rely on some form of magnetic jet launching that requires the plasma dynamics to be dominated by magnetic forces (Blandford 1976; Lovelace 1976; Blandford & Payne 1982). Thus, it is hard to understand a mechanism in which inertial energy is converted to magnetic energy during the course of the plasmoid evolution. One would expect there to be an excess of magnetic energy in the early stages as required for relativistic jet initiation. The indicated time evolution does not satisfy physical requirement (2) of Section 6.

Zoom In Zoom Out Reset image size

Another possibility that is very common in extragalactic jet models is a magnetic jet filled with an electron–positron plasma (Blandford & Königl 1979; Willott et al. 1999). In this scenario, the spherical plasmoid is a knot in the jet as indicated in Figure 9. The knot represents a region of enhanced dissipation. It might be related to a wavefront or shock signaling an increase in jet power from the source that is propagated down the conduit of the jet. In this circumstance, the following estimates of jet power represent a transient state of enhanced jet power. Alternatively, the enhanced dissipation might have arisen from interactions of the jet with the ambient medium. For example, it is likely that the jet is interacting with the very dense BAL wind (Reynolds et al. 2009). This provides a piston for launching strong dissipative waves such as shocks as well as a source of instabilities created from the jet/wind interface (Bicknell et al. 1990). In this instance, the jet power estimate is likely more indicative of the jet over a long time frame, but the interaction has allowed us to sample the quasi-steady jet (by increasing the dissipation) for a short period of time.

Zoom In Zoom Out Reset image size

The jet power in this model, Q, using Equations (13) and (18) is

$$\begin{eqnarray}&&Q=\int [{S}^{P}+{ke}(\mathrm{leptonic})]{{dA}}_{\perp }+{L}_{r},\end{eqnarray} \tag{ 22 }$$

where dA_⊥ is the cross-sectional area element normal to the jet axis and L_r is the energy flux lost to radiation. We constrain the solution space by assuming that the system approaches minimum E(lm) at late times (epoch D). This is motivated by the discrete ejections in GRS 1915+105 that approach minimum energy at late times, and it has been argued that the hot spots (the terminal jet knot) and the surrounding radio lobes are near minimum energy in extragalactic radio sources (Hardcastle et al. 2004; Croston et al. 2005; Kataoka & Stawarz 2005; Punsly 2012). A common spheroid shape for the emission region is chosen to allow direct comparison between the four classes of models. However, it is not ideal for cross-sectional integration as in Equation (22), due to its nonuniform cross-sectional area. Thus, the evaluation of Q in Equation (22) is the average over all the cross sections (i.e., $\pi {R}^{2}\to (2/3)\pi {R}^{2}$). Within the crude approximation of these simplified models, this is not a huge source of uncertainty. The basic physics of the knot evolution is shown in Figure 10. The top frame plots Q versus E(lm). Epoch D is in the minimum energy configuration. The time evolution between epochs is indicated by the dashed line, which is simply conservation of Q (radiation losses are negligible). The evolution makes sense from jet/wind launching theory, requirement (2) of Section 6. Magnetic energy flux is converted to inertial energy flux. The bottom frame of Figure 10 is a plot of the fraction of the jet power that is inertial energy flux:

$$\begin{eqnarray}&&\mathrm{Inertial}\,\ \mathrm{Fraction}=\left[\int {ke}(\mathrm{leptonic})){{dA}}_{\perp }\right]/Q.\end{eqnarray} \tag{ 23 }$$

The dashed line is the same solution as in the top frame. In epoch A, the inertial component of Q is ∼5% of Q and increases steadily to epoch D, where it is ∼40%.

Zoom In Zoom Out Reset image size

The dependence of jet power on the LOS is illustrated in the top-left frame of Figure 11. A larger LOS equates to a more powerful jet as well as a stronger discrete ejection in the leptonic model of Section 7.1. The top-right frame of Figure 11 shows that Q decreases as E_min increases. It also shows that a protonic jet (see next section) is much more powerful. The lower left-hand frame shows the time evolution of the spheroid radius. The solutions with Γ_min = 1 were used to estimate the jet-opening angle at K1 in Figure 9. The solution with Γ_min = 1 and an LOS = 145 has an expansion rate 2(dR/dt) = 0.2c. The bottom right-hand corner of Figure 11 shows that the magnetic field decreases as the knot propagates. This is expected from the conversion of magnetic to inertial energy indicated in Figure 10. Combining these results with Figure 6 and Equation (21) motivated our choice of Γ_min = 30 as an upper limit for a viable lower energy cutoff. The solution has a smooth monotonic behavior after the peak of the flare at epoch B as might be expected for a plasmoid that relaxes to the minimum energy state.

Zoom In Zoom Out Reset image size

It is possible that the jet is made of primarily protonic not positronic plasma. Then, there are two inertial contributions to Q, one from the energized leptons and another from a dominant component from the bulk motion of protons. Equation (22) needs to be modified accordingly with Equation (17),

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{\mathrm{protonic}} & = & \displaystyle \int [{S}^{P}+{ke}(\mathrm{leptonic})+{ \mathcal E }(\mathrm{proton})]{{dA}}_{\perp }\\ & & +\,{L}_{r}\approx \displaystyle \int [{S}^{P}+{ \mathcal E }(\mathrm{proton})]{{dA}}_{\perp }.\end{array}\end{eqnarray} \tag{ 24 }$$

The solutions that are plotted in the upper right-hand frame of Figure 11 are computed assuming a minimum E(lm) state in epoch D, the same state as in all the other models. There is no concern with these solutions having enough power to energize K1 in accord with requirement (1) of Section 6. However, we are unable to find any solutions that satisfy physical requirements (2) and (3) that change monotonically and smoothly as the leptonic jet solutions of the previous subsection. Thus, our methods are not likely valid quantitatively, but might provide some qualitative insight into these solutions and will be discussed below. The biggest issue is that the assumption of a constant velocity is likely grossly inaccurate.

The top frame of Figure 12 shows the solution for Γ_min = 10 and an LOS = 145 that proceeds by converting magnetic energy to inertial energy. Based on the upper right-hand frame in Figure 11, any value of Γ_min < 10 is likely too powerful (A Fanaroff–Riley II, FR II, extragalactic radio source level luminosity). In order to explore protonic jets, it is useful to introduce the protonic energy density,

$$\begin{eqnarray}&&{U}_{P}\approx {{Nm}}_{p}{c}^{2}.\end{eqnarray} \tag{ 25 }$$

The horizontal axis in the top frame of Figure 12 plots U_P/U_B, the ratio of inertial energy density to magnetic energy density. Even though there are numerous changes from our original leptonic jet models inserted by the introduction of protons, the leptons in the plasma still radiate the spectra in Figure 4 for all the solutions within Figures 11 and 12.

Zoom In Zoom Out Reset image size

There are solutions that proceed toward more magnetic energy over time, but by physical requirement (2) of Section 6, this is untenable. The dashed arrow shows the direction of time evolution, assuming a constant jet power. The end of the arrowhead is the minimum energy solution in epoch D. This fixes the solution. Without this assumption, the solution is unconstrained, and our analysis is not useful. The only solutions that satisfy physical requirements (1)–(3) of Section 6 require a dramatic change between epochs C and D. This is apparent from the large gap in U_P/U_B from epoch C to epoch D in the top frame of Figure 12, two orders of magnitude. The bottom frame of Figure 12 indicates the required dynamics. From epochs A–C, there is a smooth gradual expansion with modest entrainment of plasma. The jet power is so dominated by Poynting flux (97%) that the jet can be either positronic or protonic in the early stages without affecting the solution. However, the bottom frame of Figure 12 indicates intense protonic entrainment between epochs C and D. This is not unreasonable as deduced in Reynolds et al. (2009). At K1, there is evidence of overwhelming entrainment on parsec scales from the dense enveloping BAL wind. The solution described here is one of a powerful Poynting flux–dominated jet in which all the power is converted into proton kinetic energy on the order of 1–2 light months from the source.

There are a couple of issues with this solution that violate our basic assumptions:

• 1.  
  Clearly, large entrainment would slow down the jet. Our models are not sophisticated enough to capture this, so the models are only qualitative.
• 2.  
  In order to keep Q reasonable, we had to choose an ad hoc Γ_min ≥ 10 based on the upper right-hand frame of Figure 11. The bottom frame of Figure 12 shows an entrainment solution in which over 90% of the plasma enters the knot between epochs C and D. The plasma that surrounds the jet in the BAL wind is cold relative to the jet. This process should drastically lower Γ_min. This complexity is far beyond our simple model and involves much uncertain physics and geometry.

These shortcomings raise the question of whether our model is of any value in the analysis of this scenario. In order to explore this, we note that the model of a Poynting flux–dominated jet with a steadily increasing radius and a slow and steady conversion of magnetic energy to mechanical energy seems qualitatively reasonable from a wind theory point of view (Blandford & Payne 1982; Punsly 2008). But, can we say anything about the notion of intense entrainment between epochs C and D that is implied by the prediction of the model? The model predicts its own failure in this time frame. First, we note that the knot in the Poynting jet is a very powerful dynamic object from epochs A–C. It propagates with β_app = 0.97c, Q > 4 × 10⁴³ erg s⁻¹ of pure Poynting flux, and it is expanding (based on the top frame of Figure 12) at 2(dR/dt) = 0.16c. In a one-month period, >95% of the Poynting jet power is converted to protonic kinetic energy flux. This process also introduces 90% of the jet plasma. The Poynting flux fraction of Q changes from 97% to 3%. The radius stops expanding at 0.08c and contracts rapidly. This seems to be a very violent process involving dissipative instabilities and fast shocks (Kennel & Coroniti 1984). The introduction of new plasma and strong dissipative process should change the energy spectrum of the electrons. This should be independent of the entrainment details. Now, we look outside of the simple, spherical (likely invalid due to the strong entrainment-induced changes) model for empirical data, in particular Table 3 and Figure 4. Note that α and n in the last two columns are very similar in epochs C and D, almost the same considering the crudeness of the model. The entrainment scenario should have produced a significant change unless there is a major coincidence that the process recreates the previous particle spectrum. By contrast, note the shift of the spectral peak and the decrease of peak flux intensity implied by Table 3 and Figure 4. This is the defining characteristic of cooling by adiabatic expansion: the same spectral shape with a shift of the spectral peak to lower frequency and lower intensity (Moffet 1975). This empirical (model independent) spectral change is consistent with the steady expansion of the leptonic jet model of the last subsection and not the violent contraction indicated for an the intense entrainment scenario of the protonic jet. For this reason, we consider a knot in a leptonic jet to be more likely to represent the physical solution than a knot in a protonic jet. We also note that the leptonic jet favors Γ_min ≈ 1 and does not require the ad hoc condition Γ_min ≥ 10 needed for the protonic jet models.

It is not trivial to estimate the accretion disk luminosity. Mrk 231 has significant intrinsic extinction in the nuclear region (Lipari et al. 1994; Smith et al. 1995; Veilleux et al. 2013, 2016). This not only makes the direct observation of the accretion disk spectrum impossible, but it also makes the origin of the weak X-ray emission from 2012 October in Figure 13 very uncertain (Veilleux et al. 2013; Reynolds et al. 2017). An estimator of the accretion disk luminosity is the Hα broad emission line, and this is very prominent in Mrk 231. The intrinsic extinction has been corrected with an LMC extinction law with E(B – V) = 0.63, yielding a plausible intrinsic quasar spectrum in the optical band (Lipari et al. 1994; Bokensberg et al. 1977). This same extinction law does not work in the ultraviolet (Smith et al. 1995; Veilleux et al. 2016). This extinction correction for Hα is a factor of ∼4.5. There is so much intrinsic absorption (including that of Fe ii) that Hβ is very weak and is difficult to use for estimating the bolometric luminosity of the accretion flow, L_bol. The situation is even worse for the ultraviolet broad lines (Smith et al. 1995). The Hα line is a bona fide quasar broad line with an extinction-corrected luminosity of L_Hα = 4.7 × 10⁴³ erg s⁻¹ and an FWHM of 2800 km s⁻¹ (Smith et al. 1995). L_bol can be estimated with the formula (Greene & Ho 2007)

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{bol}} & \approx & 2.34\times {10}^{44}{\left[\displaystyle \frac{{L}_{{\rm{H}}\alpha }}{{10}^{42}\mathrm{erg}{{\rm{s}}}^{-1}}\right]}^{0.86}\,\mathrm{erg}{{\rm{s}}}^{-1}\\ & \approx & 6.45\times {10}^{45}\ \mathrm{erg}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ 26 }$$

Second, we can also estimate L_bol from the reprocessed photons in the infrared from the dusty torus, but this is not simple in Mrk 231 either. The infrared band is very luminous and is actually dominated by starburst emission. The breakdown of the different components with detailed modeling indicates that the active galactic nucleus (AGN) luminosity at wavelengths longer than 3 μm is L_IR(AGN) ≈ 4 × 10⁴⁵ erg s⁻¹, and the starburst luminosity is L_IR(SB) ≈ 10⁴⁶ erg s⁻¹ (Farrah et al. 2003). The AGN component is complicated by the existence of two strong components, a hotter component with the SED peaked at ≈4 μm and a stronger cooler component peaked at ≈60 μm. The second component is likely associated with the molecular disk imaged in the hydroxyl maser line with VLBI with an inner radius of 30 pc and an outer radius of 100 pc (Klöckner et al. 2003). Regardless of this complex nature, we use the IR-based bolometric estimators (Spinoglio et al. 1995):

$$\begin{eqnarray}&&\mathrm{log}[{L}_{\mathrm{bol}}]=0.942\mathrm{log}[\lambda {L}_{\lambda }(\lambda =12\,\mu {\rm{m}})]+3.642\approx 45.82,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\mathrm{log}[{L}_{\mathrm{bol}}]=0.837\mathrm{log}[\lambda {L}_{\lambda }(\lambda =25\,\mu {\rm{m}})]+8.263\approx 45.75,\end{eqnarray} \tag{ 28 }$$

where $\lambda {L}_{\lambda }(\lambda =12\,\mu {\rm{m}})$ and $\lambda {L}_{\lambda }(\lambda =25\,\mu {\rm{m}})$ are from the AGN component of the fit to the IR SED (Farrah et al. 2003). The three estimation methods in Equations (26)–(28) closely agree with each giving us confidence in our estimate of L_bol.

The EIC luminosity, L_EIC, and frequency, ν_EIC, are related to the synchrotron luminosity, L_synch, and frequency, ν_synch, by (Tucker 1975)

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{EIC}}}{{L}_{\mathrm{synch}}}\approx \displaystyle \frac{{U}_{\mathrm{ph}}}{{U}_{B}},\quad \displaystyle \frac{{\nu }_{\mathrm{EIC}}}{{\nu }_{\mathrm{synch}}}\sim {{\rm{\Gamma }}}^{2},\end{eqnarray} \tag{ 29 }$$

where U_ph is the energy density of the accretion disk photon field in the frame of reference of the plasmoid. The flux of the photon field of the accretion disk is diluted in the plasmoid in epoch A by geometric dilution and also experiences Doppler redshifting. For the dissipative knot in a leptonic jet model (our preferred model in Section 7.3 and Figure 10), assuming an LOS = 145 and Γ_min = 1 for illustrative purposes, the relevant parameters are

$$\begin{eqnarray}&&B=0.41{\rm{G}}\quad {U}_{B}=6.69\times {10}^{-3}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\quad \beta =0.816.\end{eqnarray} \tag{ 30 }$$

In order to compute U_ph, taking into account geometric dilation, we need to know the distance the plasmoid was from the accretion disk during the NuSTAR observation. Based on the light curve in Figure 1, the plasmoid seems to have been ejected ∼100 days earlier. We illustrate this in Figure 14, which plots both the 17.6 GHz flux density and the spectral index from 14.6 to 17.6 GHz at the beginning of the flare. The spectral index is important because as a plasmoid is ejected at the beginning of the flare, it begins very compact and has large SSA at 17.6 GHz. Therefore, the flare begins at high frequency first and evolves to lower frequency (van der Laan 1966; Moffet 1975). For strong flares, even in the presence of the ∼120 mJy steep background flux density, the spectrum flattens as the new ejection emerges and can even invert (as in Figure 14) between 14.6 and 17.6 GHz (Reynolds et al. 2017). These conditions occur near MJD 57940; thus, this is a more likely flare ejection time than the other relative minimum near MJD 57974. We estimate ∼100 days interval between the flare ejection time and the NuSTAR observation.

Zoom In Zoom Out Reset image size

Assuming the velocity of Equation (30) is constant, it would have propagated D_d−p ∼ 2.12 × 10¹⁷ cm from the disk. In order to compute the relevant Doppler redshift, note that as the plasmoid moves away from the accretion disk, the distance from the disk becomes much larger than the radius of the bulk of most of the disk luminosity, and the photons approach from behind to first approximation (Dermer & Schlickeiser 1993). Thus, the Doppler factor of the plasmoid relative to the disk would be computed from Equation (1) with θ = 180° and β from Equation (30), ${\delta }_{d-p}=0.318$. We can then compute the photon energy density in the plasmoid from Equations (8) and (Lightman et al. 1975; Tucker 1975)

$$\begin{eqnarray}&&{U}_{\mathrm{ph}}\approx {\delta }_{d-p}^{4}\displaystyle \frac{{L}_{\mathrm{bol}}}{4\pi {\rm{c}}{D}_{d-p}^{2}}=3.93\times {10}^{-3}.\end{eqnarray} \tag{ 31 }$$

From Equations (21), (29), and (30), the range of Γ that results in EIC photons in the NuSTAR 2–30 keV band is 594 < Γ < 1170. These same leptons radiate synchrotron photons in the observer's frame in the far-infrared, $1.36\times {10}^{12}\,\mathrm{Hz}\lt {\nu }_{o}\lt 5.31\times {10}^{12}\,\mathrm{Hz}$. Based on the Epoch A fit in Table 3, these leptons radiate a synchrotron luminosity of ${L}_{\mathrm{synch}}=3.34\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ as observed at Earth. From Equations (29)–(31), the expected EIC luminosity observed at Earth in the NuSTAR band from these leptons is ${L}_{\mathrm{EIC}}=1.96\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ or a flux of ${F}_{\mathrm{EIC}}\,=4.48\,\times {10}^{-14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. This value is negligible compared to the fluxes in Table 4.

Note that given the large distance from the accretion disk to the plasmoid, ${D}_{d-p}\sim 2.12\times {10}^{17}$ cm, the broad-line region is likely an order of magnitude closer to the black hole and also illuminates the plasmoid primarily from behind (Greene & Ho 2007). Thus, the EIC flux is also weak, perhaps an order of magnitude less than that induced by the accretion disk photon field.

The primary controlling factor in the determination of the EIC from the dusty torus is the geometric dilution. Thus, the closest significant component in the complex IR spectrum is the most relevant. This is the hot dust on the inner face of the torus. The hot dust is typically considered to be the region in which graphite and silicate grains are sublimated with a temperature between ∼1500 K and ∼1800 K (Mor & Netzer 2012). It creates significant flux at wavelengths 1.6–2.2 μm. There is also "hot dust" that radiates a similar luminosity at 3 μm–4 μm in Mrk 231 with a temperature of ∼750 K (Lopez-Rodriguez et al. 2017). These components have the least geometric dilution of any of the IR components and contribute about >1/3 of the total AGN IR luminosity (Farrah et al. 2003).

Besides geometric dilution, one needs to know the relevant Doppler factors and that is based on the distribution of the hot dust. The primary piece of evidence is from the 3CRR catalog in which radio galaxy and radio-loud quasar IR spectra are compared. The radio galaxy IR spectra are heavily attenuated at wavelengths short of 5 μm. This implies that the hot dust is viewed through the optically thick torus (Hönig et al. 2011). The hot dust lies on the inner face of the torus and is restricted to a solid angle less than that subtended by the main body of the torus, less than ≈45° (Barthel 1989). The covering factor, CF, of the hot dust has been estimated at 0.15–0.35 (Mor & Netzer 2012). We assume 0.15 < CF < 0.35 and use this to compute an upper limit on the EIC flux in the calculation below. We note that the distribution of hot dust at the inner edge of the torus is restricted to within 45° of the equatorial plane.

The IR spectrum from data in Lopez-Rodriguez et al. (2017) at wavelengths shorter than 4 μm was fit with two blackbody components, T = 750 K as in Lopez-Rodriguez et al. (2017) and a hotter component T = 1650 K (Mor & Netzer 2012), with a luminosity of 1.06 × 10⁴⁵ erg s⁻¹ and 1.65 × 10⁴⁵ erg s⁻¹, respectively. The radius at which these hot dust distributions reside is $\approx \sqrt{0.35/\mathrm{CF}}6.5\times {10}^{18}$ cm and $\approx \sqrt{0.35/\mathrm{CF}}1.7\times {10}^{18}$ cm, respectively. We assume that the dust is isotropically distributed relative to the equator in the range $-45^\circ \lt \theta \lt +45^\circ$. If the hot dust is concentrated closer to the equatorial plane than 45°, then the Doppler factors are smaller and this calculation is an upper limit. For the hottest component,

$$\begin{eqnarray}\begin{array}{rcl}{U}_{\mathrm{ph}} & \approx & \displaystyle \frac{1}{c}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-{\theta }_{o}}^{{\theta }_{o}}\delta {\left(\beta =0.816,\theta \right)}^{4}\\ & & \times \,{I}_{{BB}}(T=1650\,{\rm{K}})d\phi \,d\theta =\displaystyle \frac{{CF}}{0.35}(1.16\times {10}^{-2}),\end{array}\end{eqnarray} \tag{ 32 }$$

where θ_o = 45°. From Equations (29), (30), and (32), the EIC flux in the NuSTAR band from the 1650 K dust is (for $0.15\lt \mathrm{CF}\lt 0.35$)

$$\begin{eqnarray}&&5.68\times {10}^{-14}\lt {F}_{\mathrm{EIC}}(T=1650\,{\rm{K}})\lt 1.33\times {10}^{-13}.\end{eqnarray} \tag{ 33 }$$

Similarly, for the other hot dust component,

$$\begin{eqnarray}&&2.42\times {10}^{-15}\lt {F}_{\mathrm{EIC}}(T=750\,{\rm{K}})\lt 5.55\times {10}^{-15}.\end{eqnarray} \tag{ 34 }$$

Note how much more inefficient the $T=750\,{\rm{K}}$ component is at inducing EIC emission in the plasmoid, as it is farther out. For completeness, we consider the remaining ∼3 × 10⁴⁵ erg s⁻¹ of AGN IR emission from the SED that is peaked at 60 μm (Farrah et al. 2003). Assuming a blackbody, this would require a radius of 60 pc at T ≲ 40 K (Farrah et al. 2003). The rotating molecular disk resolved with VLBI is roughly this size, 30–100 pc (Klöckner et al. 2003). It seems reasonable to assume that these are the same regions. The VLBI image reveals gas that is distributed in a region that is approximately orthogonal to the plasmoid velocity in Figure 5 as was the case for the hot dust (Klöckner et al. 2003; Lonsdale et al. 2003). There would be no strong Doppler enhancement for such a distribution. The geometric dilution is 10⁴ times that of the T = 1650 K dust with only twice the luminosity. Thus, there is no evidence of a pole on distribution of luminous dust associated with the AGN. VLBI observations revealed that the AGN dominates the starburst emission within ∼100 pc of the quasar (Lonsdale et al. 2003). If we assume the most extreme upper bound of placing all the starburst emission along the jet axis 100 pc away and recompute Equation (31) with θ = 0, we get ${U}_{\mathrm{ph}}\lt 3.93\times {10}^{-3}({10}^{4}){\left[\tfrac{2\times {10}^{17}\mathrm{cm}}{3\times {10}^{20}\mathrm{cm}}\right]}^{2}\lt 3\times {10}^{-5}$. This is ∼10⁻³ of the hot dust value in Equation (32). Of course, the starburst regions are probably more isotropically distributed, and the contribution would be significantly less than this upper limit. Therefore, we do not consider starburst emission as a viable seed field for X-ray emission. In summation, the combined EIC flux from Sections 8.1 and 8.2 is <8% of the observed continuum flux in Table 4.