Probing an X-ray flare pattern in Mrk 421 induced by multiple stationary shocks: a solution to the bulk Lorentz factor crisis

TL;DR

This study analyzes 13 years of X-ray data from Mrk 421 to identify a specific variability pattern linked to stationary shocks in the jet, providing insights into jet dynamics and resolving the bulk Lorentz factor crisis.

Contribution

It introduces a novel method to detect intrinsic variability patterns in blazar jets, linking radio stationary knots to high-energy emission and estimating jet perturbation properties without direct motion observation.

Findings

Identified a significant variability pattern consistent with a shock crossing all radio knots.

Estimated the jet perturbation's apparent speed at 45 c.

Linked stationary radio knots to high-energy particle acceleration.

Abstract

The common observations of multiple radio VLBI stationary knots in high-frequency-peaked BL Lacs (HBLs) can be interpreted as multiple recollimation shocks accelerating particles along jets. This approach can resolve the so-called "bulk Lorentz factor crisis" of sources with high Lorentz factor, deduced from maximum gamma-gamma opacity and fast variability, and apparently inconsistent slow/stationary radio knots. It also suggests that a unique pattern of the non-thermal emission variability should appear after each strong flare. Taking advantage of the 13 years of observation of the HBL Mrk 421 by the X-ray Telescope on the Neil Gehrels Swift Observatory Swift-XRT, we probe for such an intrinsic variability pattern. Its significance is then statistically estimated via comparisons with numerous similar simulated lightcurves. A suggested variability pattern is identified, consistent with a main flare emission zone located in the most upstream 15.3 GHz radio knot at 0.38 mas from the core. Subsequent flux excesses in the lightcurve are consistent with a perturbation crossing all the downstream radio knots with a constant apparent speed of 45 c. The significance of the observed variability pattern not arising from stochastic processes is found above 3 standard deviations, opening a promising path for further investigations in other blazars and with other energy bands. In addition to highlight the role of stationary radio knots as high-energy particle accelerators in jets, the developed method allows estimates of the apparent speed and size of a jet perturbation without the need to directly observe any motion in jets.

Tables

Table 1. Table 1: Projected distance from the radio core of the 4 VLBI quasi-stationary radio knots referenced by MOJAVE with their different associated names.

knot #					knot #	knot #	core distance	radius
(1)	(2)	(3)	[mas] (1)	[mas] (1)
Core	-	-	-	$4.24\pm 1.62\times {10}^{-2}$
8	C7	C4	$0.38\pm 0.07$	$1.20\pm 0.22\times {10}^{-1}$
11	C6	C3	$1.03\pm 0.16$	$2.04\pm 0.63\times {10}^{-1}$
9	C5	C2	$1.76\pm 0.29$	$3.66\pm 0.66\times {10}^{-1}$
6		C1	$3.96\pm 0.28$	$9.44\pm 0.21\times {10}^{-1}$
1: MOJAVE, 2: Piner et al. (2010), 3: Lico et al. (2012)

Table 2. Table 2: Parameter boundaries applied to the lightcurve models.

			parameter	boundaries	unit
Baseline	$f_0$	$[0,\infty ]$	erg cm-2 s-1
	$f_1$	$[-\infty ,\infty ]$	erg cm-2 s-1 day-1
	$\alpha$	$[0,\infty ]$	day4 erg cm-2 s-1
Multi-	$S_{core}$	$[0,3.8]$	–
Gaussian	$S_{knot}$	$[0,1.4]$	–
	$C$	$[0,50]$	day
	${\beta }_{app}$	$[2,70]$	c
	$h$	$[0,\infty ]$	erg cm-2 s-1
EMG	$\sigma$	$[0.5,50]$	day
	$\tau$	$[0.5,50]$	day

Table 3. Table 3: Fitting results for the different proposed models. The first column is for the time range presented in Figure 6 while the second column considers the range where at least one model is above 1 % percent 1 1\% of the baseline.

			${\chi }^2/dof$	${\chi }^2/dof$
	$t\in [-20,+100]$	$t\in [-7,+70]$
