Core shift effect in blazars

TL;DR

This study investigates the core shift effect in blazars using long-term radio observations, deriving physical parameters, and estimating black hole spins, revealing insights into jet physics and accretion processes.

Contribution

The paper provides a comprehensive analysis of the core shift effect in three blazars over 40 years, estimating magnetic fields, core positions, and black hole spins using multi-frequency radio data.

Findings

Core shift effect varies among the blazars studied.

Magnetic field strengths are consistent with previous estimates.

Black hole spins are estimated to be between 0.15 and 0.9.

Abstract

We studied the pc-scale core shift effect using radio light curves for three blazars, S5 0716+714, 3C 279 and BL Lacertae, which were monitored at five frequencies ($\nu$) between 4.8 GHz and 36.8 GHz using the University of Michigan Radio Astronomical Observatory (UMRAO), the Crimean Astrophysical Observatory (CrAO), and Metsahovi Radio Observatory for over 40 years. Flares were Gaussian fitted to derive time delays between observed frequencies for each flare ($\Delta t$), peak amplitude ($A$), and their half width. Using $A \propto \nu^{\alpha}$ we infer $\alpha$ in the range $-$16.67 to 2.41 and using $\Delta t \propto \nu^{1/k_r}$, we infer $k_r \sim 1$, employed in the context of equipartition between magnetic and kinetic energy density for parameter estimation. From the estimated core position offset ($\Omega_{r \nu}$) and the core radius ($r_{\rm core}$), we infer that opacity model may not be valid in all cases. The mean magnetic field strength at 1 pc ($B_1$) and at the core ($B_{\rm core}$), are in agreement with previous estimates. We apply the magnetically arrested disk model to estimate black hole spins in the range $0.15-0.9$ for these blazars, indicating that the model is consistent with expected accretion mode in such sources. The power law shaped power spectral density has slopes $-$1.3 to $-$2.3 and is interpreted in terms of multiple shocks or magnetic instabilities.

Equations

$$\tau (r)\propto r^{(1.5-\alpha )m+n-1}{\nu }^{2.5-\alpha }$$

$$k_r=\frac{(3-2\alpha )m+2n-2}{5-2\alpha }$$

$$\Delta r=r_{{\nu }_1}-r_{{\nu }_2}=\frac{({\nu }_1-{\nu }_2)}{{\nu }_1{\nu }_2}.$$

$$\Delta \theta =\mu \Delta t.$$

$${\Omega }_{r\nu }=(1\mathrm{mas})\frac{D_L\Delta \theta }{(1+z{)}^2\left({\nu }^{-1/k_r}-{\nu }_0^{-1/k_r}\right)},$$

$$r_{\mathrm{core}}=\frac{{\Omega }_{r\nu }}{\sin i}{\nu }^{-1/k_r}.$$

$$B_{\beta }={\left(\frac{(h_P(1+z){)}^{\frac{p+4}{p+2}}}{\delta \sin i}\frac{B_{cr}}{m_e{c}^2}{\left(\frac{8{\alpha }_F}{(1\mathrm{pc})C_2(p){\sigma }_T\tan {\theta }_0}\right)}^{\frac{2}{p+2}}\right)}^{\frac{p+2}{p+6}}$$

$$\frac{1}{h_1{\beta }^{\frac{2}{p+6}}}{\left(\frac{{\Omega }_{r\nu }}{10^{-9/k_r}}\right)}^{\frac{p+4}{p+6}},$$

$$B_F=\frac{\delta (1+z{)}^7}{D_L^4{\sin }^{3}i}\frac{B_{cr}^5{h}_P^7}{c^{12}}{\left(\frac{{\alpha }_F{C}_1(p)C_3(p)\tan {\theta }_0}{(1\mathrm{Jy})24{\pi }^3{C}_2(p)}\right)}^2$$

$$\frac{1}{h_1{F}_{\nu }^2}{\left(\frac{{\Omega }_{r\nu }}{10^{-9/k_r}}\right)}^5.$$

$$B_{\mathrm{core}}=B{r}_{\mathrm{core}}^{-1},$$

$$B_{\mathrm{pj}}=\frac{Bc{\Gamma }_{\mathrm{j}}}{{\Omega }_{\mathrm{f}}{r}_{\mathrm{j}}}{\left(\frac{2\sigma }{1+\sigma }\right)}^{1/2}$$

$$B=B_{\varphi \mathrm{j}}^{\prime }{\left(\frac{1+\sigma }{2\sigma }\right)}^{1/2},$$

$${\Phi }_{\mathrm{j}}\equiv 2\pi {\int }_0^{r_{\mathrm{j}}}\mathrm{d}rr{B}_{\mathrm{p}}(r)=\frac{2}{2-\alpha }\pi B_{\mathrm{pj}}{r}_{\mathrm{j}}^2.$$

$${\Phi }_{\mathrm{j}}=\frac{\pi B_{\mathrm{pj}}{r}_{\mathrm{j}}^2}{\sigma }.$$

$$l=\frac{{\Omega }_{\mathrm{f}}}{{\Omega }_{\mathrm{H}}}$$

$${\Phi }_{\mathrm{j}}=\frac{\pi Bc{r}_{\mathrm{j}}{\Gamma }_{\mathrm{j}}}{l{\Omega }_{\mathrm{H}}\sigma }{\left(\frac{2\sigma }{1+\sigma }\right)}^{1/2}$$

$${\Omega }_{\mathrm{H}}=\frac{ac}{2r_{\mathrm{H}}},$$

$${\Phi }_{\mathrm{j}}=\frac{\pi 2^{3/2}B{r}_{\mathrm{j}}{r}_{\mathrm{H}}{\Gamma }_{\mathrm{j}}}{\sigma la}{\left(\frac{\sigma }{1+\sigma }\right)}^{1/2},$$

$$r_{\mathrm{j}}=h{\theta }_{\mathrm{j}},\mathrm{and}{\theta }_{\mathrm{j}}{\Gamma }_{\mathrm{j}}=s{\sigma }^{1/2}.$$

$${\Phi }_{\mathrm{j}}=\frac{\pi 2^{3/2}hB{r}_{\mathrm{H}}s}{la}{\left(\frac{1}{1+\sigma }\right)}^{1/2}$$

$$\frac{B}{B_{\varphi \mathrm{j}}^{\prime }}={\left(\frac{1+\sigma }{2\sigma }\right)}^{1/2}$$

