1 Introduction
Many extreme astrophysical environments, especially those associated with black holes (active galactic nuclei, especially the relativistically beamed blazars; stellar X-ray binaries, especially microquasars) and neutron stars (magnetars, pulsars and their associated nebulae), as well as gamma-ray bursts, are efficient emitters of gamma-ray radiation. Their broadband, power-law radiation spectra clearly indicate non-thermal, collisionless plasma environments that enable efficient particle acceleration to ultra-relativistic energies. In addition, many of these sources show powerful gamma-ray flaring activities with short variability time scales and high radiation efficiencies (e.g. Buehler et al. Reference Buehler, Scargle, Blandford, Baldini, Baring, Belfiore, Charles, Chiang, D’Ammando and Dermer2012; Ackermann et al. Reference Ackermann, Anantua, Asano, Baldini, Barbiellini, Bastieri, Becerra Gonzalez, Bellazzini, Bissaldi and Blandford2016), suggesting that the particle acceleration process proceeds rapidly and is likely a consequence of electromagnetic dissipation in a highly magnetised outflow from the central engine (e.g. Blandford et al. Reference Blandford, Yuan, Hoshino and Sironi2017).
One of the most promising dissipation and particle acceleration mechanisms is relativistic magnetic reconnection. Reconnection becomes relativistic when the magnetic energy density dominates the plasma rest-mass energy density, which may well be the case in the aforementioned astroplasma environments. Recently, significant progress has been made in understanding the basics of relativistic reconnection, with increasingly sophisticated fully kinetic numerical simulations employing the particle-in-cell (PIC) algorithm (e.g. Sironi & Spitkovsky Reference Sironi and Spitkovsky2014; Guo et al. Reference Guo, Li, Daughton and Liu2014; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016; Kagan, Nakar & Piran Reference Kagan, Nakar and Piran2016; Sironi, Giannios & Petropoulou Reference Sironi, Giannios and Petropoulou2016).
Kinetic numerical simulations of relativistic magnetic reconnection are mostly local, starting from an initial configuration involving highly symmetric predefined Harris-type current layers. In the context of global astrophysical systems like pulsar magnetospheres or relativistic jets of active galaxies, reconnection can be expected to be an intermittent process induced by dynamical collisions of misaligned magnetic domains or global current-driven instabilities (Mizuno et al. Reference Mizuno, Lyubarsky, Nishikawa and Hardee2009; Lyutikov et al. Reference Lyutikov, Komissarov, Sironi and Porth2018). Realistic (three-dimensional) global kinetic simulations of pulsar magnetospheres (Philippov, Spitkovsky & Cerutti Reference Philippov, Spitkovsky and Cerutti2015) or relativistic jets (Nishikawa Reference Nishikawa, Frederiksen, Nordlund, Mizuno, Hardee, Niemiec, Gómez, Pe’er, Duţan and Sol2016) are computationally very challenging. However, intermittent reconnection with dynamical formation of kinetically thin current layers can now be investigated with a novel local magnetostatic configuration referred to as ‘ABC fields’ (see § 2).
It has also been possible to investigate the effects of radiation reaction to obtain self-consistent radiative signatures from particles accelerated by relativistic magnetic reconnection. The primary radiative mechanism investigated so far is synchrotron radiation (SYN) with direct application to the problem of rapid gamma-ray flares from the Crab Nebula (e.g. Cerutti et al. Reference Cerutti, Werner, Uzdensky and Begelman2013, Reference Cerutti, Werner, Uzdensky and Begelman2014), and synchro-curvature radiation with application to the gamma-ray emission of pulsars (e.g. Cerutti, Philippov & Spitkovsky Reference Cerutti, Philippov and Spitkovsky2016).
In other cases, especially in blazars and microquasars, inverse Compton radiation can be equally important, or even dominate over the synchrotron mechanism. Since blazars explain the bulk of cosmic gamma-ray radiation in the GeV and TeV bands (Ajello et al.
Reference Ajello, Gasparrini, Sánchez-Conde, Zaharijas, Gustafsson, Cohen-Tanugi, Dermer, Inoue, Hartmann and Ackermann2015), this radiation is produced mainly in the inverse Compton (IC) process. The simplest radiative models of blazars assume that the observed synchrotron and IC spectral components are produced by the same population of energetic electrons. For a given blazar, the ratio of IC to synchrotron luminosities, known as the Compton dominance
$\text{CD}=L_{\text{IC}}/L_{\text{syn}}$
, is an indirect probe of plasma magnetisation in relativistic jets. If the IC scattering proceeds in the Thomson regime, one can find that
$\text{CD}\sim U_{\text{rad}}^{\prime }/U_{\text{B}}^{\prime }$
, where
$U_{\text{rad}}^{\prime }$
is the energy density of soft radiation fields (external or synchrotron) and
$U_{\text{B}}^{\prime }$
is the magnetic energy density, both evaluated in the jet co-moving frame (Janiak, Sikora & Moderski Reference Janiak, Sikora and Moderski2015).
In this work, we present the first results from two- and three-dimensional PIC simulations of relativistic magnetic reconnection initiated from the ‘ABC fields’, taking into account radiation reaction from both synchrotron and IC processes.Footnote 1 The IC radiation is calculated as upscattering of fixed monoenergetic uniform isotropic radiation field, which corresponds to the external radiation Compton (ERC) mechanism (see the Appendix) operating in relativistic jets of the most luminous blazars. We investigate the spectral and temporal distributions of simulated radiative signatures of particle acceleration in dynamical current layers that form between highly magnetised coalescing magnetic domains. We also study the value of Compton dominance in relation to the adopted initial ratio of radiative and magnetic energy densities, as well as the perpendicular profiles of the current layers.
Details of our numerical configuration are provided in § 2. Section 3 presents the results on the total energy transformations and radiative efficiency. Section 4 describes the spectral distribution of simulated synchrotron and IC emission. Section 5 presents the results on time variability of simulated radiation signals. In § 6 we investigate the effect of strong radiative cooling on the perpendicular structure of current layers. Our results are discussed in § 7. In the Appendix, we provide a short summary of the observational properties of blazars and the characteristic physical parameters of their relativistic jets.
2 Numerical set-up
We perform particle-in-cell numerical simulations of periodic plasma configurations referred to as ABC fields (from Arnold–Beltrami–Childress; see Dombre et al. Reference Dombre, Frisch, Henon, Greene and Soward1986; East et al. Reference East, Zrake, Yuan and Blandford2015; Lyutikov et al. Reference Lyutikov, Sironi, Komissarov and Porth2017) using a modified version of the Zeltron code (Cerutti et al. Reference Cerutti, Werner, Uzdensky and Begelman2013). In three dimensions, the magnetic fields are defined as:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}B_{x}(x,y,z)=B_{0}[\sin (\unicode[STIX]{x1D6FC}z)+\cos (\unicode[STIX]{x1D6FC}y)],\\ B_{y}(x,y,z)=B_{0}[\sin (\unicode[STIX]{x1D6FC}x)+\cos (\unicode[STIX]{x1D6FC}z)],\\ B_{z}(x,y,z)=B_{0}[\sin (\unicode[STIX]{x1D6FC}y)+\cos (\unicode[STIX]{x1D6FC}x)],\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$B_{0}$
is the nominal magnetic field strength and
$\unicode[STIX]{x1D6FC}=2\unicode[STIX]{x03C0}k/L$
with
$k$
the wavenumber and
$L$
the linear size of the simulation domain:
$x,y,z\in [0:L]$
. The case of
$k=1$
corresponds to the lowest-energy stable Taylor state. The case of
$k=2$
is the lowest-energy unstable configuration. In two dimensions there exists a smaller unstable configuration obtained by rotating the coordinate system (Nalewajko et al.
Reference Nalewajko, Zrake, Yuan, East and Blandford2016; Yuan et al.
Reference Yuan, Nalewajko, Zrake, East and Blandford2016).
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}B_{x}(x,y)=B_{0}[\sin (\unicode[STIX]{x1D6FC}(x+y))+\sin (\unicode[STIX]{x1D6FC}(x-y))]/\sqrt{2},\\ B_{y}(x,y)=B_{0}[\sin (\unicode[STIX]{x1D6FC}(x-y))-\sin (\unicode[STIX]{x1D6FC}(x+y))]/\sqrt{2},\\ B_{z}(x,y)=B_{0}[\cos (\unicode[STIX]{x1D6FC}(x+y))-\cos (\unicode[STIX]{x1D6FC}(x-y))].\end{array}\right\} & & \displaystyle\end{eqnarray}$$
We will refer to this configuration as ‘diagonal ABC fields’. We note that this initial magnetic field distribution is characterised by
$\max (B)=2B_{0}$
and
$\langle B^{2}\rangle =2B_{0}^{2}$
.
In order to obtain a rough equilibrium, current density
$\boldsymbol{j}=(kc/L)\boldsymbol{B}$
, with
$c$
the speed of light, is provided by a population of relativistic electrons and positrons characterised by uniform total number density
$n=3\sqrt{2}kB_{0}/(2e\tilde{a}_{1}L)$
, with
$e$
the elementary charge and
$\tilde{a}_{1}=B_{0}/\max (B)=0.5$
a constant; by the Maxwell–Jüttner energy distribution with relativistic temperature
$\unicode[STIX]{x1D6E9}_{e}=k_{\text{B}}T_{e}/(m_{e}c^{2})$
, with
$k_{\text{B}}$
the Boltzmann constant,
$T_{e}$
the temperature and
$m_{e}$
the electron mass; and by the dipole moment of the local angular distribution of particle momenta
$a_{1}=(B/B_{0})\tilde{a}_{1}$
. As we work here in the limit of highly relativistic electron temperatures
$\unicode[STIX]{x1D6E9}_{e}\gg 1$
, we use the following statistics of the Maxwell–Jüttner distribution
$\langle \unicode[STIX]{x1D6FE}\rangle \simeq 3\unicode[STIX]{x1D6E9}_{e}$
and
$\langle \unicode[STIX]{x1D6FE}^{2}\rangle \simeq 12\unicode[STIX]{x1D6E9}_{e}^{2}$
, with
$\unicode[STIX]{x1D6FE}$
the particle Lorentz factor. The corresponding mean hot magnetisation value is given by
$\langle \unicode[STIX]{x1D70E}_{\text{hot}}\rangle =\langle B^{2}\rangle /(4\unicode[STIX]{x03C0}w)$
, where
$w\simeq 4\unicode[STIX]{x1D6E9}_{e}nm_{e}c^{2}$
is the ultra-relativistic specific enthalpy:
$$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D70E}_{\text{hot}}\rangle \simeq \frac{\tilde{a}_{1}}{12\sqrt{2}\unicode[STIX]{x03C0}k}\left(\frac{L}{\unicode[STIX]{x1D70C}_{0}}\right), & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}_{0}=\unicode[STIX]{x1D6E9}_{e}m_{e}c^{2}/(eB_{0})$
is the nominal particle gyroradius. The characteristic property of ABC fields is that magnetisation scales linearly with the scale separation between the magnetic field coherence scale (of order
$L$
) and the kinetic gyration scale (
$\unicode[STIX]{x1D70C}_{0}$
) (Nalewajko et al.
Reference Nalewajko, Zrake, Yuan, East and Blandford2016). This is because a minimum particle number density is required in order to realise a smoothly distributed current density.
Radiation reaction from synchrotron and inverse Compton processes was included in advancing the particle momenta. For the synchrotron process, local values of magnetic and electric fields
$\boldsymbol{B},\boldsymbol{E}$
are used for an electron with velocity
$\boldsymbol{v}=\unicode[STIX]{x1D737}c$
, Lorentz factor
$\unicode[STIX]{x1D6FE}=(1-\unicode[STIX]{x1D6FD}^{2})^{-1/2}$
and dimensionless 4-velocity
$\boldsymbol{u}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D737}$
(Cerutti et al.
Reference Cerutti, Werner, Uzdensky and Begelman2013):
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t} & = & \displaystyle -\frac{P_{\text{syn}}\boldsymbol{u}}{\unicode[STIX]{x1D6FE}m_{e}c^{2}},\quad P_{\text{syn}}=\frac{\unicode[STIX]{x1D70E}_{\text{T}}c}{4\unicode[STIX]{x03C0}}[(\unicode[STIX]{x1D6FE}\boldsymbol{E}+\boldsymbol{u}\times \boldsymbol{B})^{2}-(\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{u})^{2}].\end{eqnarray}$$
For the inverse Compton process, we assume a uniform isotropic external radiation field parametrised by energy density
$U_{\text{ext}}$
and photon energy
$E_{\text{ext}}$
. In order to account for the Klein–Nishina cross-section in the Zeltron code, we have updated the radiation reaction formula of Cerutti et al. (Reference Cerutti, Werner, Uzdensky and Begelman2013) (Jones Reference Jones1968):
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\frac{P_{\text{IC}}\boldsymbol{u}}{\unicode[STIX]{x1D6FE}m_{e}c^{2}},\quad P_{\text{IC}}=\frac{4}{3}\unicode[STIX]{x1D70E}_{T}cU_{\text{ext}}u^{2}f_{\text{KN}}(b),\quad b=\frac{4\unicode[STIX]{x1D6FE}E_{\text{ext}}}{m_{e}c^{2}}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle f_{\text{KN}}(b) & = & \displaystyle \frac{9}{b^{3}}\left[\left(\frac{b}{2}+6+\frac{6}{b}\right)\ln (1+b)-\frac{1}{(1+b)^{2}}\left(\frac{11}{12}b^{3}+6b^{2}+9b+4\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,2+2\text{Li}_{2}(-b)\right],\end{eqnarray}$$
where
$\unicode[STIX]{x1D70E}_{T}$
is the Thomson cross section.
The synchrotron radiation spectra are calculated by summing spectral contributions from all individual macroparticles (electrons and positrons) of velocity
$\boldsymbol{v}$
and Lorentz factor
$\unicode[STIX]{x1D6FE}$
(Blumenthal & Gould Reference Blumenthal and Gould1970):
$$\begin{eqnarray}\displaystyle & \displaystyle L_{\text{syn}}(\unicode[STIX]{x1D708})=\frac{\sqrt{3}e^{2}}{c}\mathop{\sum }_{e^{+}e^{-}}N_{e,1}F(\unicode[STIX]{x1D709})\unicode[STIX]{x1D6FA}_{\text{syn}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D709}\int _{\unicode[STIX]{x1D709}}^{\infty }K_{5/3}(x)\,\text{d}x, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D709}=\frac{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D708}}{3\unicode[STIX]{x1D6FE}^{2}\unicode[STIX]{x1D6FA}_{\text{syn}}},\quad \unicode[STIX]{x1D6FA}_{\text{syn}}=\frac{e}{m_{e}c}|(\boldsymbol{E}+\boldsymbol{n}\times \boldsymbol{B})\times \boldsymbol{n}|, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{n}=\boldsymbol{v}/|\boldsymbol{v}|$
is the unit vector along the particle velocity,
$N_{e,1}$
is the number of electrons represented by individual macroparticle and
$K_{a}$
is the modified Bessel function of the second kind. Because we consider ultra-relativistic particles with
$\unicode[STIX]{x1D6E9}_{e}\gg 1$
, our approach is self-consistent for
$\unicode[STIX]{x1D709}=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}>(\unicode[STIX]{x1D714}_{c}\unicode[STIX]{x0394}t)^{-1}=(5\unicode[STIX]{x1D6E9}_{e}^{3})^{-1}$
, where
$\unicode[STIX]{x1D714}_{c}=(3/2)\unicode[STIX]{x1D6FE}^{2}\unicode[STIX]{x1D6FA}_{\text{syn}}\simeq 18\unicode[STIX]{x1D6E9}_{e}^{3}(c/\unicode[STIX]{x1D70C}_{0})$
and
$\unicode[STIX]{x0394}t=0.99\unicode[STIX]{x0394}x/(\sqrt{2}c)$
is the simulation time step, and hence there is no need for a more general Fourier transformation method of Hededal & Nordlund (Reference Hededal and Nordlund2005).
The IC radiation spectra are calculated in an analogous way, using the general Klein–Nishina kernel (Blumenthal & Gould Reference Blumenthal and Gould1970):
$$\begin{eqnarray}\displaystyle & \displaystyle L_{\text{IC}}(\unicode[STIX]{x1D708})=\frac{12h\unicode[STIX]{x1D70E}_{\text{T}}U_{\text{ext}}}{m_{e}c}\mathop{\sum }_{e^{+}e^{-}}N_{e,1}\frac{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D709}}{b^{2}}K(\unicode[STIX]{x1D709},b), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle K(\unicode[STIX]{x1D709},b)=1+q-2q^{2}+2q\log q+q(1-q)\frac{b\unicode[STIX]{x1D709}}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle q=\frac{\unicode[STIX]{x1D709}}{b(1-\unicode[STIX]{x1D709})}<1,\quad \unicode[STIX]{x1D709}=\frac{h\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D6FE}m_{e}c^{2}}<1,\quad b=\frac{4\unicode[STIX]{x1D6FE}E_{\text{ext}}}{m_{e}c^{2}}. & & \displaystyle\end{eqnarray}$$
We performed two series of two-dimensional (2-D) simulations initiated from diagonal ABC fields with
$k=1$
on numerical grid of size
$N_{x}=N_{y}=2048$
at resolution
$\unicode[STIX]{x0394}x^{i}=\unicode[STIX]{x1D70C}_{0}/2.56$
, and with 128 particles per cell. In all cases, the energy density of external radiation fields was set at
$U_{\text{ext}}=\langle B^{2}\rangle /8\unicode[STIX]{x03C0}=B_{0}^{2}/4\unicode[STIX]{x03C0}=2U_{0}$
, where
$U_{0}=B_{0}^{2}/8\unicode[STIX]{x03C0}$
is the nominal magnetic energy density. The domain size of
$L=800\unicode[STIX]{x1D70C}_{0}$
corresponds to the mean hot magnetisation value of
$\langle \unicode[STIX]{x1D70E}_{\text{hot}}\rangle \simeq 7.5$
. Parameters of all PIC simulations presented in this work are compared in table 1.
The first series, denoted as ‘2Da’, was performed for the same value of
$B_{0}=1~\text{G}$
and for different values of initial particle temperature:
$\unicode[STIX]{x1D6E9}_{e}=10^{5},3\times 10^{5},10^{6},3\times 10^{6},10^{7}$
. Increasing the particle temperature leads to increasing the efficiency of radiative cooling, and hence a shorter cooling length. The nominal cooling length with IC cooling in the Thomson limit is given by:
$$\begin{eqnarray}\displaystyle l_{\text{cool}}=c\unicode[STIX]{x1D70F}_{\text{cool}}=\frac{\langle \unicode[STIX]{x1D6FE}\rangle }{\langle |\text{d}\unicode[STIX]{x1D6FE}/c\text{d}t|\rangle }=\frac{\langle \unicode[STIX]{x1D6FE}\rangle }{\langle \unicode[STIX]{x1D6FE}^{2}\rangle }\frac{3m_{e}c^{2}}{4\unicode[STIX]{x1D70E}_{T}U_{\text{cool}}}\simeq \frac{(3\unicode[STIX]{x03C0}/8)e}{\unicode[STIX]{x1D70E}_{T}\unicode[STIX]{x1D6E9}_{e}^{2}B_{0}}\unicode[STIX]{x1D70C}_{0}, & & \displaystyle\end{eqnarray}$$
where
$U_{\text{cool}}=\langle U_{B}\rangle +U_{\text{ext}}=4U_{0}$
is the effective cooling energy density. The values of
$l_{\text{cool}}/\unicode[STIX]{x1D70C}_{0}$
are reported in table 1.
The second series, denoted as ‘2Db’, was performed for the same value of initial particle temperature:
$\unicode[STIX]{x1D6E9}_{e}=10^{4}$
, and for different values of
$B_{0}$
and
$E_{\text{ext}}$
. Two of these simulations were performed in the fast-cooling (fc) regime by setting
$B_{0}=2\times 10^{4}~\text{G}$
, and two simulations were performed in the slow-cooling regime (sc) by setting
$B_{0}=2\times 10^{2}~\text{G}$
. We note that such high magnetic field strengths have been suggested in the context of the most extreme GeV gamma-ray variability of blazar 3C 279 observed by Fermi/LAT on suborbital time scales of a few minutes (Ackermann et al.
Reference Ackermann, Anantua, Asano, Baldini, Barbiellini, Bastieri, Becerra Gonzalez, Bellazzini, Bissaldi and Blandford2016). Furthermore, two simulations were performed with Compton scattering in the Klein–Nishina (kn) regime by setting
$E_{\text{ext}}=40~\text{eV}$
, and two simulations were performed with Compton scattering in the Thomson (th) regime by setting
$E_{\text{ext}}=0.4~\text{eV}$
.
Components of the total energy as functions of simulation time, normalised to the total energy including the total radiation losses. (a–d) The ‘2Da’ series; (e–h) the ‘2Db’ series; (i–l) the 3-D simulation. In (c,g,k), the solid lines show the contribution of positrons, and the dashed lines show the contribution of electrons. In (d,h,l), the solid lines show integrated synchrotron energy losses, and the dashed lines show the integrated IC energy losses.
List of PIC simulations described in this work.
$L/\unicode[STIX]{x1D70C}_{0}$
is the physical size of the simulation domain in units of nominal gyroradius
$\unicode[STIX]{x1D70C}_{0}=\unicode[STIX]{x1D6E9}_{e}m_{e}c^{2}/(eB_{0})$
,
$B_{0}$
is the nominal magnetic field strength in units of Gauss,
$\unicode[STIX]{x1D6E9}_{e}=k_{\text{B}}T_{e}/(m_{e}c^{2})$
is the dimensionless relativistic temperature of the initial particle energy distribution,
$E_{\text{ext}}$
is the photon energy of soft radiation in units of eV,
$l_{\text{cool}}$
is the nominal cooling length defined in (2.13) and
$\unicode[STIX]{x1D70F}_{\text{E}}$
is the exponential growth rate of total electric energy in units of light-crossing time scale
$L/c$
.
All of our 2-D simulations show the same basic behaviour as that reported in earlier works (Nalewajko et al.
Reference Nalewajko, Zrake, Yuan, East and Blandford2016; Yuan et al.
Reference Yuan, Nalewajko, Zrake, East and Blandford2016). Initial equilibrium evolves due to coalescence instability, which leads to the formation of dynamical current layers over a single light-crossing time scale
$L/c$
. This is followed by slowly damped nonlinear global oscillations. Particles are accelerated most efficiently by electric fields present in linear current layers.
We are also reporting partial results from a single 3-D simulation for ‘standard’ ABC fields with
$k=2$
and
$\tilde{a}_{1}=0.4$
, that was performed with the same prescription for radiative losses. This simulation was performed on numerical grid of size
$N_{x}=N_{y}=N_{z}=1152$
at resolution
$\unicode[STIX]{x0394}x^{i}=\unicode[STIX]{x1D70C}_{0}/1.28$
with 16 particles per cell. Detailed results of this simulation will be presented in another publication.
3 Global energy transformations and radiative efficiency
Figure 1 compares the time evolutions of total energy components, which include the magnetic, electric, kinetic energies, as well as the global energy radiated away in the synchrotron and inverse Compton processes.
3.1 Results of the 2-D simulations
In all our 2-D simulations, the total energy of the system is dominated by magnetic energy, with the initial share of 83 %, decreasing to
$\simeq 68\,\%$
by the moment of saturation of the linear instability (
$ct/L\simeq 1.9$
), followed by a further slow decrease to
$\simeq 63\,\%$
by
$ct/L\simeq 7$
. Evolution of the magnetic energy fraction is very similar in all cases. The electric energy reaches a peak of
$\simeq 10\,\%$
at the saturation point, with slightly higher values in the fast-cooling regimes. Major differences are seen in the evolutions of energy shared by the particles. In the ‘2Da’ series, the initial fraction of
$\simeq 16\,\%$
(including both electrons and positrons) increases to
${\sim}30\,\%$
for
$\unicode[STIX]{x1D6E9}_{e}=10^{5}$
, but it decreases down to
${<}1\,\%$
for
$\unicode[STIX]{x1D6E9}_{e}=10^{7}$
. In the ‘2Db’ series, the kinetic energy fraction reaches between
$\simeq (30{-}32)\,\%$
in the slow-cooling cases, and
${<}10\,\%$
in the fast-cooling cases. These differences are reflected in the radiative efficiencies, which account for cumulative radiative energy losses of all particles. In the slow-cooling cases, the radiative efficiencies of our simulations are of order
${\sim}(3{-}5)\,\%$
, counting both synchrotron and IC channels.
The most interesting result here is the radiative efficiency in the fast-cooling regime. In the ‘2Da’ series, by
$ct/L=3$
, the synchrotron radiative efficiency reaches 10 % for
$\unicode[STIX]{x1D6E9}_{e}=10^{7}$
, and the IC radiative efficiency becomes even higher, up to
$\simeq 17\,\%$
. We also find that the IC radiative efficiency dominates the synchrotron radiative efficiency for
$\unicode[STIX]{x1D6E9}_{e}\geqslant 10^{6}$
, while they are comparable to each other for
$\unicode[STIX]{x1D6E9}_{e}<10^{6}$
. In the ‘2Db’ series, by
$ct/L=7$
, with IC cooling in the Klein–Nishina regime (fc_kn), we achieved a synchrotron radiative efficiency of
$\simeq 23\,\%$
, while the IC radiative efficiency was only
$\simeq 2\,\%$
. However, in the Thomson regime (fc_th), the IC radiative efficiency of
$\simeq 23\,\%$
dominates the synchrotron efficiency of
$\simeq 8\,\%$
, and this is despite these efficiencies being initially (
$ct/L<1$
) comparable due to our choice of
$U_{\text{ext}}/U_{\text{B0}}$
. While
$U_{\text{ext}}$
is constant, the magnetic energy density
$U_{\text{B}}$
decreases during the simulation by
$\simeq 25\,\%$
. On the other hand, the ratio of cumulative radiative efficiencies is
$f_{\text{IC}}/f_{\text{syn}}\sim 3$
(we also note that this ratio is less than unity in the sc_th case), which means that the local magnetic field component perpendicular to the velocities of emitting particles is systematically reduced. There can be two reasons for this: (i) that particles are radiating preferentially in regions of low magnetic field strength, (ii) that particles are radiating preferentially along local magnetic field lines. Both possibilities are consistent with emission being produced primarily within the reconnection current layers and their outflowing nozzles (Yuan et al.
Reference Yuan, Nalewajko, Zrake, East and Blandford2016).
We also measure growth time scales of the linear coalescence instability defined as the
$e$
-folding time scale
$\unicode[STIX]{x1D70F}$
of the total electric energy
$E^{2}$
, normalised to the light-crossing time scale
$L/c$
. In the ‘2Da’ series, we find systematic decrease of growth time scale with increasing cooling efficiency, from
$c\unicode[STIX]{x1D70F}/L\simeq 0.18$
for
$\unicode[STIX]{x1D6E9}_{e}\leqslant 3\times 10^{5}$
to
$c\unicode[STIX]{x1D70F}/L\simeq 0.145$
for
$\unicode[STIX]{x1D6E9}_{e}=10^{7}$
. In the ‘2Db’ series, the longest growth time scale
$c\unicode[STIX]{x1D70F}/L\simeq 0.165$
is measured in the ‘fc_th’ case, and the shortest growth time scales of
$c\unicode[STIX]{x1D70F}/L\simeq 0.18$
are measured in the ‘sc’ cases. We thus find a consistent picture that the instability growth time scale is enhanced by strong radiative cooling, which also results in removing the gas pressure and increasing effective magnetisation
$\unicode[STIX]{x1D70E}$
. This is consistent with the previous result showing that the growth time scale decreases systematically with
$\unicode[STIX]{x1D70E}$
and is longer than the force-free limit of
$c\unicode[STIX]{x1D70F}_{\text{ff}}/L\simeq 0.13$
(Nalewajko et al.
Reference Nalewajko, Zrake, Yuan, East and Blandford2016).
3.2 Results of the 3-D simulation
Figure 1 shows that the initial state is characterised by lower magnetisation, with magnetic fields containing
${\sim}71\,\%$
of the total energy. Correspondingly, more energy (
${\sim}29\,\%$
; summing electrons and positrons) is contained in the particles. The overall cooling rate is very fast, so that before magnetic reconnection is able to compete with the radiative cooling, the kinetic energy content of particles decreases to
${\sim}16\,\%$
. The exponential growth time scale of electric energy is
$c\unicode[STIX]{x1D70F}/L\simeq 0.19$
, but it saturates at the level of
${\sim}5\,\%$
, just half of the 2-D level. However, noting that our 3-D configuration is characterised by an ABC wavelength shorter by factor
$\sqrt{2}$
compared with the 2-D configuration, the effective growth time scale of the 3-D simulation would amount to 0.27. By
$ct/L=6.5$
, the radiative efficiency of the synchrotron mechanism is
$\simeq 24\,\%$
, dominating the radiative efficiency of the IC mechanism of
$\simeq 12\,\%$
. The Klein–Nishina parameter
$b\simeq 0.25$
is somewhat higher than in the case 2Db_fc_th (
$b\simeq 0.1$
). Although still technically in the Thomson regime, the difference in
$b$
values can partly explain the difference between synchrotron-dominated radiative losses in three dimensions and the IC-dominated losses in two dimensions.
4 Spectral energy distribution of radiation
Figure 2 shows the time evolution of the
$\unicode[STIX]{x1D708}F_{\unicode[STIX]{x1D708}}$
radiation spectral energy distributions (SED), showing separately the synchrotron and IC components. All panels are scaled in the same way with respect to the initial synchrotron luminosity
$F_{\text{SYN},0}$
, and to the respective peak frequencies. Corresponding to these SEDs, figure 3 compares the positions
$(\unicode[STIX]{x1D708}_{p},F_{p})$
of respective SED peaks.
Time evolution of radiation spectra – synchrotron (a) and IC (b) – compared for simulations in the ‘2Db’ series and for the 3-D simulation. The line colours indicate the simulation time progressing from deep blue (
$ct/L\leqslant 1$
) to brown (
$ct/L\geqslant 4$
). The spectra are normalised to the initial peak of the synchrotron component. The grey stripes indicate the frequency ranges from which lightcurves presented in figure 4 are extracted.
4.1 Results of the 2Db simulations
In the slow-cooling regime, we can see how a high-frequency synchrotron component develops due to energetic particles emerging at the sites of magnetic reconnection. This component has peak frequency of
$\unicode[STIX]{x1D708}_{\text{syn},p}\sim 200\unicode[STIX]{x1D708}_{\text{syn},0}$
, and peak luminosity of
$F_{\text{syn},p}\sim 3F_{\text{syn},0}$
. This component does not depend on whether the IC cooling proceeds in the Thomson or Klein–Nishina regime. However, a major difference can be seen in the IC component: as can be expected, it is significantly suppressed in the ‘sc_kn’ case, and the SED remains dominated by low-energy electrons. In the ‘sc_th’ case, a high-frequency peak is obtained at
$\unicode[STIX]{x1D708}_{\text{IC},p}\sim 100\unicode[STIX]{x1D708}_{\text{IC},0}$
and
$F_{\text{IC},p}\sim F_{\text{IC},0}$
. The Compton dominance in the ‘sc_th’ case is in the range
$\text{CD}\in [0.5:1]$
, while it decreases to
${\sim}0.02$
in the ‘kn’ case. The peak frequency ratio is of order
$\unicode[STIX]{x1D708}_{\text{IC},p}/\unicode[STIX]{x1D708}_{\text{syn},p}\sim 0.01(\unicode[STIX]{x1D708}_{\text{IC},0}/\unicode[STIX]{x1D708}_{\text{syn},0})$
.
The situation is different in the fast-cooling regime, where strong radiative losses force a rapid evolution of the initial SEDs. Initially, the synchrotron luminosity decreases more strongly than the IC luminosity, and hence the Compton dominance increases. In the ‘fc_th’ case, this increase is by one order of magnitude, and we have
$\text{CD}\sim 10$
through the end of simulation. In the ‘fc_kn’ case, the increase is by factor 4, but it is reversed due to rebrightening of the synchrotron component. The latter case is also characterised by minor evolution of the IC component.
4.2 Results of the 3-D simulation
Figure 2 shows fairly similar time evolution of the synchrotron and IC SED components. Figure 3 demonstrates similar initial behaviours of the synchrotron peaks between three dimensions and 2Db_fc_th, and a notably stable location of the inverse Compton/synchrotron (IC/SYN) peak ratio with Compton dominance of the order of unity.
Time evolution of radiation SED peaks – synchrotron (a), IC (b) and Compton dominance (IC/SYN; c) – compared for simulations in the ‘2Db’ series and for the 3-D simulation.
5 Lightcurves
We calculate lightcurves of synchrotron and inverse Compton radiation for a specific observer and for several frequency bands, for all simulations in the ‘2Db’ series and for the 3-D one. Figure 2 shows the frequency bands used for calculating the lightcurves presented in figure 4. We focus our attention on the spectral bands dominated by contribution from particles accelerated by magnetic reconnection.
In figure 4, we compare directly synchrotron and IC lightcurves normalised independently to unity, as would be observed simultaneously by the same observer. Major difference can be seen between the fast- and slow-cooling regimes. In the former (including the 3-D case), a brief flash lasting a little over a single light-crossing time
$\unicode[STIX]{x0394}(ct/L)\sim 1$
would be seen. In the ‘fc_th’ case, we also see a strong signal of synchrotron radiation for
$ct/L<1$
produced by the initial population of very hot particles. In the slow-cooling regime, the lightcurves are more extended in time, and they are expected by decay at a very slow rate. A quasi-periodic pattern at the period of
${\sim}L/(2c)$
is apparent in both synchrotron and IC lightcurves. No major difference is seen between the Thomson and Klein–Nishina regimes of the IC process.
(a) Lightcurves calculated for simulations in the ‘2Db’ series and for the 3-D one in the frequency ranges indicated in figure 2. These lightcurves are smoothed and normalised independently to their respective maximum values. (b) Autocorrelation functions (ACF) and discrete correlation functions (DCF) calculated for the sections of lightcurves (without smoothing) indicated by shaded stripes in the respective top panels. (c) Power spectral densities (PSD) calculated for the sections of lightcurves (without smoothing nor normalisation) indicated by shaded stripes in the respective top panels. The dashed black lines indicate power-law indices
$-1$
and
$-2$
for reference.
In figure 4, we also present the temporal correlation analysis of the calculated lightcurves, using the method of discrete correlation function (DCF; Edelson & Krolik Reference Edelson and Krolik1988). The DCF method takes two series of measurements
$(t_{1,i},f_{1,i})$
and
$(t_{2,i},f_{2,i})$
, where time sampling can be irregular, bins all pairs of measurements according to time delay
$\unicode[STIX]{x0394}t=t_{2,j}-t_{1,i}$
, and evaluates the average discrete correlation
$\langle \tilde{f}_{1,i}\tilde{f}_{2,j}\rangle$
for each time delay bin using normalised measurements
$\tilde{f}_{i}=(f_{i}-\langle f\rangle )/\unicode[STIX]{x1D70E}_{f}$
, where
$\unicode[STIX]{x1D70E}_{f}$
is a dispersion of
$f_{i}$
. The result of DCF analysis is correlation coefficient as function of time delay
$\text{DCF}(\unicode[STIX]{x0394}t)\in [-1:1]$
. We calculate DCF(SYN,IC) between the synchrotron and IC lightcurves, as well autocorrelation functions
$\text{ACF}(X)=\text{DCF}(X,X)$
for either the synchrotron or IC lightcurves. The ACFs reveal characteristic variability time scales (measured between two zero points bracketing the main peak) of
$\simeq 0.8(L/c)$
in the 2Db_fc_kn and three-dimensional cases, and
$\simeq 0.4(L/c)$
in the 2Db_sc cases. The case 2Db_fc_th is more ambiguous, as the lightcurves in this case are barely resolved in time. No major differences are seen between the ACFs calculated for the synchrotron or IC lightcurves. With this insight from the ACFs, the DCF between synchrotron and IC lightcurves can tell us two things: the peak correlation coefficient, and the presence of any time lag between the two lightcurves. We find that the correlation coefficient for the cases with longer variability time scales (2Db_fc_kn and three-dimensional) is
${\geqslant}0.75$
, while in the cases with shorter variability time scales (2Db_sc) it is
${\leqslant}0.5$
. A minor time lag of
$+0.1(L/c)$
(synchrotron lagging the IC) is indicated only for the 3-D case, however, it is very unlikely to be statistically significant, as
$\text{DCF}(\unicode[STIX]{x0394}t=0)=0.87$
.
Also in figure 4, we present the variability statistics in the form of power spectral density (PSD) calculated separately for the synchrotron and the IC lightcurves. We demonstrate that in all ‘2Db’ cases, the PSD is consistent with a power law with index between
$-2$
(red noise) and
$-1$
(pink noise). However, in the 3-D case, the PSD is significantly steeper, with the index close to
$-3$
.
6 Perpendicular profiles of the current layers
Figure 5 shows perpendicular profiles across the current layers at the moment of saturation of coalescence instability, compared for simulations in the ‘2Da’ series with initial plasma temperatures spanning two orders of magnitude. We show the profiles of three parameters: number density
$n$
, mean particle energy (‘temperature’)
$\langle \unicode[STIX]{x1D6FE}\rangle$
and non-ideal electric field scalar
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
. Consistent with the findings of Nalewajko et al. (Reference Nalewajko, Zrake, Yuan, East and Blandford2016), the profiles of
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
are much thicker than the profiles of density or temperature. Here, we find that radiative cooling has no effect on the profiles of
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
. On the other hand, its effect on the temperature profile is seen for
$\unicode[STIX]{x1D6E9}>10^{6}$
. In the cases of very strong radiative cooling, the temperature profiles are flattened, approaching the thickness scale of the
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
. The profiles of plasma density are not consistent between the different cooling rates (i.e. the density peak values do not vary systematically with
$\unicode[STIX]{x1D6E9}_{e}$
), although they seem to preserve the same thickness scale.
Perpendicular profiles of current layers for simulations in the ‘2Da’ series with different particle temperatures
$\unicode[STIX]{x1D6E9}$
, and hence different radiative cooling time scales. The panels present: the particle number density
$n$
(a), the mean particle energy
$\langle \unicode[STIX]{x1D6FE}\rangle$
(b) and the
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
(c). For each
$\unicode[STIX]{x1D6E9}$
value, the profiles were extracted at the simulation time that maximises the amplitude of
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
.
Decomposition of the perpendicular profile of
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
according to the Vlasov momentum equation (6.1) compared for two simulations of the ‘2Db’ series: ‘fc_th’ (a) and ‘sc_th’ (b).
Figure 6 shows the decomposition of the perpendicular profile of
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
into components of so-called Vlasov momentum equation. We consider a phase-space distribution of particles (either electrons or positrons)
$f(\boldsymbol{x},\boldsymbol{p})$
that satisfies the Vlasov equation
$\text{d}f/\text{d}t=\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}t+v^{i}(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}x^{i})+F^{i}(\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}p^{i})=0$
, where
$F^{i}$
is any force acting on individual particles (in this case, it consists of Lorentz force and radiation reaction). A first-moment equation
$\int p^{i}(\text{d}f/\text{d}t)\text{d}^{3}p=0$
in scalar product with the magnetic field leads to the following form:
$$\begin{eqnarray}\displaystyle \left(E^{i}-\frac{\unicode[STIX]{x2202}u^{i}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}P^{ij}}{\unicode[STIX]{x2202}x^{j}}-Q_{\text{rad}}^{i}\right)B^{i}=0. & & \displaystyle\end{eqnarray}$$
Here,
$\boldsymbol{u}=n\langle \boldsymbol{p}\rangle$
is the momentum density,
$P^{ij}=n\langle p^{i}v^{j}\rangle$
is the pressure tensor and
$\boldsymbol{Q}_{\text{rad}}$
is a vector related to the radiation reaction force (detailed derivation will be presented elsewhere). It has been understood previously that gradients of the off-diagonal terms of the pressure tensor are important in supporting the non-ideal electric field in collisionless pair-plasma reconnection (Bessho & Bhattacharjee Reference Bessho and Bhattacharjee2007). Figure 6 shows that gradients of the pressure tensor balance
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
in the outer regions of the current layer. However, in the central regions of the current layer in the fast-cooling ‘fc_th’ case, there is a gap between the electric field and pressure terms, that is only partially filled with the momentum density term. The residual
$(E^{i}-\unicode[STIX]{x2202}_{t}u^{i}-\unicode[STIX]{x2202}_{j}P^{ij})B^{i}$
appears on the thickness scale of the temperature profiles (cf. figure 5). Since this residual is very small in the slow-cooling ‘sc_th’ case, we associate it with radiation reaction.
7 Discussion
In this work, we consider an IC process with the soft photons represented as a fixed uniform isotropic monoenergetic radiation field, corresponding to the ERC scenario characteristic for the luminous flat spectrum radio quasar (FSRQ)-type blazars (see the Appendix). We note that in order to include the synchrotron self-Compton (SSC) mechanism, we would need to calculate the production and propagation of multi-directional synchrotron radiation field, which would require much more complicated and expensive algorithms. Therefore, we focus here on a very specific application to the broadband emission from FSRQs.
7.1 Compton dominance
In all of our simulations, we have set the energy density of soft radiation fields at the level of mean magnetic energy density
$U_{\text{ext}}=\langle U_{\text{B}}\rangle$
. The idea was to make radiative cooling due to the synchrotron and IC processes equally important, at least as long as the IC scattering proceeds in the Thomson regime, and so the resulting luminosities should also be comparable, i.e. the Compton dominance parameter
$\text{CD}=L_{\text{IC}}/L_{\text{syn}}$
should be close to unity by design. In blazars, the Compton dominance is generally observed to range between 0.1 and 100 (see the Appendix).
It is hoped that the observed values of Compton dominance can be used to constrain the physical parameters of blazar jets, especially in the context of luminous FSRQ blazars where gamma rays are produced by Comptonisation of external radiation fields. With independent estimates of
$U_{\text{ext}}$
(from observations of broad emission lines or thermal emission of accretion disks), the value of
$\text{CD}$
provides an independent constrain on magnetic field strength (that determines
$L_{\text{syn}}$
) and hence the jet magnetisation (with total jet power constrained by
$L_{\text{IC}}$
) (Janiak et al.
Reference Janiak, Sikora and Moderski2015). However, this reasoning relies on the assumption that
$\text{CD}\simeq U_{\text{ext}}/U_{\text{B}}$
, which we are now testing in the case of electrons accelerated in relativistic magnetic reconnection.
Our simulated synchrotron and IC spectra show a range of behaviours leading to departures of Compton dominance from unity in both directions. In the case of IC scattering in the Klein–Nishina regime,
$\text{CD}<1$
from the beginning due to suppression of the scattering cross-section. As the IC spectral peak is dominated by contribution from cold electrons, the evolution of
$\text{CD}$
is driven primarily by evolution of the synchrotron peak. However, with IC scattering in the Thomson regime, the results are more nuanced. In the slow-cooling case (sc_th), energetic electrons produce a more prominent synchrotron excess, slightly suppressing the value of
$\text{CD}$
. On the other hand, in the fast-cooling case (fc_th), the high-frequency IC excess is more prominent, and the value of
$\text{CD}$
increases to
${\sim}10$
. We think that the main reason for this is that in the fast-cooling regime, energetic particles are strongly concentrated within current layers where (perpendicular) magnetic fields are much weaker and synchrotron losses are suppressed until the particles are able to exit the current layers. This separation between acceleration and radiation phases has been seen previously in Yuan et al. (Reference Yuan, Nalewajko, Zrake, East and Blandford2016), where only the synchrotron cooling was considered. However, if IC cooling is efficient (in the Thomson regime), it can easily dominate within the current layers. We should also note that this effect is not seen in our 3-D simulation, where cooling efficiency is comparable and IC scattering is more strictly Thomson. We have indication that the acceleration regions are less regular and they may not be able to trap energetic particles as well as in the 2-D case.
Variations of Compton dominance are observed even for individual sources when compared at different flux states, the highest values are typically observed during the most spectacular gamma-ray flares (Hayashida et al. Reference Hayashida, Nalewajko, Madejski, Sikora, Itoh, Ajello, Blandford, Buson, Chiang and Fukazawa2015). It has been previously suggested that the highest observed values of Compton dominance could result from localised destruction of magnetic fields by relativistic reconnection. Our results suggest that this could indeed happen (at least in the 2-D approximation) if energetic particles are confined to their acceleration regions where magnetic fields are decreased due to reconnection or forced to propagate along the magnetic field lines by parallel electric fields.
7.2 Variability of lightcurves
We have simulated synchrotron and IC radiation signals (lightcurves) produced by the same population of energetic electrons and calculated the corresponding discrete correlation functions (DCF). In the fast-cooling regime, we produce short simultaneous (time lags
${<}0.1L/c$
) flares of synchrotron and IC radiation. In the slow-cooling regimes, we obtain quasi-periodic oscillations with the period of
$\simeq L/(2c)$
, which reflects the periodic nature of the simulated system and the fact that particles crossing the periodic boundary normal to the observer line produce radiation signals separated by a delay of
$L/c$
. In the fast-cooling 2-D cases, the emitting regions can be to some degree isolated from the periodic boundaries. This is more difficult in the 3-D case, which features a larger number of magnetic reconnection sites that contribute to the observed signals. Hence, we do not find evidence for systematic time delay between synchrotron and IC signals produced by the same electrons.
The power spectral density (PSD) analysis reveals that in 2-D simulations the simulated signals are consistent with power laws with indices between
$-1$
and
$-2$
, as is generally observed in blazars (Appendix). However, the 3-D simulation produced smoother lightcurves with steep PSD of index
$\simeq -3$
. This signal could result from averaging comparable contributions from multiple emitting regions located in a relatively small simulation volume. The 3-D ABC configuration with
$k=2$
may simply be too regular to produce realistic radiation signals, as no individual reconnection site can dominate the observed signal.
7.3 Effects of radiative losses on development of current layers
Simulations of magnetic dissipation initialised with the ‘ABC fields’ are particularly valuable for studies of ab initio formation of dynamically evolving current layers. In the previous study (Nalewajko et al.
Reference Nalewajko, Zrake, Yuan, East and Blandford2016), we have found that such current layers develop two separate perpendicular thickness scales: a thinner density scale corresponding to the plasma skin depth, and a thicker
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
scale corresponding to the gyroradius of particles heated by the non-ideal electric field within the layer.
Our new results shown in figure 5 indicate that these separate thickness scales persist even under severe radiative cooling. We have also found evidence that the temperature thickness scale can be modified by strong radiative cooling, and hence it is not fundamentally related to the density thickness scale. Dynamical effects of strong radiative cooling are also found in the decomposition of the
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
profile across the current layer (figure 6), using the Vlasov momentum equation that anticipates the contribution of radiation reaction.
Effects of intense radiative cooling on the structure of reconnection current layers were investigated analytically by Uzdensky & McKinney (Reference Uzdensky and McKinney2011) in the framework of the classical Sweet–Parker model of resistive reconnection. Their main prediction is that of systematic compression of plasma within the current layer,
$A=n/n_{0}$
, which would effectively enhance the reconnection rate like
$v_{\text{rec}}\propto A^{1/2}$
. Our results obtained for the ‘2Da’ series of simulations and presented in figure 5 reveal neither a systematic increase in the plasma density with increasing
$\unicode[STIX]{x1D6E9}_{e}$
or decreasing
$l_{\text{cool}}$
, nor a systematic increase in the reconnection rate, an evidence for which would be revealed by the peak levels of the
$\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$
profiles.
We can formally estimate the value of predicted compression ratio
$A$
. Using equations (45)–(46) of Uzdensky & McKinney (Reference Uzdensky and McKinney2011) with our (2.13), the following expressions can be found for the volumetric energy conversion rates:
$Q_{0}=v_{\text{A0}}B_{0}^{2}/(4\unicode[STIX]{x03C0}L_{\text{layer}})$
and
$Q_{\text{rad}}=(4/3)\unicode[STIX]{x1D70E}_{\text{T}}cu^{2}U_{0}n=\unicode[STIX]{x1D6E9}_{e}nm_{e}c^{3}/l_{\text{cool}}\simeq cB_{0}^{2}/(8\unicode[STIX]{x03C0}l_{\text{cool}})$
, and hence
$A=Q_{\text{rad}}/Q_{0}\simeq (1/2)(c/v_{\text{A0}})(L_{\text{layer}}/l_{\text{cool}})$
, where
$v_{\text{A0}}\sim c$
is the Alfvén wave velocity in the relativistically magnetised background plasma. In our 2-D simulations initiated with diagonal ABC fields, typical half-length of the saturated current layer is
$L_{\text{layer}}\sim 300\unicode[STIX]{x1D70C}_{0}$
, hence we have at least nominally achieved conditions for
$A>1$
only for 2 simulations in the ‘2Da’ series for
$\unicode[STIX]{x1D6E9}_{e}\geqslant 3\times 10^{6}$
, for which
$l_{\text{cool}}<300\unicode[STIX]{x1D70C}_{0}$
(see table 1). Indeed, these are the two simulations, for which we find modified temperature profiles across the current layers.
Acknowledgements
We thank the anonymous reviewers for their helpful comments and suggestions. Zeltron is an open explicit parallelised particle-in-cell numerical code created by Benoît Cerutti and co-developed by Gregory Werner at the University of Colorado Boulder (http://benoit.cerutti.free.fr/Zeltron/). These results are based on numerical simulations performed at supercomputers: Mira at the Argonne National Laboratory, USA (2016 INCITE allocation; PI: D. Uzdensky) and Prometheus at Cyfronet AGH, Poland (PLGrid grant abcpic17; PI: K. Nalewajko). This work was supported by the Polish National Science Centre grant 2015/18/E/ST9/00580. Y.Y. gratefully acknowledges the support through a Lyman Spitzer, Jr. Postdoctoral Fellowship from the Department of Astrophysical Sciences, Princeton University.
Appendix. Synchrotron and IC radiation of blazars
The most direct application of this work concerns the broadband non-thermal emission from extragalactic astrophysical sources called blazars. Here we provide a brief summary of the physical parameters of their emission regions (see Madejski & Sikora Reference Madejski and Sikora2016 for a recent review).
Blazars are a subclass of active galactic nuclei, which are accreting supermassive black holes (of typical mass
$M_{\text{bh}}\sim 10^{8}{-}10^{9}M_{\odot }$
) residing in the centres of certain galaxies. The emission of blazars is strongly dominated by broadband non-thermal continua forming two main components: (i) the low-energy component extends from radio waves to UV/X-ray and is universally attributed to synchrotron radiation; (ii) the high-energy component peaks in the MeV/GeV/TeV gamma-ray band, and is generally attributed to leptonic IC radiation, although hadronic radiative processes are also being proposed (Sikora et al.
Reference Sikora, Stawarz, Moderski, Nalewajko and Madejski2009; Böttcher et al.
Reference Böttcher, Reimer, Sweeney and Prakash2013). The bolometric apparent luminosities of blazars range from
$10^{43}~\text{erg}~\text{s}^{-1}$
to
$10^{49}~\text{erg}~\text{s}^{-1}$
, exceeding the luminosities of their host galaxies (up to
$10^{45}~\text{erg}~\text{s}^{-1}$
). The luminosity ratio
$L_{\text{IC}}/L_{\text{syn}}$
of the IC and synchrotron spectral components is known as the Compton dominance. The observed values of Compton dominance range from 0.1–1 in the case of low-luminosity (BL Lac type) blazars up to 100 in the case of high-luminosity (flat spectrum radio quasars; FSRQs) blazars (Finke Reference Finke2013; Nalewajko & Gupta Reference Nalewajko and Gupta2017). At the same time, the non-thermal continua of blazars show persistent chaotic variability over time scales ranging from decades to minutes. The power spectral densities (PSD) of actual blazar lightcurves are consistent with featureless power laws with indices between
$-1$
and
$-2$
(Abdo et al.
Reference Abdo, Ackermann, Ajello, Antolini, Baldini, Ballet, Barbiellini, Bastieri, Bechtol and Bellazzini2010; Goyal et al.
Reference Goyal, Stawarz, Ostrowski, Larionov, Gopal-Krishna, Joshi, Soida and Agudo2017). Variability of the high-energy IC component (gamma-rays) is in general well correlated with variability of the low-energy synchrotron component (X-ray/optical), especially if one compares signals produced by electrons in the similar energy range. In some cases, the correlations are strictly simultaneous, as confirmed by the discrete correlation function (DCF) analysis (Wehrle et al.
Reference Wehrle, Marscher, Jorstad, Gurwell, Joshi, MacDonald, Williamson, Agudo and Grupe2012). For the above reasons, this non-thermal emission is attributed to relativistic jets streaming from the supermassive black holes with typical bulk Lorentz factors
$\unicode[STIX]{x1D6E4}_{j}\sim 5{-}40$
and (by chance) pointing towards the observer. Radiation produced isotropically within relativistic jets is very strongly beamed (the apparent bolometric luminosity scales roughly like
$L_{\text{app}}\sim \unicode[STIX]{x1D6E4}_{j}^{4}L_{\text{int}}^{\prime }$
), at the same time the variability time scale is shortened like
$T_{\text{app}}\sim T_{\text{int}}/\unicode[STIX]{x1D6E4}_{j}$
.
The extreme broadness of non-thermal blazar spectra requires a population of ultra-relativistic electrons reaching typical isotropic Lorentz factors (in the rest frame of bulk jet flow) of
$\unicode[STIX]{x1D6FE}_{e}\sim 10^{3}{-}10^{5}$
. A soft radiation field of co-moving photon energy
$E_{\text{soft}}^{\prime }$
will be upscattered to a gamma-ray photon of typical energy
$E_{\unicode[STIX]{x1D6FE}}^{\prime }\sim \unicode[STIX]{x1D6FE}_{e}^{2}E_{\text{soft}}^{\prime }$
. This means that the observed gamma rays are typically produced in a single scattering of a UV/optical/IR soft photon. The soft photons can be local synchrotron photons – such IC process is called synchrotron self-Compton (SSC) and is thought to dominate in the low-luminosity (BL Lac type) blazars. Alternatively, the soft photons may originate outside the relativistic jet – such an IC process is called external radiation Compton (ERC) and is thought to dominate in the high-luminosity (FSRQs) blazars. Luminous blazars are known to have a very rich radiative environment including direct thermal (UV/optical) emission of the accretion disk (quasar), broad emission lines (e.g. Ly
$\unicode[STIX]{x1D6FC}$
at
$\simeq 10~\text{eV}$
) and thermal IR emission of the dusty torus (
${\sim}0.3~\text{eV}$
).
It is typically inferred from detailed modelling of observed FSRQ spectra that they should be produced at distance scales in the range
$r\sim 0.1{-}10~\text{pc}$
from the black hole (Nalewajko, Begelman & Sikora Reference Nalewajko, Begelman and Sikora2014). The total inferred power of the underlying relativistic jet can reach Eddington luminosity
$P_{j}\sim L_{\text{Edd}}\sim 10^{46}~\text{erg}~\text{s}^{-1}$
. We can adopt a conical geometry with opening angle
$\unicode[STIX]{x1D6E9}_{j}\sim 0.3/\unicode[STIX]{x1D6E4}_{j}$
. The actual value of jet hot magnetisation
$\unicode[STIX]{x1D70E}_{\text{hot},j}=B_{j}^{\prime 2}/(4\unicode[STIX]{x03C0}w_{j})$
(where
$B_{j}^{\prime }$
is the co-moving jet magnetic field strength and
$w_{j}$
is the co-moving relativistic enthalpy density) in the dissipation/emission regions is uncertain. In fact, it is still being debated whether dissipation should be provided by shock waves at low magnetisations or by magnetic reconnection at high magnetisations, in any case one needs to distinguish between background/upstream and processed/downstream (see a discussion in Sironi, Petropoulou & Giannios Reference Sironi, Petropoulou and Giannios2015). We can assume that
$\unicode[STIX]{x1D70E}_{\text{hot},j}\sim 1$
and estimate the co-moving magnetic field strength from the magnetic jet power
$L_{B}=cR_{j}^{2}\unicode[STIX]{x1D6E4}_{j}^{2}B_{j}^{\prime 2}/8\sim L_{\text{Edd}}/4$
, with jet radius
$R_{j}=\unicode[STIX]{x1D6E9}_{j}r$
, hence
$B_{j}^{\prime }\sim (2L_{\text{Edd}}/c)^{1/2}/(0.3r)\sim 1~\text{G}(r/1~\text{pc})^{-1}$
. With such magnetic field strength, the characteristic electron gyroradius is
$\unicode[STIX]{x1D70C}_{e}=1.7\times 10^{7}(r/1~\text{pc})(\unicode[STIX]{x1D6FE}_{e}/10^{4})~\text{cm}$
, approximately 10 orders of magnitude smaller than the characteristic jet radius
$R_{j}\sim 10^{17}(r/1~\text{pc})~\text{cm}$
.
While we can estimate the relativistic specific enthalpy of the jet as
$w_{j}=B_{j}^{\prime 2}/(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}_{\text{hot},j})$
, the number density of electrons depends on the poorly understood composition of jet plasma. Relativistic jets can be composed of protons (ions) and ions and electrons with significant addition of electron–positron pairs (Sikora et al.
Reference Sikora, Stawarz, Moderski, Nalewajko and Madejski2009). If the energetic coupling between protons and leptons is strong, they can both participate in non-thermal particle acceleration and share the dissipated magnetic/kinetic energy. We can effectively parametrise this unknown microphysics with
$g_{p}=\langle \unicode[STIX]{x1D6FE}_{e}\rangle +\langle \unicode[STIX]{x1D6FE}_{p}\rangle (n_{p}m_{p}/n_{e}m_{e})$
, and then find the electron number density from relation
$w_{j}\simeq 4g_{p}n_{e}m_{e}c^{2}$
. The electron number density scales like
$n_{e}=B_{j}^{\prime 2}/(16\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}_{\text{hot},j}g_{p}m_{e}c^{2})\sim 5\times 10^{4}~\text{cm}^{-3}g_{p}^{-1}(r/1~\text{pc})^{-2}$
. The Thomson optical depth across the jet is
$\unicode[STIX]{x1D70F}_{\text{T},j}=\unicode[STIX]{x1D70E}_{\text{T}}n_{e}R_{j}^{\prime }\sim 0.003g_{p}^{-1}(r/1~\text{pc})^{-1}$
. Given that most likely
$g_{p}\gg 1$
, relativistic jets of blazars are optically thin at parsec scales.