$EMG$	$7.86\times {10}^5/325$	$6.05\times {10}^5/233$
$G_m,core$	$7.16\times {10}^5/324$	$5.26\times {10}^5/232$
$G_m,knot$	$6.74\times {10}^5/324$	$4.83\times {10}^5/232$

Table 4. Table 4: Cut sets applied to select flares in the X-ray lightcurve. Systematic uncertainties are estimated from the loose and hard cuts.

Cuts					default	loose	hard	unit
$FT$	$90\%$	$80\%$	$95\%$	flux percentile
${\Delta }_{Flares}$	100	75	150	days
$\Delta t_{data}$	10	20	5	days

Table 5. Table 5: Parameter values for the different studied models with associated systematic uncertainties.

parameter			value	uncertainty
EMG
$f_0$	$\left(7.13\pm 0.27\right)\times {10}^{-1}$	${}_{-0.32}{}^{+0.40}\times {10}^{-1}$
$f_1$	$\left(-6.06\pm 0.90\right)\times {10}^{-4}$	${}_{-1.7}{}^{+2.5}\times {10}^{-4}$
$\tau$	$1.03\pm 0.42$	${}_{-0.45}{}^{}{}_{}{}^{+0}$
$\sigma$	$1.63\pm 0.21$	${}_{-0.25}{}^{}{}_{}{}^{+0}$
$h$	$3.77\pm 0.55$	${}_{-0.71}{}^{}{}_{}{}^{+0.64}$
core-flare model
$f_0$	$\left(6.67\pm 0.28\right)\times {10}^{-1}$	${}_{-0.25}{}^{+0.36}\times {10}^{-1}$
$f_1$	$\left(-4.80\pm 0.91\right)\times {10}^{-4}$	${}_{-1.4}{}^{+2.1}\times {10}^{-4}$
$\alpha$	$\left(-1.09\pm 0.25\right)\times {10}^2$	${}_{-0}{}^{+0.48}\times {10}^2$
$S$	$\left(4.1\pm 0.9\right)\times {10}^{-1}$	${}_{-1.3}{}^{+1.3}\times {10}^{-1}$
$C$	$1.56\pm 0.13$	${}_{-0.25}{}^{}{}_{}{}^{+0}$
${\beta }_{app}$	$30.3\pm 1.6$	${}_{-1.2}{}^{}{}_{}{}^{+3.9}$
knot-flare model
$f_0$	$\left(6.63\pm 0.27\right)\times {10}^{-1}$	${}_{-0.23}{}^{+0.36}\times {10}^{-1}$
$f_1$	$\left(-4.69\pm 0.90\right)\times {10}^{-4}$	${}_{-0.7}{}^{+1.9}\times {10}^{-4}$
$\alpha$	$\left(-1.22\pm 0.26\right)\times {10}^2$	${}_{-0}{}^{+0.53}\times {10}^2$
$S$	$\left(2.4\pm 0.6\right)\times {10}^{-1}$	${}_{-0.5}{}^{+1.1}\times {10}^{-1}$
$C$	$1.57\pm 0.15$	${}_{-0.39}{}^{}{}_{}{}^{+0.03}$
${\beta }_{app}$	$44.6\pm 1.2$	${}_{-0.3}{}^{}{}_{}{}^{+3.8}$

Table 6. Table 6: Post-simulation cuts to select lightcurves with enough selected flares and sufficiently high average fluxes M F , f ​ l ​ a ​ r ​ e ​ s subscript 𝑀 𝐹 𝑓 𝑙 𝑎 𝑟 𝑒 𝑠 M_{F,flares} to ensure a fair comparison with the Swift -XRT dataset. These cuts are adjusted for the three flares selections described in Section 5.4 .

				Default	Loose	Hard
Original dataset
nb flares	6	13	5
$M_{F,flares}$*	3.68	2.74	4.02
Simulations
Cut (nb flares)	$\ge 4$	$\ge 11$	$\ge 4$
Cut ($M_{F,flares}$)*	$\ge 3.3$	$\ge 2.5$	$\ge 3.5$
* Fluxes in ${10}^{-9}$ erg cm-2 s-1

Equations

$$\Delta t_i=(1+z)\frac{\Delta x_i}{c{\beta }_{app}}.$$