$$B_{\mathrm{pj}}=\frac{B_{\varphi \mathrm{j}}^{\prime }c{\Gamma }_{\mathrm{j}}}{{\Omega }_{\mathrm{f}}{r}_{\mathrm{j}}}.$$

$${\Phi }_{\mathrm{BH}}={\varphi }_{\mathrm{BH}}(\dot{M}c{)}^{1/2}{r}_{\mathrm{g}},$$

$$\frac{\pi 2^{3/2}hBs}{la}{\left(\frac{1}{1+\sigma }\right)}^{1/2}{r}_{\mathrm{g}}[1+(1-a^2{)}^{1/2}]={\varphi }_{\mathrm{BH}}(\dot{M}c{)}^{1/2}{r}_{\mathrm{g}}.$$

$$[\frac{1+(1-a^2{)}^{1/2}}{a}{]}^2{h}^2{B}^{2}c\left(\frac{1}{1+\sigma }\right)=qL,$$

$$q=\frac{l^2{\varphi }_{\mathrm{BH}}^2(1+\sigma )}{s^2ϵ2^3{\pi }^2}$$

$$w(a)=\frac{qL}{h^2{B}^{2}c}=q{L}_{46}(\frac{0.2\mathrm{Gpc}}{hB}{)}^{-2},$$

$$w(a)=[\frac{1+(1-a^2{)}^{1/2}}{a}{]}^2,$$

$$a(w)=\frac{2\sqrt{w}}{w^2+1}.$$

$$y=A{e}^{[-(t-\overline{m}{)}^2/(2{\sigma }^2)]},$$

Full text

††pagerange: Core shift effect in blazars–References

1 Introduction

Flat spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lac) together known as blazars are a subset of core dominated Active Galactic Nuclei (AGNs) characterized by a luminous core, rapid variability over entire electro magnetic (EM) spectrum, high radio to optical polarization, apparently superluminal jet components and non thermal emission from a Doppler boosted relativistic jet pointing $\leq$ 10*∘* to the observer line of sight (LOS). Core-jet morphology is a common characteristic of most of these AGN in very long baseline interferometry (VLBI) images, where core is the optically thick base of the jet (e.g. Blandford & Königl 1979).

The apparent systematic outward shift of the VLBI core position with decreasing observation frequency is the core shift effect, and can be attributed to the synchrotron self absorption process (SSA) operating and causing the transition from optically thick to thin jet emission, the optical depth given by,

[TABLE]

where $r$ is the radial distance from the center, $m$ is the power law index in the magnetic field decay ($B=B_{1}(r/1\rm{pc})^{-m}$) while $n$ is the power law index in the particle number density decay ($N=N_{1}(r/1\rm{pc})^{-n}$) with distance. Here $N_{1}$ and $B_{1}$ refer to $N$ and $B$ at $r=1$ pc where it is expected that the jet shape changes from conical to parabolic closer to the true jet origin (e.g., Beskin & Nokhrina 2006; McKinney 2006; Krichbaum et al 2006; Nakamura & Asada 2013). In the conical jet model, the VLBI core of blazars are compact, stationary and bright features located at one end of the jet at sub-mas–mas scales. They are characterized by superluminal jet components which emerge and separate and whose kinematics is governed by propagating disturbances or shock waves. The core is partially optically thick to SSA since the spectral slope is $\leq$ 1, which can be attributed to magnetic field geometry and the relativistic electron density gradients. Cawthorne (2006) suggests that in some cases, the core could correspond to a stationary feature like a conical shock, inferred in some long term studies of blazar jets (e.g. Marscher 2009).

Core shift measurements can provide information on the physical parameters of the VLBI jet including core magnetic field strength, distance from the core to the base of the jet (e.g. Lobanov 1998; Hirotani 2005), spectral index, inputs to rotation measure maps, and clarify the nature of the absorbing material in the SSA jet. There are several ways to measure apparent core-shifts. One of them is the phase referencing technique where the telescope is switched between the target source and a nearby phase calibrator, the switching time being shorter compared to the coherence time. However as this technique is complex, resource intensive and relatively expensive in terms of the observing time, it has been applied only for some AGN e.g., 3C 390.1, 3C 395, 1038+528, M81, 4C 39.25 (Guirado et al. 1995; Marcaide & Shapiro 1984; Lara et al. 1994; Bietenholz, Bartel, & Rupen 2004). Other methods require a brightness distribution model where the core position is estimated by identifying individual optically thin jet features in the VLBI images which do not change their position with frequency (e.g. Kovalev et al. 2008; Sokolovsky et al. 2011). Though, this method requires simultaneous multiband VLBI observations. Another method involves calculation of the time delay between the lower and higher frequency radio emission which can be attributed to the opacity effects. It is applied assuming that at a given frequency, a flare is associated with jet component ejection from the core. The observation time of the flare and hence the shift increases with decreasing frequency (e.g. Kudryavtseva et al. 2011; Mohan et al. 2015). Magnitudes of core shifts are time dependent as they depend on the present state of an AGN and are resolution limited for the VLBI observations. Relative advantages of this method include the use of only single dish measurements with quasi-simultaneous observations, the ability to use VLBI core flux density measurements if available, and the available methodology for the timing analysis (Mohan et al. 2015). Shortcomings include a possible inability to capture a flare due to non-simultaneous observations, the difficulty in identifying a flare in two frequencies as the same (opacity effect) or distinct, and the requirement of more than two frequencies to confirm an increasing time delay and hence, the opacity based apparent core shift. As multiple quasi-simultaneous observations are employed in our study, these shortcomings are met thus making it ideally suited.

We discuss the time delay based core shift estimation and the resulting jet parameters for three blazars, S5 0716+714, 3C 279, and BL Lac using single dish observations at five frequency bands from University of Michigan Radio Astronomical Observatory (UMRAO), Crimean Astrophysical Observatory (CrAO) and Metsähovi Radio Observatory, observed between 1970 and 2015. The estimated parameters are compared with estimates from previous studies. We adopt a standard cold dark matter cosmological model with Hubble constant $H_{0}=71~{}$km/s/Mpc, matter energy density $\Omega_{m}=0.27$, and dark energy density $\Omega_{\Lambda}=0.73$. In Section 2, we briefly discuss the source properties relevant to the current study; in Section 3, the parameters that are calculated including core position offset, magnetic field strength, and size of the emitting core are presented, the estimates of which can be employed in the estimation of black hole spin in the context of the magnetically arrested disk (MAD) scenario; a description of observations and data reduction for our single dish observations along with the analysis procedure is then presented in Section 4. This includes a piecewise Gaussian fitting and the variability analysis using the Fourier periodogram technique. In Section 5, we present the results of our analysis and their comparison with reports in literature. The implications and a brief summary of our results is then presented in Section 6.