$$\text{If}(F_{i,j}-m)$$

$$F_{norm,i,j}$$

$$\Omega ={\left(\frac{G{M}_{\bullet }}{r}\right)}^{1/2}\frac{1}{r-r_g},$$

$$B(t)=f_0+f_{1}t.$$

$$t_i=7.26\times 10^2\frac{(1+z)x_i}{{\beta }_{app}}\mathrm{days}.$$

$${\sigma }_{G,i}=\frac{7.26\times 10^2}{\sqrt{2\log (2)}}\frac{(1+z)R_i}{{\beta }_{app}}S+C\mathrm{days},$$

$$P_i(t)=\frac{1}{{\sigma }_{G,i}\sqrt{2\pi }}\exp \left[\frac{-(t-t_i{)}^2}{2{\sigma }_{G,i}^2}\right],$$

$$A_i=\alpha /{\sigma }_{G,i}^3.$$

$$G_m(t)=\sum _{i=n}^5\left[A_i{P}_i(t)\right]+B(t),$$

$$EMG(t)=\frac{h\sigma }{\tau }\sqrt{\frac{\pi }{2}}\exp \left(\frac{{\sigma }^2}{2{\tau }^2}-\frac{t-\mu }{\tau }\right)$$

$$\times \mathrm{erfc}\left[\frac{1}{\sqrt{2}}\left(\frac{\sigma }{\tau }-\frac{t-\mu }{\sigma }\right)\right]+B(t),$$

$${\Phi }_{sim,LN}(t)=\exp \left(({\Phi }_{sim,norm}(t)\times a)+b\right),$$

$$F_{var}=\frac{\sqrt{V\left(\Phi (t)\right)-{\Delta }^2}}{⟨\Phi (t)⟩},$$

$$⟨\Phi (t)⟩=\exp (\mu +{\sigma }^2/2)$$

$$V\left(\Phi (t)\right)=(e^{{\sigma }^2}-1)⟨\Phi (t){⟩}^2$$

$$a=\sqrt{\left(\log (F_{var,obs}^2)+1\right)}.$$

$$b=⟨{\Phi }_{obs}(t)⟩-a^2/2.$$

$$median(\Phi (t))=e^{\mu }.$$

$${\sigma }_{result}={\mathrm{erf}}^{-1}(f_{sim})\sqrt{2}.$$

$$\theta$$

$$\delta$$

$$\Gamma$$

$$\Gamma =\frac{2\rho }{{\alpha }_{app}\sin \theta }.$$

$${\sigma }_{G,i}=\sqrt{{\sigma }_p^2+{\sigma }_{s,i}^2}.$$

$$W_p=2\sqrt{2\log (2)}\sqrt{{\sigma }_{G,i}^2-{\sigma }_{s,i}^2}\frac{c{\beta }_{app}}{1+z}.$$

$$0\le {\sigma }_{s,i}\le \frac{R_i}{\sqrt{2\log (2)}}\frac{1+z}{c{\beta }_{app}}.$$

$$W_{p,int}=\frac{W_p}{\Gamma \sin \theta }.$$

Full text

(galaxies:) BL Lacertae objects: individual (Markarian 421) — galaxies: jets — radiation mechanisms: non-thermal — acceleration of particles

††journal: ApJ††software: XSpec (Arnaud, 1996), astroML (Vanderplas et al., 2012), Astropy (Astropy Collaboration et al., 2013), SciPy (Jones et al., 2001–), NumPy (Walt et al., 2011), Matplotlib (Hunter, 2007))

1 Introduction

Multiwavelength studies of the variability and modeling of radio-loud AGN broadband SEDs attest to a compact emission zone moving with a high Lorentz factor close to the central engine. The particle individual Lorentz factors are often estimated to be above $10^{6}$ for the most energetic blazars, implying long-standing and powerful particle acceleration mechanisms. While the scenario of magnetic reconnection has received considerable attention during recent years, due to recent progress with MHD simulations (Sironi & Spitkovsky, 2014), the scenario of acceleration by shocks remains the most studied and the most accepted for the typical activity state of radio-loud AGN and their common variability (Marscher & Gear, 1985; Spada et al., 2001; Fromm et al., 2011).

The shock scenario is supported by multiple observations of gamma-ray flares in coincidence with the emergence of a jet perturbation (or overdensity) in or close to the radio core, mainly seen in flat-spectrum radio quasars (FSRQs) and some low- or intermediate-frequency-peaked BL Lacs (LBLs and IBLs) (Jorstad et al., 2001; Marscher et al., 2008; Abeysekara et al., 2018). The formation of recollimation shocks (also referenced as conical standing shock or reconfinement shock) in jets is also a phenomenon naturally observed in hydrodynamic and magnetohydrodynamic jet simulations as soon as a supersonic, or super-Alfvenic, non-pressured matched flow propagates through an external medium. This pressure mismatch at the interface between the jet inlet and the external medium generates two conical waves, namely a shock wave and a rarefaction wave. The shock wave propagates toward the external medium, and is reflected toward the jet axis as it reaches equilibrium with the external medium pressure. The rarefaction wave propagates toward the jet axis, locally dropping the jet pressure and accelerating the flow. The flow is then significantly slowed down after it reaches the reflection point of the conical waves at the jet axis. This process repeats and can produce a string of recollimation shocks until the full dissipation of energy carried out by the waves (e.g. Falle, 1991; van Putten, 1996; Gómez et al., 1997; Mizuno et al., 2015; Hervet et al., 2017).

Contrary to other blazar types, high-frequency-peaked BL Lacs (HBLs) show mainly stationary or low-speed VLBI radio features (radio knots) in their jets, in stark contrast to the high Lorentz factor values deduced from their variability or SED modeling (Hervet et al., 2016; Piner & Edwards, 2018). Most of the interpretations of this issue imply two distinct regions between radio knots and high-energy emission zones. Slow/stationary radio knots are assumed to come from a slower and wider jet part than the high energy emission zone. It can be understood as a strong jet deceleration very close to the core (Georganopoulos & Kazanas, 2003), or a stratified jet with differential speeds, as non-steady outflows (Lyutikov & Lister, 2010) or spine-layer structure (Ghisellini et al., 2005; Piner & Edwards, 2018). We adopt the interpretation of slow/stationary radio knots as a multiple recollimation shock structure, very stable for these sources due to their lower outer-jet kinetic power (Hervet et al., 2017).

Following the shock-in-jet model developed by Marscher & Gear (1985), a flare should happen when a perturbation (or moving shock) passes trough a recollimation shock. This scenario was adapted and improved by many further works and is quite successful as a picture of the the general broadband blazar flaring behavior (e.g. Komissarov & Falle, 1997; Türler et al., 2000; Türler, 2011; Nalewajko & Sikora, 2009; Nalewajko et al., 2012; Fromm et al., 2011, 2016; Marscher, 2014). Successive flares are then believed to be triggered by a stochastic injection from the central engine. However, while this approach assumes one shock at the base of the jet is responsible for the main dissipation process, it does not consider the other potential flares produced by downstream shocks. We investigate here the possibility of successive flares associated with successive recollimation shocks in relativistic jets. If we relate stationary radio knots to recollimation shocks, we can predict a distinct pattern of variability based on inter-knot gaps. Thus, after each strong flare occurring at the base of the jet, one should detect several other flares in accordance with the VLBI radio knot distribution in the jet, for a given velocity of the flow. The confirmation of such a pattern in HBL lightcurves would validate the role of stationary radio knots as high-energy particle accelerators, and characterize the apparent speed and size of underlying perturbations, extremely valuable for constraining the modeling parameters.

In Section 2 we introduce the basic concept of the proposed scenario and the ideal source for its application, Mrk 421. In Section 3 and 4 we describe how X-ray long-term lightcurves are handled in view of having the most efficient probe to detect a possible intrinsic post-flare variability pattern. The theoretical models used to check our scenario are developed in Section 5. In Section 6 we describe the method used to create simulated lightcurves as similar as possible to the real dataset, and also discuss biases induced by these simulations. Results and a general discussion are in Section 7.