2 Notes on Individual Sources

2.1 S5 0716+714

S5 0716+714 is a bright, high declination ($\alpha_{2000.0}$ = 07h 21m 53.4s, $\delta_{2000.0}$ = $+71^{\circ}20^{{}^{\prime}}36.4^{\prime\prime}$) BL Lacertae object. Based on the marginal detection of the host galaxy, a redshift of 0.31 $\pm$ 0.08 was reported by Nilsson et al. (2008). Recently, Paiano et al. (2017) obtained a lower limit of z $>$ 0.10 during the high optical state of the source. S5 0716+714 is a prominent blazar in optical bands and has shown variability on diverse timescales from minutes through hours to days (e.g. Heidt & Wagner 1996; Nesci et al. 1998; Giommi et al. 1999; Raiteri et al. 1999; Raiteri et al. 2003; Gupta et al. 2008; Gupta et al. 2012; Agarwal et al. 2016) in radio to X ray bands (e.g. Wagner et al. 1996; Raiteri et al. 2003 and references therein). It was included in the S5 catalog after its discovery in the Bonn-NRAO Radio 5 GHz Survey (Kühr et al. 1981).

A lower limit on the apparent brightness temperature, $T_{B}$ $\geq$ 2 $\times$ 1012 K was obtained using 5 GHz VLBI observations and suggests a minimum equipartition Doppler factor of $\delta$ $\geq$ 4 (Bach et al. 2006). In addition, high bulk Lorentz factors of $\geq$ 16$-$21 were estimated based on the inferred jet speeds (Jorstad et al. 2001). Radio maps from VLBI reveal a compact core-jet structure as has been suggested for BL Lacertae with different components moving at different velocities (Bach et al. 2003). Jet speed $>$ 11-15 h*-1c were estimated from the component proper motions Jorstad et al. (2001), who also found quasi-periodic ejection of components every $\sim$ 0.7 yr. An analysis of the above 26 VLBI observations in 5-22 GHz frequency band gives an apparent velocity in range 5-10 h-1*c. Radio flux density and linear polarization light curves obtained from UMRAO at 4.8, 8.0, and 14.5 GHz were presented by Aller et al. (1985) for the 1981-1984 period and those during 1985-1992 were studied by Wagner et al. (1996). Later, Teraesranta et al. (1998) presented this target’s radio data till 1999 which included 22 and 37 GHz observations from MRO.

2.2 3C 279

The FSRQ 3C 279 ($\alpha_{2000.0}$ = 12h 56m 11.17s, $\delta_{2000.0}$ = $-05^{\circ}47^{{}^{\prime}}21.5^{\prime\prime}$) at a redshift of 0.536 (Burbidge & Rosenburg 1965) was the first extragalactic radio source indicating superluminal motion (Cohen et al. 1971). It displays violent multiwavelength flux variability and is well studied by various observation programs (e.g. Maraschi et al. 1994; Wehrle et al. 1998; Lindfors et al. 2006, Collmar et al. 2007). Flux variability at X-rays, R band and 14.5 GHz frequency band have been found to be significantly correlated, suggesting a possible co-spatial origin. The radio morphology reveals a bright, stationary core with a jet extending to $\sim$ 5$\aas@@fstack{\prime\prime}$ along a position angle of 205*∘* (de Pater & Perley 1983). Superluminal motion of the jet with apparent speed in range 5$c$ - 17$c$ was found to be associated with the target using very long baseline array (VLBA) radio observations at 43 GHz between 1998 March and 2001 April (Jorstad et al. 2004). This source has been extensively monitored with the high resolution VLBI revealing a one sided jet extending southwest on pc-scale characterized by bright knots ejected from the core region, along with multiple projected apparent speeds and polarization angles (e.g. Jorstad et al. 2005; Unwin eta al. 1989; Chatterjee et al. 2008; Wehrle et al. 2001). Using VLBA polarimetry studies, electric field vector was found to be aligned with the jet direction on pc to kpc-scales implying that the magnetic field is pre-dominantly perpendicular to the relativistic jets on above length scales. From VLBA radio observations at 43 GHz between 1998 March - 2001 April, a bulk Lorentz factor of 15.5 $\pm$ 2.5, viewing angle of 2*∘.1 $\pm$ 1∘*.1 and Doppler factor of 24.1 $\pm$ 6.5 are obtained (Jorstad et al. 2004, 2005).

2.3 BL Lacertae

BL Lacertae ($\alpha_{2000.0}$ = 22h 02m 43.29s $\delta_{2000.0}$ = $+42^{\circ}16^{{}^{\prime}}39.98^{\prime\prime}$) also known as IES 2200+420, is a blazar archetype hosted by an elliptical galaxy at a redshift of $\sim$ 0.07 (Miller, French, & Hawley 1978; Hyvönen et al. 2007). Emission lines are either absent or extremely weak in BL Lacertae class of blazars. But broad Hα and Hβ lines have been found in the spectrum of this object rendering its classification in BL Lacertae class under delima (Vermeulen et al. 1995). Since its association with radio source VRO 42.22.01 (Schmitt, 1968), it has been well studied across the entire EM spectrum. It has been extensively observed in the optical band at diverse timescales i.e. from minutes to years timescales (e.g. Villata et al. 2004; Agarwal & Gupta 2015).

From $\sim$ 15 years of variability observations (June 1971 to January 1985), Webb et al. (1988) inferred recurrent variations after every 0.31, 0.60, and 0.88 yr using the Fourier periodogram. Marchenko et al. (1996) found a statistically significant long term component with $P$ = 7.8 $\pm$ 0.2 yr in the 20 year duration light curve using the whitening of the time series method. Smith & Nair (1995) found a period of 7.7 yr using a 20 yr duration light curve using Fourier analysis thus, consistent with Marchenko et al. (1996). Later, Fan et al. (1998) obtained a long term period of 14 yr using the Jurkevich method on the optical light curve. The inclination angle for the Doppler boosted approaching jet of the BL Lacertae is $\sim$ 6*∘-10∘* with a flow speed of 0.981-0.994 c and a bulk Lorentz factor of 7.0 $\pm$ 1.8 (Jorstad et al. 2005).

3 Frequency-dependent core shifts and tests of the magnetically arrested disk (MAD) scenario

3.1 Core shift and jet parameters

The optically thick radio core is not completely resolved with VLBI and thus the emission from the core region is mainly from the $\tau$ = 1 surface at a particular frequency. Due to its frequency dependence, the radio core position follows $r\propto\nu^{-1/k_{r}}$ (e.g. Konigl 1981; O’Sullivan & Gabuzda 2009) where $r$ is the distance from the central region, $\alpha$ is the spectral index in the power law flux density dependence $S_{\nu}\propto\nu^{\alpha}$, and the index $k_{r}$ is

[TABLE]

assuming that the ambient medium pressure on the jet is negligible and external pressure is non-negligible with non-zero gradient along the jet (Lobanov 1998).

The quantity $k_{r}$ becomes independent of $\alpha$ if there is equipartition between the jet particle density and magnetic field energy density for which $m=1$ and $n=2$ (e.g. Hutter & Mufson 1986; Lobanov 1998). The above combination of $m$ and $n$ is valid for synchrotron emission from compact VLBI cores (Konigl 1981) for a constant jet speed and half opening angle. With $k_{r}=1$, the core shift between two frequencies $\nu_{1}$ and $\nu_{2}$ is

[TABLE]

In the models proposed by Lobanov (1998) and Hirotani (2005), magnetic field at some given distance along the jet frame is derived without using the information about the radiative flux of the jet. The magnetic field calculated then would depend on the normalization of electron distribution. These works assumed equipartition between the relativistic electrons and the magnetic field. For S5 0716+714, 3C 279, and BL Lacertae we use the modified equations following Zdziarski et al. (2015) which improves on the method previously followed in the application to 3C 454.3 in (Mohan et al. (2015); Paper 1). The angular core separation $\Delta\theta$ (mas) for a mean component proper motion $\mu$ (mas/y) is

[TABLE]

The core offset parameter $\Omega_{r\nu}$ (pc GHz${}^{1/k_{r}}$) which is the core offset distance per unit observation frequency $\nu^{1/k_{r}}$ (GHz) interval is

[TABLE]

where $D_{L}$ is the luminosity distance. The core offset distance $r_{\rm core}$ (pc) is the projected distance of the actual emitting core (at a frequency $\nu_{0}\to\infty$) at the observation frequency $\nu$ (GHz) and for a jet inclination angle $i$ is

[TABLE]

The magnetic field strength $B_{\beta}$ (G) at a distance $h_{1}=h/(1{\rm pc})$ is

[TABLE]

where $p=2\alpha+1$ is the particle energy distribution index, $h_{P}$ is Planck’s constant, $m_{e}$ is the electron mass, $e$ is its charge, $B_{cr}=2\pi m^{2}_{e}c^{3}/(eh_{P})$, $\alpha_{F}=1/137$ is the fine structure constant, $\beta$ is the ratio between the particle kinetic energy density and the magnetic field energy density and the constants $C_{1}$, $C_{2}$ and $C_{3}$ are defined in terms of $p$ (Zdziarski et al. 2012). The magnetic field strength $B_{F}$ (G) at the distance $h_{1}$, accounting for an optically thin jet based flux $F_{\nu}$ (includes jet and core emission) is,

[TABLE]

The magnetic field strength $B_{\rm core}$ (G) at the emitting core at the distance $r_{\rm core}$ (pc) is

[TABLE]

where $B$ can be calculated by determining $B_{\beta}$, $B_{F}$ and the ratio $B_{\beta}/B_{F}$ from their means, used as a normalizing multiplicative factor.

3.2 Application of the MAD model

Magnetic fields are expected to play a prominent role in structuring the accretion flow, the disk-jet connection, and the collimated outflow at the sub–pc to pc-scales and in possible helical signatures at the kpc-scales. Further, using the derived field strength, the bolometric and jet luminosity one can explore electrodynamical and hybrid jet models and place constraints on the spin and mass of the black hole.

Zamaninasab et al. (2014) considered samples of blazars and radio galaxies and obtained core shift and magnetic field strengths. They considered a model of jet formation from black hole (BH) spin-energy extraction (Blandford & Znajek, 1977). Zdziarski et al. (2015) considered models with the accretion being magnetically arrested (MAD; Narayan, Igumenshchev, & Abramowicz (2003)). Such flows have dragged so much magnetic flux to the BH that it becomes dynamically important owing to obstruction of accretion, and is given by $\Phi_{BH}\simeq 50(\dot{M}c)^{1/2}=10^{-4}M_{8}^{3/2}\dot{m}^{1/2}$ pc2 (Tchekhovskoy, Narayan, & McKinney, 2011; McKinney, Tchekhovskoy, & Blandford, 2012) where the accretion rate is in units of the Eddington rate. The surface value of magnetic flux depending on $h$ is (Tchekhovskoy et al. 2009),

[TABLE]

where $B$ is the transverse average magnetic field strength. To compare with the observations, we use the average value of $B^{\prime}_{\phi}$ (Zdiarski et al. 2015),

[TABLE]

where $\sigma$ is the ratio of the magnetic to particle energy and the magnetic flux in terms of the mean toroidal field $B$ is

[TABLE]

Using $\alpha=2/(1+\sigma)$ where $\alpha$ is the radial index of the variation of the poloidal field,

[TABLE]

The ratio of the angular frequency of the field lines to the BH angular frequency is defined by a parameter $l$ as

[TABLE]

where $r_{\rm H}=[1+(1-a^{2})^{1/2}]r_{\rm g}$ is the BH horizon radius. Now using (10) and (14) we find

[TABLE]

and using

[TABLE]

we obtain

[TABLE]

where

[TABLE]

Applying eqn (18) in (17), we derive

[TABLE]

where $a$ is the dimensionless spin parameter. Now,

[TABLE]

and inserting (20) in (10) we get,

[TABLE]

According to Narayan et al. (2003) the poloidal flux treading the BH on one hemisphere is given as:

[TABLE]

If we equate $\Phi_{\rm j}$ with $\Phi_{\rm BH}$ applying the relation $r_{\rm H}=[1+(1-a^{2})^{1/2}]r_{\rm g}$, we arrive at

[TABLE]

Putting $L=\epsilon\dot{M}c^{2}$, we find