Throughout this paper, a flat $\mathrm{\Lambda CDM}$ cosmology is adopted with $\mathrm{H_{0}}=69.7$ km s*-1* Mpc*-1*, $\Omega_{M}=0.286$, and $\Omega_{v}=0.714$ (Bennett et al., 2014). It leads to a projected scale of 1 mas = 0.603 pc at the redshift $z=0.030$ of Mrk 421.

2 Method and application to Mrk 421

2.1 Concept of the method

The core of the method is to probe flares associated with the flow passing through the knots, assuming they are stationary shocks. For a given apparent speed $\beta_{app}$, the time delay of the secondary flares can be set knowing the radio knot positions, as shown in Figure 1.

Considering a constant speed of the flow through a straight jet, the time gap $\Delta t$ between each successive flare in the lightcurve should be directly proportional to the observed inter-knot gap $\Delta x$. We have the relation

[TABLE]

Considering the association of radio knots with recollimation shocks, the underlying flow is expected to accelerate upstream of each shock due to the presence of rarefaction waves locally decreasing the pressure. The speed should then decrease after the shock. The realistic speed profile would be an oscillation, likely with a slower acceleration due to global conical opening of the jet (Komissarov & Falle, 1997; Gómez et al., 1997; Mizuno et al., 2015; Hervet et al., 2017). Throughout this paper we consider the approximation of an average constant speed of the underlying flow valid, with the main motivation keeping the lightcurve model developed in Section 5 as simple as possible. This approximation can be supported with the observed motions in radio jets, which in the majority are well fitted by a constant-speed motion (Lister et al., 2016). As further discussed in Section 5, the theoretical model developed also considers the width of the peaks from the size of the radio knots and a damping factor between successive flares.

2.2 Mrk 421: the ideal candidate

Mrk 421 is the brightest X-ray and gamma-ray HBL in the sky in its flaring and average state (Stroh & Falcone, 2013). It is one of the most monitored blazars in all wavelengths and shows frequent giant flares (e.g. Aleksić et al., 2015; Abeysekara et al., 2017; Fraija et al., 2017). Mrk 421 is perfectly adapted for this study by also presenting 4 well-defined VLBI quasi-stationary knots within 5 mas of the radio core at 15.3 GHz, as shown in Figure 2 (Left) from the MOJAVE collaboration.111http://www.physics.purdue.edu/MOJAVE All the observed knots show either non-radial or downward motions. Such motions would be very challenging to be described with a ballistic model, but can naturally match low amplitude shifts/oscillations of quasi-stationary recollimation shocks. The fastest measured knot measured in VLBI ($\#$6) displays an apparent speed of $0.217\pm 0.026$ c, roughly perpendicular to the jet direction (Lister et al., 2016), and the usual Doppler factor deduced from broadband spectral energy distribution (SED) modeling is about 20-25 (Błażejowski et al., 2005; Baloković et al., 2016; Carnerero et al., 2017; Kapanadze et al., 2018a, b), which can be seen as a lower limit, since the Doppler factor is usually constrained from the shortest variability timescale observed and from the maximum possible photon-photon opacity within the emitting region. For a canonical blazar angle with the line of sight of 2 deg, the SED models lead to a Lorentz factor $\Gamma_{model}\gtrsim 14$, which should be related to apparent downstream speed of $\beta_{app}\gtrsim 11$ c. Mrk 421 is then strongly affected by the bulk Lorentz factor crisis, which is ideal for our study.

For this study we consider these 4 knots as stationary recollimation shocks with their distance to the radio core given by the mean value of the measured distances from the MOJAVE Collaboration. The uncertainty on their distance to the core and radius are given by the standard deviation of the dataset. The Mrk 421 knot string follows a conical expansion well, as shown in Figure 2 (Right). The knots’ radius is fitted by a linear function $f(x)=(0.195\pm 0.015)x+(3.94\pm 0.76)\times 10^{-2}$ mas, with a reduced $\chi^{2}$ of 0.28. The radio knot positions of Mrk 421 were measured in several other studies for different frequencies and epochs. Although the MOJAVE dataset is the one the most simultaneous with the lightcurve in our study, it remains relevant to check the consistency of these measurements with the previous observations described in Piner et al. (2010) (with extended dataset from Piner et al. (1999); Piner & Edwards (2005)), and Lico et al. (2012).