[TABLE]

where

[TABLE]

is a dimensionless parameter and from equation (24),

[TABLE]

where, the quantity

[TABLE]

and the bolometric luminosity $L_{46}$ is in units of $10^{46}$ erg/s. For typical values of the parameters, $s=0.5,l=0.5,\epsilon=0.4,\phi_{\rm BH}=50,\sigma=0.5$, we find that $q\simeq 120$. Inverting the eqn (27), we obtain

[TABLE]

The function $a(w)$ is plotted in Fig. 13. It predicts that the values of $0.296<w<1$ are disallowed since they yield spin values $a>1$. The relation given in Figure 2 of Zamaninasab et al. (2014) and Zdziarski et al. (2015) also validate a variant of eqn. (26) with the choice of $\phi_{BH}=50$ and $\epsilon=0.4$. We have adopted these values and computed $w$ and $q$ to obtain approximate values of spin consistent with known luminosity and variability in these sources. For each source, we have calculated a linear fit of log(Bcore) against log(rcore). The corresponding fit values of the intercept were used to calculate the spin values for each source. The luminosity is expected to be peaked at an intermediate $a$ due to the fact that in the BZ model, the horizon traps smaller flux for larger spin as the horizon shrinks with $a$ but at lower spin, rate of change of flux which is proportional to the spin reduces. The spin in combination with the varying $Bh$ in the sources, can explain the jet luminosity within the context of the MAD model. We test the model predictions using the core shift measurements for the three sources in section 5.2.

4 Observations and analysis

The long term light curves of blazars are obtained from the University of Michigan Radio Astronomical Observatory (UMRAO) which has monitored compact variable radio sources spanning over $\sim$ 40 years with a 26 m paraboloid dish and provides multifrequency (4.8, 8.0 and 14.5 GHz), high time resolution (better than a month) light curves (Aller et al. 1985, 1999). Observations at 4.8 GHz started in 1978 on regular basis, at 8.0 GHz in 1965 and at 14.5 GHz in 1974. Monitoring at 22.2 GHz and 36.8 GHz were carried out with the 22-m radio telescope (RT-22) of the Crimean Astrophysical Observatory (CrAO; Volvach 2006). We use modulated radiometers in combination with the registration regime “ON-ON” for collecting data from the telescope (Nesterov, Volvach, & Strepka 2000). The 14 m radio telescope of Aalto University Metsähovi Radio Observatory in Finland was used for observations at 37.0 GHz. Data obtained at Metsähovi and RT-22 were combined in a single array to supplement each other. A detailed description of the data reduction and analysis of Metsähovi data is presented in Teraesranta et al. (1998).

Fig. 1 shows the total flux density light curves for S5 0716+714 at five frequencies from 4.8 GHz to 36.8 GHz while Fig. 2 displays the same for 3C 279 and Fig. 3 for BL Lacertae, spanning more than four decades of observations.

The method employed here was applied to the FSRQ, 3C 454.3, in paper 1 with the results being consistent with previous studies. The radio light curves of S5 0716+714, 3C 279 and BL Lacertae are piecewise Gaussian fitted as described in paper 1 using

[TABLE]

where $y$ is the typical flare profile, $A$ is the flare amplitude, $\overline{m}$ is the peak position and $\sigma$ is the flare width. Following the procedure described in Paper 1. Initially pre-processing was performed to locate the maxima and minima of the flaring regions. Gaussian filtering was applied from the baseline at the minimum amplitude of the LC with some initial values of maxima value, maxima position and half the difference between consecutive minima ($A_{0}$, $\overline{m}_{0}$ and $\sigma_{0}$) respectively. Trial values of $A$, $\overline{m}$ and $\sigma$ were generated from the regions around the initial values. These trial values account for possible values of the flare peak amplitude, position, and difference in the original LC. Single Gaussian fit procedure described in Paper 1 is cycled through the entire LC covering all flares till the best $\chi^{2}$ fit value was obtained.

The light curves, best fit Gaussians and the residuals for S5 0716+714 are displayed in Figure 4, while for BL Lacertae they are presented in Fig 7. For 3C 279, we split the light curve into two segments due to the difference in the flaring amplitudes. Segment 1 is from the beginning of observations till 1997 covering low amplitude flares compared to larger amplitude flares in Segment 2 from 1998 to the end of the observation. Due to flaring behaviour being different during these epochs, applying the Gaussian fits to the entire light curve can lead to loss in accounting for prominent features (even though low amplitude) from the Segment 1 thus causing incorrect results. For the light curves of Segment 1 of 3C 279, the Gaussian fitting and the residuals are presented in Fig. 5. The Gaussian fits along with the LC and residual (in the bottom panel) for the second segment of 3C 279 are shown in Figure 6 for five frequency bands. In all these figures, black points show the data points while the solid blue curve is the sum of all fitted Gaussians, covering the main features in the light curve. Residuals plotted in the lower panel of these figures are calculated as $[y(A_{i},\overline{m}_{i},\sigma_{i})-x_{i}(t_{i})]/{\rm Standard~{}Deviation(x_{i}(t_{i}))}$.

The parameters derived from the Gaussian fit procedure are given in Tables 1$-$4. Columns 1–7 give the flare nomenclature, observing frequency, maximum amplitude of the flare ($A$), epoch of maximum flux, full width at half maximum (FWHM) from the Gaussian fit, time delays and spectral index calculated using the convention $A\propto\nu^{\alpha}$ for those flares present in at least three frequencies.

The periodogram analysis employed here is similar to that carried out for the analysis of X-ray, radio and optical light curves from Seyfert galaxies and blazars in (Mohan & Mangalam 2014; Mohan et al. 2015; Mohan et al. 2016). It is employed in the determination of the power spectral density which is emergent from various physical process causing the variability in the light curve (jet and accretion disk based) and effects due to the sampling process during the accumulation of the light curve. It is also used to infer any possible statistically significant quasi-periodic features in the light curves originating from the above physical processes. Here, the 4.8 - 36.8 GHz light curves were interpolated and re-sampled at regular intervals of $\Delta t$ = 0.1 d for the time series analysis with the normalized Fourier periodogram,

[TABLE]

where $\Delta t$ is the sampling time step for the evenly sampled light curve $x(t_{n})=x(n\Delta t)$, and $|F(f_{j})|$ is its discrete Fourier transform evaluated at frequencies $f_{j}=j/(N\Delta t)$ with $j=1,2,..,(N/2-1)$. The analysis procedure includes the periodogram estimation, parametric model fits to determine the PSD shape, model selection and statistical significance testing and is discussed in more detail in the above mentioned work and references therein.

5 Results and discussion

5.1 Core shift effect

The highest frequency at which the flare occurred has been used as the reference frequency for measuring time lags ($\Delta t$) and are plotted against observation frequency in Figure 8 for S5 0716+714, in Figs. 9 and 10 for 3C 279 and Figs. 11 and 12 for BL Lacertae. They generally indicate a power law and were thus fitted with a function $\Delta t=a\nu^{-1/k_{r}}+b$ using a weighted non-linear least squares method to calculate $a$, $b$, and $k_{r}$. The median of each parameter and the associated median deviation are

For S5 0716+714

b = $-$0.03 $\pm$ 0.02   yr

a = 1.01 $\pm$ 0.08   yr (GHz$)^{-1/k_{r}}$

kr = 1.03 $\pm$ 0.09 . 2. 2.

First segment of 3C 279: beginning of observation – 1997.0

b = $-$0.30 $\pm$ 0.04   yr

a = 340.2 $\pm$ 48.3   yr (GHz)${}^{-1/k_{r}}$

kr = 1.84 $\pm$ 0.13 ; 3. 3.

Second segment of 3C 279: 1998.0 – end of observation

b = $-$0.09 $\pm$ 0.01   yr

a = 1.86 $\pm$ 0.28   yr (GHz$)^{-1/k_{r}}$

kr = 1.02 $\pm$ 0.07 . 4. 4.

For BL Lacertae

b = $-$0.13 $\pm$ 0.05   yr

a = 1.52 $\pm$ 0.16   yr (GHz)${}^{-1/k_{r}}$

kr = 1.06 $\pm$ 0.22 ;

From the time lag vs. frequency information (Figs. 8 - 12), the inferred $k_{r}$ are consistent with the equipartition between the magnetic and particle energy densities. The weighted mean $k_{r}$ in all sources $\sim$ 1 except for Segment 1 of 3C 279 where k${}_{r}\sim 1.8$ indicating either a real effect or an anomaly in the data. Large time lags ($>$ 0.5 yr) dominate its 4.8 GHz light curve indicating that opacity effects may not be linear with frequency, resulting in differences compared to expectations from a single shock propagating downstream causing the flare. We assumed $k=1$ in this case to remain consistent with the results for further parameter estimates. From the above fits, we have information on the amplitude and hence the peak flux for each flare at the observed frequencies.

As $r_{core}$ is dependent on $k_{r}$, the inferred values can help in understanding the jet geometry and also provides information on the distribution of magnetic and electron number densities along the jet. As $S_{\nu}\propto\nu^{\alpha}$, the amplitude versus frequency data is fit with a function $A=N~{}\nu^{\alpha}$ to determine the spectral index $\alpha$ for each flare. We find $\alpha$in the range $-$16.45 and 2.41 for S5 0716+714. These are positive in 20 out of 24 cases, and negative for the remaining implying that amplitude increase with frequency is the dominant trend in S5 0716+714 as given in Table 1. The $\alpha$ for 3C 279 Segment 1 are in the range 0.02 to 1.08 and are given in the last column of Table 2 while for Segment 2, $\alpha$ are in the range 0.00 to 1.04 and are given in Column 7 of Table 3. The $\alpha$ for 3C 279 are positive in all cases with a trend of increase in amplitude with the observing frequency. For BL Lacertae, we found $\alpha$ in the range of $-$16.67 and 1.20 as presented in Column 7 of Table 4 and the amplitude was found to increase with observing frequency, with positive $\alpha$ for 22 out of 35 cases while the opposite trend with negative $\alpha$ was found in the remaining cases.

Using the time lags derived from the Gaussian fit procedure, we calculated the core position offset for each flare as it is further used to calculate the deprojected distance of the VLBI core from the base of the jet at a given frequency ($r_{core}$). Using above two quantities, along with Equations 6 to 9, we then obtained magnetic field strength at 1 pc from the base of the jet and also weighted mean value of the magnetic field in the core. All above parameters are given in Table 5 for S5 0716+714, Tables 6 and 7 for Segment 1 and 2 of 3C 279, respectively and in Table 8 for BL Lacertae. The columns of Tables 5$-$8 give the (1) flare nomenclature, (2) observing frequency ($\nu$), (3) core position offset (mas), (4) core position offset in pcGHz, (5) distance of the emitting core from the base of the jet, (6) magnetic field strength at 1 pc distance (in G), (7) equipartition magnetic field strength in mG at the emitting region.

Using the frequency dependent time delays obtainable from the core shift, associated jet parameters can be estimated independent of full VLBI observations. Proper motions of the jet components can be inferred from VLBI observations (e.g. Lister et al. 2013). Pushkarev et al. (2009) measured proper motions for S5 0716+714 jet components in range 0.6 mas yr*-1* to 2.3 mas yr*-1* from which we adopted a median value of 1.4 mas yr*-1*. For 3C 279, they found proper motions in the range 0.005 mas yr*-1* to 0.651 mas yr*-1* between 1995 to 2010 from which we adopted a median value of $\mu=0.35$ mas yr*-1*. Similarly, for BL Lacertae we adopted a median proper motion of 1.08 mas yr*-1*.

In addition, we adopted $\theta=5^{\circ}.2$, $\phi=1^{\circ}.6$, $\delta=10.9$, $z=0.30~{}$ and $D_{L}=1.60$ Mpc for S5 0716+714; $\theta=2^{\circ}.4$, $\phi=0^{\circ}.4$, $\delta=24.1$, $z=0.536~{}$ and $D_{L}=3.08$ Gpc for 3C 279; and $\theta=7^{\circ}.7$, $\phi=1^{\circ}.9$, $\delta=7.2$, $z=0.069~{}$ and $D=3.07$ Gpc for BL Lacertae (Jorstad et al. 2010; Pushkarev et al. 2009; Algaba, Gabuzda, & Smith 2012). The core position offset ($\Omega_{r\nu}$) for S5 0716+714 is $19.81\pm 10.49$ pcGHz; for Segment 1 of 3C 279 is $6.69\pm 2.77$ pcGHz and for Segment 2 is $8.24\pm 2.93$ pcGHz; for BL Lacertae, Mutel et al. (1990) measured an $\Omega_{r\nu}$ of 3.7 pcGHz based on 5.0 and 10.6 GHz observations while O’Sullivan & Gabuzda (2009) measured an $\Omega_{r\nu}$ of 3.4 $\pm$ 0.3 pcGHz. Our estimate of $\Omega_{r\nu}=4.14\pm 2.58$ pcGHz is consistent with these studies.

It is expected for $r_{core}$ to decrease with increasing frequency as evident from Equation (6). For 4 out of 24 flares (G, J, O, R) of S5 0716+714, above trend was inferred with removal of 14.5 GHz point from the flare R while the opposite trend was dominant in 3 out of the remaining (M, Q, S) after removal of 4.8 GHz point from flare M and 8.0 GHz data point from flare S. But no conclusion could be drawn for remaining 17 cases (A, B, C, D, E, F, H, I, K, L, N, P, T, U, V, W, X) as number of data points are at most just two. In case of 3C 279, the expected trend was found in 8 out of 23 (D, E, F, G, L, O, U, V) flares after removing 1 of the data point in some of these such as: 22.2 GHz from flares G, L and U while 14.5 GHz from flare D and F whereas 4.8 GHz from flare E. The opposite trend was also found in 4 out of remaining 15 (B, K, R, W) with removal of 22.2 GHz point from K and W flares while 4.8 GHz from flare R. No conclusion could be drawn for remaining 11 flares, namely A, C, H, I, J, M, N, P, Q, S, T, due to one or two data points. Similarly for BL Lacertae, $r_{core}$ was found to decrease with increasing frequency for 6 out of 36 flares (10, 14, 16, 20, 25, 31) with removal of 22.2 GHz data point from flares 16, 20 and 31, while 14.5 GHz from flare 25. On the other hand, $r_{core}$ increased with frequency in 9 of the remaining cases (4, 7, 13, 15, 17, 23, 24, 27, 32) after ignoring 22.2 GHz point from flare R and 14.5 GHz data point from flares 4, 7, 17, 23, 24, and 27. The remaining cases were not consistent with either of the above patterns and are thus inconclusive owing to lesser number of data points. In the synchrotron opacity model, $r_{\rm core}$ is the opacity $\tau=1$ surface where the core transitions from optically thick to thin due to synchrotron self absorption. This transition is expected to occur at smaller $r_{\rm core}$ for increasing observation frequency, closer to the true core distance in the context of a single shock propagating downstream causing the flaring core. Large time delays ($>$ 0.5 yr) dominating some lower frequency light curves (e.g. 4.8 GHz in 3C 279) result in non-uniformly distributed $r_{\rm core}$ estimates (uncertainty in the location of the $\tau=1$ surface) thus indicating that opacity effects may not be linear with frequency, and implying that the synchrotron opacity model may not be valid in all cases.

For S5 0716+714, 3C 279 (Segments 1 and 2) and BL Lacertae, the inferred $B$ are 0.22 $\pm$ 0.36 G, 0.0005 $\pm$ 0.0004 G, 0.35 $\pm$ 0.45 G and 0.02 $\pm$ 0.06 G respectively. We obtained a weighted mean core magnetic field strength $B_{core}=15\pm 5mG$ for S5 0716+714, for 3C 279 $B_{core}=0.014\pm 0.03mG$ (Segment 1), $B_{core}=17\pm 6mG$ (Segment 2) and $B_{core}=8\pm 2mG$ for BL Lacertae. For 3C 279, our Segment 2 results are consistent with Hirotani (2005) who obtained B = 12 $\pm$ 8 mG.

The 4.8 - 36.8 GHz light curves of S5 0716+714, 3C 279 (segments 1 and 2), and BL Lacertae are analyzed to infer the PSD shape and the results are summarized in Tables 10 - 12 respectively while are plotted in Figures 14 - 17. No statistically significant QPO is detected in any case. For the light curves of S5 0716+714, the power law PSD shape is the better model with a mean slope of $-1.3\pm 0.2$. For the segment 1 light curves of 3C 279, the power law PSD shape is the better model with a mean slope of $-2.2\pm 0.4$. For the segment 2 light curves of 3C 279, the power law PSD model is the better shape for 3/5 light curves (4.8 GHz, 22.2 GHz and 36.8 GHz) with a mean slope of $-2.3\pm 0.4$, indicating that the PSD slope has not changed much between segments 1 and 2. For the remaining segment 2 light curves (8 GHz, 14.5 GHz), either both models could describe the PSD shape or both may not suitably describe them which then necessitates the use of other PSD models such as the generalized bending power law (McHardy et al., 2004), the damped random walk (Kelly, Bechtold, & Siemiginowska, 2009; MacLeod et al., 2010; Sobolewska et al., 2014) or the continuous auto regressive moving average process (Kelly et al., 2014). For the light curves of BL Lacertae, the power law PSD shape is the better model with a mean slope of $-1.8\pm 0.3$. The PSD from AGN light curves are typically power law shaped with higher power density at low frequencies due to long term, higher amplitude variations. If smaller observation frequencies (e.g. 4.8 GHz) are associated with variability mechanisms operational downstream in the jet, the slope is expected to steepen with decreasing observation frequency. We infer similar slopes for a given blazar at all observation frequencies thus indicating no such regular trend. The flares can then arise from multiple propagating shocks originating possibly in distinct events of injecting relativistic electrons into the jet with similar energies. As the inferred $k_{r}\sim 1$ indicates equipartition, a strong magnetic field energy density can also participate in the origin of flares through magnetic instabilities (e.g. Komissarov, 2011) including reconnections or kink instabilities (e.g. Lyutikov, 2003; Uzdensky, 2016) and is supported by simulations which suggest that power law electron distributions can be produced with efficient acceleration (e.g. Sironi et al., 2015).

5.2 Tests of the MAD model

We apply the formalism developed earlier to test the MAD scenario for S5 0716+714, 3C 279 and BL Lacertae. The $B~{}h$ in Gpc values are obtained from the best fit slopes using the data of $B_{\rm core}(r_{\rm core})$ including all the flares for each source.

The data in the log-log space is fit with the model $\log B_{\rm core}=m\log r_{\rm core}+\log c$ resulting in

BL Lacertae: ($m$, $c$) = (-0.094, 8.53 pc mG) 2. 2.

S5 0716+714: ($m$, $c$) = (-0.127, 31.44 pc mG) 3. 3.

3C 279: ($m$, $c$) = (-0.194, 61.86 pc mG)

The resulting $Bh$ at $r_{\rm core}=1$ pc, $w(Bh)$ and $a$ are presented in Table 9. The eqn. (28) predicts that the values of $0.296<w<1$ are disallowed since they yield spin values $a>1$ (see Fig. 13). We find $w$ outside the disallowed region $0.296-1$; the fact that $w$ values are of unity and reasonable values of spin are obtained is interesting and it implies that the magnetic arrest may be operating at some level in these systems. The higher values of spin correlate well with the luminosities which in addition is associated with the strong magnetic fields. Both these facts corroborate the MAD model with a black hole spin assisted relativistic jet (Blandford-Znajek mechanism) possibly operational in this class of objects.

6 Conclusions

A prominent constituent of radio loud AGNs is their core-jet morphology. We studied the core shift effect in the blazars S5 0716+714, 3C 279 and BL Lac using their 4.8 GHz - 36.8 GHz radio light curves obtained from more than three decades of continuous monitoring.

From the time lag and frequency fits, we found the weighted mean $k_{r}$ values as 1.03 for S5 0716+714, 1.84 for segment 1 of 3C 279, 1.02 for segment 2 of 3C 279, while for BL Lacertae it is 1.06 thus being consistent with the equipartition between the magnetic field and particle energy densities. 2. 2.

In case of the BL Lacertae, S5 0716+714, we found $\Omega_{r\nu}$ = $19.81\pm 10.49$ pcGHz for $\mu=1.4$ mas/yr. For $\mu=0.35$ mas/yr in case of 3C 279 we calculated $\Omega_{r\nu}$ = $6.69\pm 2.77$ pcGHz for Segment 1 and $\Omega_{r\nu}$ = $8.24\pm 2.93$ pcGHz for Segment 2. Using $\mu=1.08$ mas/yr in BL Lacertae we found the resultant $\Omega_{r\nu}$ to be $6.29\pm 0.45$ pcGHz. 3. 3.

For S5 0716+714 we found weighted mean $B$ = 0.22 $\pm$ 0.36 G and $B_{core}=15\pm 5mG$. In case of the FSRQ 3C 279, We infer weighted mean $B$ = 0.0005 $\pm$ 0.0004 G and $B_{core}=0.014\pm 0.03mG$ for segment 1 while $B$ = 0.35 $\pm$ 0.45 G and $B_{core}=17\pm 6mG$ for segment 2. For BL Lacertae, we estimated $B$ = 0.02 $\pm$ 0.06 G and $B_{core}=8\pm 2mG$. 4. 4.

The light curves of all three sources were analyzed to infer the PSD shape. No statistically significant QPO was detected in either of the sources. The power law PSD shape is the best fit model for our sources with mean slope of $-1.3$ for S5 0716+714, $-2.3$ for 3C 279 and $-1.8$ for BL Lacertae. 5. 5.

We have estimated the spin from the derived field strength based on reasonable fiducial values and are 0.16 for BL Lacertae, 0.17 for S5 0716+714 and 0.9 for 3C 279. Independent estimates of spin from the Iron line or other techniques if available, can greatly help in reducing uncertainties. Using more accurate measurements of the field strength, bolometric and the jet luminosity, one can explore other electrodynamical and hybrid jet models and place constraints on the spin and mass of the black hole. 6. 6.

The methodology can be applied to a large sample of quasars to estimate black hole spin, compare with estimates from independent methods, study the distribution of inferred spins, and relate this to source properties including radio loudness, jet luminosity, redshift, amongst others.

These studies give the evidence of frequency dependent core shifts and also encourage us further to calculate jet parameters in the region unresolved by VLBI technique. The linear polarization based electric vector position angle (EVPA) follows a bimodal distribution with a tendency of polarization orthogonal to the jet in quasars and along the jet in BL Lac objects (e.g. Cawthorne et al., 1993; Gabuzda & Cawthorne, 2003). The shock mechanism produces magnetic fields compressed in the plane of the shock, resulting in polarization along the jet for a transverse shock. As shocks are transient events, the resultant jet polarization orientation may not be retained over long durations. Further, as internal shocks can be oblique, with a range of orientations, a natural bimodal distribution of the relative EVPAs cannot be expected. The observation of a bimodal distribution has been inferred in the computation of the linear polarization and the EVPA variability in the context of a relativistic jet flow carrying large-scale helical magnetic fields (e.g. Lyutikov et al., 2005). The measurement of the field strength and its variation along the length of the pc-scale jet using the core shift method can aid detailed polarization models. A particular model (Mangalam, 2017; in preparation) invokes the assumption of helical magnetic structure and special relativistic kinematics of a blob transiting the jet to compute the instantaneous Doppler factor and hence predict the varying flux density, linear polarization and EVPA. Radio continuum surveys with the SKA can provide AGN imaging and timing (flux density, polarization) data across a range of frequencies simultaneously, ideally suited for the application of the time delay method. The SKA-VLBI synergy can provide key and unique contributions to the current study and additionally, in the context of studies of magnetic fields in AGN jets (Roy et al., 2016). Some of these include (1) verification of the SSA model using high resolution sub–pc-scale images, as the SSA opacity is related to the magnetic field strength and core distance (e.g. Lobanov, 1998), (2) offering a sizeable sample for statistical and comparative studies of pc-scale jet kinematics and energetics which can serve as crucial inputs to fully relativistic radiation magneto-hydrodynamic simulations, and (3) exploiting the higher sensitivity to understand evolution of magnetic fields in these AGNs over cosmic time.

Acknowledgments

We thank the referee for a constructive discussion and suggestions which have improved the manuscript. PM is supported by the Chinese Academy of Sciences President’s International Fellowship Initiative (CAS-PIFI; grant no. 2016PM024) post-doctoral fellowship and the National Natural Science Foundation of China (NSFC) Research Fund for International Young Scientists (grant no. 11650110438). ACG is partially supported by CAS President’s International Fellowship Initiative (PIFI; grant no. 2016VMB073). The work at UMRAO was supported in part by a series of grants from the NSF and NASA, most recently AST-0607523 and NASA Fermi GI grants NNX09AU16G, NNX10AP16G, NNX11AO13G, and NNX13AP18G. MFG is supported by the National Science Foundation of China (grants 11473054 and U1531245) and by the Science and Technology Commission of Shanghai Municipality (grant 14ZR1447100).