Probing the Emission Mechanism and Magnetic Field of Neutrino Blazars with Multi-Wavelength Polarization Signatures

TL;DR

This paper investigates the emission mechanisms and magnetic fields in neutrino blazars by analyzing multi-wavelength polarization signatures and performing simulations to distinguish leptonic and hadronic models.

Contribution

It introduces analytical models and RMHD simulations to connect polarization observations with blazar emission mechanisms and magnetic field properties.

Findings

High-energy polarization degree can differentiate emission mechanisms.

Neutrino flux correlates with magnetic field strength and radiation processes.

Polarization signatures reveal jet energy composition and magnetic field importance.

Abstract

The characteristic two-component blazar spectral energy distribution (SED) can be of either leptonic and/or hadronic origins. The potential association of the high-energy neutrino event IceCube-170922A with the flaring blazar TXS~0506+056 indicates that hadronic processes may operate in a blazar jet. Despite multi-wavelength follow-ups of the event and extensive theoretical modelings, the radiation mechanisms and the underlying magnetic field strength and configuration remain poorly understood. In this paper, we consider generic leptonic and hadronic blazar spectral models with distinct magnetic field strengths and radiation mechanisms. We analytically reproduce the SEDs and the neutrino flux of hadronic models, and predict their X-ray to $\gamma$-ray polarization degrees. Furthermore, by performing relativistic magnetohydrodynamic (RMHD) simulations taking into account the polarization-dependent radiation transfer, we study the time-dependent multi-wavelength polarization variability of the proton synchrotron model under a shock scenario. Our results suggest that the high-energy polarization degree and the neutrino flux can be jointly used to pinpoint the leptonic and/or hadronic blazar radiation mechanisms in the X-ray and $\gamma$-ray bands, and to infer the magnetic field strength in the emission region. Additionally, the temporal multi-wavelength polarization signatures in the proton synchrotron model shed light on the jet energy composition and the dynamical importance of magnetic fields in the blazar emission region. Future multi-wavelength polarimetry facilities such as {\it IXPE} and {\it AMEGO} together with neutrino telescopes such as {\it IceCube} can provide unprecedented observational constraints to probe the blazar radiation mechanisms and jet dynamics.

Tables

Table 1. Table 1: Model parameters and derived quantities of the four spectral fitting models shown in Figure 1 . The redshift measurement of TXS 0506+056 is taken from   Ajello et al. ( 2014 ) .

Model Parameters
Fitting Parameters	Model 1	Model 2	Model 3	Model 4
Redshift ($z$)	0.3365
Bulk Lorentz factor (${\Gamma }_{obs}$)	10
Viewing angle (${\theta }_{obs}$)	0
Blob size ($R^{\prime }$, cm)	$5\times {10}^{15}$
Escape time ($t_{esc}^{\prime }$, s)	$6.67\times {10}^5$
Magnetic field ($B^{\prime }$, G)	1.3	1.2	50	100
Electron minimal Lorentz factor (${\gamma }_{e,1}^{\prime }$)	2000	4000	400	250
Electron maximal Lorentz factor (${\gamma }_{e,2}^{\prime }$)	40000	40000	20000	10000
Electron power-law index ($p_e$)	2.0	2.0	3.2	3.2
Electron kinetic luminosity ($L_{e,kin}$, ${10}^{43}\mathrm{erg}{s}^{-1}$)	40	32	6	1.5
Proton minimal Lorentz factor (${\gamma }_{p,1}^{\prime }$)	–	1
Proton turnover Lorentz factor (${\gamma }_{p,b}^{\prime }$)	–	$2\times {10}^8$	$2\times {10}^8$	–
Proton maximal Lorentz factor (${\gamma }_{p,2}^{\prime }$)	–	$2\times {10}^9$
Proton power-law index before turnover ($p_{p,1}$)	–	2.1	2.1	2.2
Proton power-law index after turnover ($p_{p,2}$)	–	3.0	3.0	–
Proton kinetic luminosity ($L_{p,kin}$, ${10}^{46}\mathrm{erg}{s}^{-1}$)	–	3.2	8	9
Derived Quantities
Jet Power ($P_{jet}^{obs}=L_{e,kin}+L_{p,kin}+L_B$, ${10}^{46}\mathrm{erg}{s}^{-1}$)	0.0072	3.2	10.3	18.4
${\sigma }^{\prime }$	0.28	0.0004	0.29	1.04
$L_{e,kin}/L_{p,kin}$	–	$1.2\times {10}^{-3}$	$8.1\times {10}^{-6}$	$2.7\times {10}^{-6}$

Table 2. Table 2: Summary of parameters. Top: Conversion between the RMHD code units and the physical value. Bottom: Additional parameters used in the 3DPol simulation. All parameters in this table are in the comoving frame of the emission region, except for the bulk Lorentz factor.

RMHD Parameters	Code Value	Physical Value
Length	$L_0^{\prime }$	$2.0\times {10}^{16}\mathrm{cm}$
Time	$L_0^{\prime }/c$	$3.33\times {10}^{5}s$
Velocity	$c$	$3\times {10}^{10}\mathrm{cm}{s}^{-1}$
Magnetic Field	$B_0^{\prime }$	$100G$
Thermal Pressure	$B_0^{\prime 2}/\left(4\pi \right)$	$7.96\mathrm{erg}{\mathrm{cm}}^{-3}$
Plasma Density	$B_0^{\prime 2}/\left(4\pi c^2\right)$	$8.84\times {10}^{-19}g{\mathrm{cm}}^{-3}$

Equations

$${\gamma }_p^{\prime }\sim 10^6{\delta }_1^{-1},$$

$${ϵ}_{p\gamma }^{\prime }\approx \frac{{\overline{ϵ}}_{\Delta }}{2{\gamma }_p^{\prime }}=150{\gamma }_{p,6}^{\prime -1}\mathrm{eV},$$

$$\begin{array}[]{l}\dot{\gamma}^{\prime}_{p,p\gamma}=-c\sigma_{p\gamma}\frac{u^{\prime}_{rad}}{m_{e}c^{2}}\gamma_{p}^{\prime}\\ \dot{\gamma}^{\prime}_{p,syn}=-\frac{4}{3}c\sigma_{T}\frac{u^{\prime}_{B}}{m_{e}c^{2}}(\frac{m_{e}}{m_{p}})^{3}\gamma_{p}^{\prime 2}\end{array}~{}~{},$$

$$\frac{{\dot{\gamma }}_{p,p\gamma }^{\prime }}{{\dot{\gamma }}_{p,syn}^{\prime }}\sim 7\times 10^5\frac{u_{rad}^{\prime }}{u_B^{\prime }}{\gamma }_p^{\prime -1}.$$

$$f_{p\gamma }^{\prime }$$

$$\left(\epsilon_{\nu}^{2}F_{\epsilon_{\nu}}\right)^{\rm obs}\approx\frac{\delta^{4}}{4\pi d_{L}^{2}}\frac{3}{8}f_{p\gamma}^{\prime}\gamma_{p}^{\prime 2}m_{p}c^{2}\frac{d\dot{N^{\prime}}}{d\gamma_{p}^{\prime}}\big{|}_{E_{p}^{\prime}\approx 20\,\epsilon_{\nu}(1+z)/\delta}$$

$$f_{p\gamma }^{\prime }=\frac{t_{\mathrm{dyn}}^{\prime }}{t_{p\gamma }^{\prime }}=1.4\times 10^{-2}{R}_{\mathrm{BLR},16}$$

$$\Pi ({\nu }^{\prime })=Z_m\frac{{\Pi }_{max}({\nu }^{\prime })P_{pol}^{\prime }({\nu }^{\prime })}{P_{pol}^{\prime }({\nu }^{\prime })+P_{un}^{\prime }({\nu }^{\prime })},$$

$$\Pi ({\nu }^{\prime })=10\%\frac{{\Pi }_{max}({\nu }^{\prime })}{{\Pi }_{max,o}},$$

$$\begin{array}[]{l}B^{\prime}_{z}=B^{\prime}_{0}\times J_{0}(kr^{\prime})\\ B^{\prime}_{\phi}=B^{\prime}_{0}\times J_{1}(kr^{\prime})\end{array}~{}~{},$$

$$\begin{array}[]{l}n^{\prime}_{e}(\gamma^{\prime}_{e})=n^{\prime}_{e,0}\gamma_{e}^{\prime-2.0}e^{-\frac{\gamma^{\prime}_{e}}{10^{4}}}\\ n^{\prime}_{p}(\gamma^{\prime}_{p})=n^{\prime}_{p,0}\gamma_{p}^{\prime-2.0}e^{-\frac{\gamma^{\prime}_{p}}{5\times 10^{8}}}\end{array}~{}~{},$$

$$L_{\mathrm{Edd}}=4\pi GM{m}_{p}c/{\sigma }_T\sim 1.2\times 10^{47}{M}_{9}erg{s}^{-1}.$$

$$L_B=4\pi R^{\prime 2}c{u}_B^{\prime }{\Gamma }_{obs}^2\sim 4\times 10^{49}\mathrm{erg}{\mathrm{s}}^{-1},$$

Full text

galaxies: active — galaxies: jets — radiation mechanisms: non-thermal — relativistic processes

1 Introduction

Blazars are active galactic nuclei (AGN) whose jet is pointing very close to our line of sight (LOS). They are among the most powerful astrophysical objects in the universe, and exhibit highly variable emission from radio up to TeV $\gamma$-rays (Ackermann et al., 2016), indicating extreme particle acceleration. As a result, blazars have long been suspected as the extragalactic sources of ultra-high-energy cosmic rays (UHECRs). Blazar SEDs are typically composed of two parts: a low-energy component from radio up to soft X-ray and a high-energy component from X-ray up to TeV $\gamma$-ray. The low-energy component is generally believed to be the synchrotron emission of primary electrons, as it usually shows a high level of polarization (Scarpa & Falomo, 1997). The origin of the high-energy component however remains debated. In a leptonic scenario, the high-energy emission is produced by the same primary electrons upscattering seed photons through the inverse Compton (IC) process. The seed photons can be the low-energy synchrotron emission by the same population of electrons (synchrotron-self Compton, SSC, e.g., Marscher & Gear, 1985; Maraschi et al., 1992) or external photon fields such as the thermal radiation from the accretion disk, the broad line region (BLR), and the dusty torus (external Compton, EC, e.g., Dermer et al., 1992; Sikora et al., 1994). Alternatively, in a hadronic scenario, the high-energy component is produced by the primary proton synchrotron (PS) and/or the synchrotron of secondary charged particles from hadronic interactions (Mannheim & Biermann, 1992; Mücke & Protheroe, 2001; Petropoulou et al., 2015). Usually the hadronic model infers the acceleration of UHECRs in the blazar zone. So far the two models result in indistinguishable SED fittings (Böttcher et al., 2013; Cerruti et al., 2015). Time-dependent studies of the leptonic and hadronic models have shown some moderate differences in the multi-wavelength light curves (Li & Kusunose, 2000; Chen et al., 2014; Diltz et al., 2015; Zhang et al., 2016b). However, these features are hard to be quantified in practice even with simultaneous multi-wavelength coverage. Therefore, we need additional constraints to understand the blazar emission.

The production of high-energy neutrinos is the smoking gun of cosmic hadronic interactions. If blazars accelerate a considerable amount of high-energy protons, they may interact with the local strong blazar radiation field and produce charged pions, which then decay and emit neutrinos. The recent IceCube-170922A event was reported to be coincident with the blazar TXS 0506+056 flaring state (IceCube Collaboration, 2018a, b). Following the neutrino alert, the blazar was detected in multi-wavelength campaigns including very-high-energy $\gamma$-rays (IceCube Collaboration, 2018a; Ansoldi et al., 2018; Abeysekara et al., 2018). All wavelengths exhibit strong variability, with a $\sim 7\%$ optical polarization degree during flares reported by the Kanata Telescope (IceCube Collaboration, 2018a). Theoretically, many models have been proposed to explain the SED of the TXS 0506+056 flare and the IceCube-170922A event (e.g., Cerruti et al., 2019; Gao et al., 2019; Keivani et al., 2018; Ansoldi et al., 2018; Reimer et al., 2018). These models can be categorized into two generic groups: a leptonic setup with a subdominant hadronic component, where the IC dominates the high-energy emission but the subdominant hadronic component produces the neutrinos and makes considerable X-rays through hadronic cascading synchrotron; and a hadronic setup where the X-ray consists of both PS and cascading synchrotron but the $\gamma$-ray is dominated by the PS. Nonetheless, even though the two scenarios imply very different physical conditions (such as magnetic field strength) in the emission region, the above theoretical studies cannot unambiguously pinpoint the high-energy radiation mechanisms. In addition, these studies generally find that a simple one-zone model cannot fully explain the temporal behavior of the radiation and neutrino emission. Thus the underlying particle acceleration and the related magnetic field evolution are poorly constrained.

Magnetic field strength and evolution hold the key to disentangle various neutrino blazar models and shed light on the particle evolution. Polarimetry has been a standard probe of the blazar magnetic field evolution. In particular, the optical polarization signatures have been found highly variable alongside multi-wavelength flares, including significant polarization degree drops and $\gtrsim 180^{\circ}$ polarization angle swings (e.g., Abdo et al., 2010; Angelakis et al., 2016). Numerical simulations have shown that the amplitude of the optical polarization variability may indicate the magnetic energy percentage (often called the magnetization) in the blazar zone (Zhang et al., 2016a, 2017). In the high-energy bands, the PS scenario predicts a high polarization degree in X-rays and $\gamma$-rays (Zhang & Böttcher, 2013; Paliya et al., 2018), which may be detected by future high-energy polarimeters such as IXPE and AMEGO (Weisskopf et al., 2016; McEnery, 2017). The multi-wavelength polarimetry along with neutrino detection may unveil the details of the jet dynamics, particle acceleration, and radiation mechanisms.

In this paper, we make the first attempt to predict the multi-wavelength polarization signatures of both leptonic and hadronic blazar models, via semi-analytical calculation of broadband spectral polarization degree and first-principle-integrated polarized radiation transfer simulation of temporal polarization signatures. We aim to identify generic observable features that can be directly confronted with future multi-messenger and multi-wavelength polarimetry observations of blazars. In section 2, we take the TXS 0506+056 SED as an example, and apply one-zone spectral modeling of four benchmark leptonic and leptohadronic scenarios that correspond to drastically different jet physical conditions. Based on the fitting parameters, we predict the X-ray to MeV $\gamma$-ray spectral polarization degree that can be directly confronted with future high-energy polarimeters. In section 3, we perform 3D RMHD simulations with polarized radiation transfer to study the effects of the blazar zone magnetization on the temporal behavior of multi-wavelength polarization signatures of a PS blazar model. We demonstrate that the optical and MeV $\gamma$-ray polarization patterns can be very similar, and drastic temporal polarization variations, such as a polarization angle swing, indicate strong shock compression in a weakly magnetized emission region. We discuss the implications of our results on current and future observations in section 4, and summarize in section 5.

2 High-Energy Spectral Polarization Degree

In this section, we compute the high-energy spectral polarization degrees of four generic leptonic and leptohadronic blazar spectral models. Previous works have shown that the magnetic field is a key parameter that affects the contribution of different radiation mechanisms to the blazar high-energy spectral component (e.g., Böttcher et al., 2013; Cerruti et al., 2019; Reimer et al., 2018). We thus pick four benchmark parameter regimes with different radiation mechanisms and physical properties, especially the magnetic field strength and magnetization. To facilitate realistic predictions of the high-energy polarization, we take the recent TXS 0506+056 observation (IceCube Collaboration, 2018a) to derive the general physical parameters and particle distributions via one-zone spectral modeling. While we do not aim at detailed fitting, the derived parameters can help to illustrate their effects on the high-energy SED and spectral polarization. In the following, we first describe our models, then compute the broadband SEDs and neutrino flux, finally evaluate the spectral polarization degree.

2.1 Model Description

The key difference of the IC-dominated and PS-dominated blazar models generally lies in the magnetic field strength. In the former case, the magnetic field is usually on the order of $0.1~{}\rm{G}$, so that for most blazars the seed photon energy density for IC is comparable to or even larger than the magnetic energy density. Additionally, even if there exists an energetic proton population, they can only cool through pion production processes. Therefore, the high-energy emission is mostly dominated by the IC. In contrast, PS-dominated models usually require a magnetic field of at least $10~{}\rm{G}$. The IC seed photon energy density is then negligible compared to the magnetic energy density in most blazars, and protons can emit synchrotron efficiently. Given the above considerations, we set up the following four scenarios.

A pure IC-dominated model, where the high-energy spectral component is dominated by electron SSC. In principle, there could be a negligible hadronic component. 2. 2.

An IC-dominated model, where the high-energy spectral component is dominated by SSC, but it contains a small contribution from the hadronic component. 3. 3.

A PS-dominated model, where the high-energy spectral component is dominated by synchrotron of both protons and hadronic cascading pairs. There could be a small contribution from SSC of electrons as well. 4. 4.

A pure PS-dominated model, where the high-energy spectral component is dominated by proton synchrotron, with a minor contribution from the cascading pair synchrotron. Primary electron SSC is negligible.

Since TXS 0506+056 does not present a strong thermal component in the SED, we will not include a thermal background in our spectral fitting. However, in principle a small thermal flux could hide below the primary electron synchrotron, and may contribute to the photopion processes and the external Compton emission (EC).

One can immediately derive some general conclusions from the above four models. Model 1 is unlikely to produce any detectable neutrinos. Several previous works have found that Model 2 provides the best fit to the electromagnetic and neutrino emission from TXS 0506+056, although the modelled neutrino flux is significantly lower than the uppler limit inferred by IceCube-170922A (e.g., Cerruti et al., 2019). The major difference between Model 3 and 4 lies in the magnetic field strength, which results in different cooling rates.

In the following all quantities in the comoving frame of the blazar jet are marked with a prime, those in a frame that is stationary to the black hole are unscripted, and those in the observer’s frame are marked with a subscript “obs”.

Assuming that the blazar jet is relativistic with a bulk Lorentz factor $\Gamma_{obs}=10$, and the viewing angle is at $\theta_{\rm obs}=1/\Gamma_{obs}$, one can find that the observed neutrino energy of $E_{\nu}\sim 300~{}\rm{TeV}$ requires a parent proton energy of

[TABLE]

where $\delta_{obs}=[\Gamma_{\rm obs}(1-\beta_{\rm obs}\cos\theta_{\rm obs})]^{-1}=10\delta_{1}$ is the Dopper factor, with $\beta_{\rm obs}$ being the speed of the jet in unit of $c$. Then the target photon energy that corresponds to the $\Delta$-resonance of the photopion production is

[TABLE]

where we take $\gamma^{\prime}_{p}=10^{6}\gamma^{\prime}_{p,6}$ and $\bar{\epsilon}_{\Delta}\sim 0.3$ GeV is the resonance energy of the photopion production. The cooling rates of the photopion and synchrotron process are given by

[TABLE]

respectively, where $\sigma_{p\gamma}\sim 70\,\mu$b is the inelasticity-weighted average photopion cross section, $\sigma_{T}$ is the Thomson cross section, $u^{\prime}_{\rm rad}$ and $u^{\prime}_{B}$ are the photopion target photon energy density and magnetic energy density, respectively. Then one can find that

[TABLE]

The target photon energy density at $\epsilon^{\prime}_{p\gamma}$ is not well constrained by observations, but previous models have shown that an external radiation field is necessary in order to produce an observable amount of neutrinos (Cerruti et al., 2019; Reimer et al., 2018). If we adopt $u^{\prime}_{rad}\sim 100~{}\rm{erg\,cm^{-3}}$ as suggested by Reimer et al. (2018), then we can easily find that for moderate magnetic field strength, $B^{\prime}\lesssim 50~{}\rm{G}$, the photopion cooling rate is faster than the synchrotron for $\gamma^{\prime}_{p}\sim 10^{6}$, leading to efficient neutrino production.

2.2 Broadband SED

Our spectral model is based on a stationary lepto-hadronic radiation code developed by Böttcher et al. (2013), which semi-analytically treats radiative, photomeson, and adiabatic cooling as well as particle escaping. To facilitate direct comparisons, we choose the same bulk Lorentz factor, viewing angle, and emission region size for all four models. The default particle escaping time scale is chosen as 4 times of the light crossing time scale. The derived SED is corrected for the extragalactic background light attenuation. Figure 1 shows our fitting results with parameters listed in Table 1. We remind the readers that we do not aim at accurate fitting of the TXS 0506+056 observation. Instead, we are interested in placing general constraints on the parameter space of the four models, so that their key differences can stand out.

From Figure 1, we can identify the dominating radiation mechanisms in the high-energy bands for the four models. In Model 1, the primary electron SSC emission dominates the entire high-energy spectral component. Model 2 is similar to Model 1 but it has a significant synchrotron contribution by pairs in the X-ray band. These pairs are the byproducts of the photopion process, where the energetic protons interact with the primary electron synchrotron and produce neutrinos. We remind the readers that in the broadband SED calculation we do not include any external photon field for Compton scattering or hadronic interactions, because TXS 0506+056 is often considered as a BL Lac (however, some other works suggest that it could be an FSRQ, most recent one refer to Padovani et al., 2019). Generally speaking, however, the EC contribution can be significant in the $\gamma$-ray bands, especially for FSRQs. Both Model 3 and 4 have primary proton synchrotron dominating the high-energy bands, but in Model 3 the magnetic field is weaker. Since we do not include an external photon field, the target photon field for photomeson processes in Model 2-4 is the primary electron synchrotron emission.

2.3 Neutrino Flux

High-energy neutrinos are produced when primary protons interact with ambient photons in the jet. Using the measurement of the optical to X-ray flux of TXS 0506+056 we can derive the energy density of the target photons and the effective optical depth for the $p\gamma$ interaction (also referred as the pion production efficiency, $f^{\prime}_{p\gamma}$ below).

If neutrinos and $\gamma$ rays are produced in the same region, and target photons come from an optically thin area such as the emission region itself, we can obtain its density from the broadband SED. The pion production efficiency can be approximated as

[TABLE]

The fiducial value $\epsilon^{\prime}=200$ eV is chosen to be the target photon energy that yields the neutrino emission in the energy range observed by IceCube from TXS 0506+056. Assuming this target photon energies, the photopion production via the $\Delta$-resonance implies an emission of neutrinos with energies of $\sim 290$ TeV in the observer frame. $R^{\prime}$ and $t^{\prime}_{\rm esc}$ are the size of the emission region and the escape time of protons respectively. Their values are presented in Table 1.

Neutrinos carry $\sim 3/4$ of the energy of the charged pions, and about $1/2$ of the $p\gamma$ interactions lead to the production of charged pions. The neutrino flux can thus be estimated by

[TABLE]

where $d\dot{N}^{\prime}/d\gamma^{\prime}_{p}$ is the proton injection rate derived from Tabel 1. The neutrino flux of the four models is shown in Figure 1. Apparently, if the neutrinos are produced within the multi-wavelength flaring region, their flux should be much lower than the observed level, consistent with similar calculations done for example by Keivani et al. 2018.

To explain both the neutrinos and electromagnetic flaring emission observed at all wavelengths, below we also consider a two-zone model. Previous works have suggested that TXS 0506+056 may not be a typical BL Lac (e.g., Padovani et al., 2019), but it may have a weak external thermal component from the broad line region (BLR). Here the two zones refers to within the BLR (including the thermal emission from BLR, zone 1) and beyond BLR (without thermal emission from BLR, zone 2). Notice that this is the only place in the paper that we consider the contribution of external photon field, so as to explain neutrino production. We envision that the relativistic jet, containing high-energy protons and electrons, passes through the BLR and continues to move away from the central engine. In zone 1, the highly energetic protons interact with the dense BLR photon field and produce neutrinos through photomeson processes. The $\sim 100~{}\rm{GeV}$ $\gamma$-rays, however, can hardly escape during the neutrino production phase due to the large optical depth. Thus we do not expect to observe the $\sim 100~{}\rm{GeV}$ $\gamma$-ray emission accompayning the neutrino production. When these electrons and protons move to zone 2, the BLR photon field becomes negligible, and $\sim 100~{}\rm{GeV}$ photons start to escape. However, in zone 2 neutrino production is greatly reduced due to the lack of target photon field. Therefore, this transition through the two zones predicts a delay of the $100~{}\rm{GeV}$ $\gamma$-ray emission with respect to the neutrino event, while the radiation emitted at energies lower than $\sim 10~{}\rm{GeV}$ should appear simultaneous with the neutrino event.

Such two-zone model is appealing, even though it introduces several free parameters and diminish the model prediction power. Here we choose to estimate the neutrino production rate of this model by using the same physical parameters for both zones as in the above one-zone model, and introduce a BLR photon field in zone 1. The energy density of the line photons in BLR is $u_{\rm BLR}\approx 0.3\,\rm erg\,cm^{-3}$ in the black hole frame (Ghisellini & Tavecchio, 2008). The most prominent contribution comes from the Ly$\alpha$ at $\epsilon_{\rm UV}=2\times 10^{15}$ Hz. The photopion interaction rate in the comoving frame is ${t^{\prime}}^{-1}_{p\gamma}=\sigma_{p\gamma}\,c\,\Gamma^{2}\,u_{\rm BLR}/(\Gamma\,\epsilon_{\rm UV})$. The dynamical time for the jet to travel through the external field is $t^{\prime}_{\rm dyn}=R/\Gamma\,c$, thus the interaction efficiency is

[TABLE]

where $R_{\rm BLR}$ is the size of the BLR. This shows that the introduction of a BLR photon field indeed allows a higher photopion production efficiency that may explain the IceCube event.

To summarize, if neutrinos and VHE $\gamma$-rays are produced co-spatially, we expect that the actual neutrino flux level should be much lower than the estimated value suggested by IceCube Collaboration (2018a, b). But if they are produced in different regions, a strong external photon background would be allowed. The pion production efficiency could be greatly enhanced, leading to an average neutrino flux that is comparable to the IceCube measurement. However, in this model, the VHE $\gamma$-ray flare should appear later than the neutrino emission.

2.4 X-ray to $\gamma$-ray Polarization Degree

Given the fitting parameters and the particle distributions, we calculate the spectral polarization degree in the high-energy spectral component. We follow the procedure developed by Zhang & Böttcher (2013) to evaluate the spectral polarization. The calculation assumes a mono-directional magnetic field that is perpendicular to the LOS in the comoving frame of the emission region, and calculate the polarization-dependent radiation power parallel and perpendicular to the magnetic field for primary electron SSC, as well as synchrotron by primary electrons, primary protons, and secondary pairs from hadronic interactions. The evaluation is done by integration of the derived particle spectra from the SED fitting, thus it naturally includes all spectral changes and breaks derived from the fitting. Generally speaking, since the electrons are distributed isotropically in the comoving frame, the EC is unpolarized. In the end, by taking the EC component as unpolarized, the procedure obtains the maximal spectral polarization degree in all high-energy bands.

In practice, the magnetic field in the emission region is partially ordered, thus we need to correct the spectral polarization degree under a partially ordered magnetic field. Zhang et al. (2016b) have simulated the multi-wavelength polarization signatures in a helical magnetic field, concluding that the partially ordered magnetic field should lower the polarization degree in the optical and high-energy bands by the same ratio. Paliya et al. (2018) have used a simple approximation to include this effect, and we will follow it here. The general formalism for the observed polarization degree is

[TABLE]

where $Z_{m}$ is the correction factor for the partially ordered magnetic field, $\Pi_{max}$ is the maximal polarization degree, $P^{\prime}_{pol}$ and $P^{\prime}_{un}$ are the radiation powers of polarized components (such as SSC and synchrotron) and unpolarized components (such as EC and thermal), respectively, in a mono-directional magnetic field. Since we ignore all thermal and EC components in our spectral modeling, then $\Pi(\nu^{\prime})=Z_{m}\Pi_{max}(\nu^{\prime})$ for both optical and high-energy bands. Assuming a default observed optical polarization degree of $10\%$, then we find

[TABLE]

where $\Pi_{max,o}$ is the maximal optical polarization degree, which can be calculated based on the fitting parameters.

Figure 2 shows the predicted spectral polarization degree from X-ray to $\gamma$-ray energies for all four models. Since the SSC generally lowers the seed synchrotron polarization degree by half, Model 1 predicts a low polarization degree, $\sim 7\%$, throughout the high-energy spectral component. The sharp rise of the polarization degree in Model 2 is due to the strong secondary pair synchrotron in the X-ray bands. The polarization degree gradually drops towards higher energies, as the dominating radiation process transitions into the SSC. We want to emphasize that if the blazar has a strong EC component in the $\gamma$-ray bands, we expect that Model 1 and 2 should appear more weakly polarized or even unpolarized towards higher energies, as the EC component is unpolarized. Although Model 3 has a strong secondary pair synchrotron in the X-ray to soft $\gamma$-ray bands, 3 and 4 predict very similar high-energy polarization degree at $\sim 15\%$ level throughout the high-energy spectral component, except for the rising part in the X-ray bands. This rising part marks the transition where the proton emission processes become dominating in the high-energy spectral component, which is dependent on the exact fitting parameters for the PS component. The $\sim 15\%$ polarization degree for both models is easily understandable, as typically all synchrotron emission is similarly polarized, regardless of the type of the charged particles. We notice that generally all spectral polarization curves tend to rise towards higher energies. This is because the underlying particle spectra become softer at higher energies due to cooling effects, which boost the polarization degree. Clearly, PS-dominated models generally show higher polarization degree than IC-dominated models.

3 Temporal Multi-Wavelength Polarization Signatures

In this section, we study the effects of magnetization of the blazar emission region on the temporal behavior of multi-wavelength polarization signatures. As shown in the previous section, the PS models predict considerably higher X-ray and $\gamma$-ray polarization, which may be more interesting and easier to detect with future high-energy polarimeters such as IXPE and AMEGO. We limit our simulations in the PS setups. We consider the archetypal blazar flaring model where a shock propagates through the emission region and triggers strong particle acceleration. We choose two emission region setups, namely, a low magnetization environment with the magnetization factor $\sigma^{\prime}\sim 0.1$ (this factor is defined as the ratio of magnetic energy density over enthalpy in the fluid frame), and a high magnetization environment with $\sigma^{\prime}\sim 1$; everything else remain the same. To evaluate the time-dependent polarization signatures, we employ 3D RMHD simulations to self-consistently evolve the magnetic field during the shock propagation. In addition, we use 3D polarized radiation transfer simulations to calculate the observed polarization signatures. In the following, we first describe our simulation setup, then present the RMHD simulation results, finally calculate the temporal polarization signatures.

3.1 Simulation Setup

Our general simulation setup is very similar to Zhang et al. (2016a), except that now we include the PS in the radiation processes. Here we summarize the key model setups that are directly related to our simulations. Our model assumes that the blazar jet is traveling with a bulk Lorentz factor of $\Gamma_{obs}$ in the observer’s frame. In the comoving frame of the jet, the blazar emission region is initially cold (thermal pressure is trivial), pervaded by a helical magnetic field with a magnetization factor of $\sigma^{\prime}$. For simplicity, we assume that the initial magnetic field is under force balance with no pre-existing turbulence, and the plasma density and thermal pressure are uniform in the emission region. A flat shock of bulk Lorentz factor $\Gamma^{\prime}$ then propagates through the jet and encounters the emission region, which changes the plasma physical conditions through its interaction with the emission region and triggers local particle acceleration. After the shock leaves the emission region, we continue to simulate the jet evolution until when it reaches a quasi-stationary state. In this way, all the plasma physical quantities, in particular the evolution of the magnetic field and jet energy composition, are self-consistently simulated under first-principle magnetized fluid dynamics.

To produce the aforementioned model we use the 3D RMHD code LA-COMPASS (Li & Li, 2003). We perform the RMHD simulation in the jet comoving frame, with uniform Cartesian grids. The emission region is set to be a cubic box of length $5L^{\prime}_{0}$ at the center of the simulation domain (thus, $x,~{}y,~{}z$ range from $-2.5L^{\prime}_{0}$ to $2.5L^{\prime}_{0}$). The shock is set to be a flat layer (thickness $0.5L^{\prime}_{0}$) of fast-moving plasma initially with $\Gamma^{\prime}=5$. To avoid any artificial effects at the injection of this fast-moving layer, we put it at $-7.5L^{\prime}_{0}$ and let it propagate towards the emission region, thus its entering speed into the emission region can be slower. We pick the outflow boundary condition for the simulation domain. However, to get rid of any boundary effects, we choose a sufficiently large domain, $x,~{}y$ go from $-5L^{\prime}_{0}$ to $5L^{\prime}_{0}$, and $z$ goes from $-20L^{\prime}_{0}$ to $20L^{\prime}_{0}$, so that throughout our simulation, no signal from the emission region can arrive at the boundary. Within the emission region, we take a force-free helical magnetic field morphology in the form of

[TABLE]

where $J_{0}$ and $J_{1}$ are the Bessel functions of the first kind, $r^{\prime}$ is the radius measured from the $z$-axis, $B^{\prime}_{0}$ and $k$ are normalization factors of the magnetic field strength and the radius, respectively. The initial $\sigma^{\prime}$ is defined as the volume average magnetization factor in the emission region, with the low magnetization setup as $\sigma^{\prime}\sim 0.1$ and high magnetization setup as $\sigma^{\prime}\sim 1$. Given the above magnetic field morphology, the local magnetization factor is larger than the initial $\sigma^{\prime}$ at the central $z$-axis, and gradually drops outwards. In order to capture the detailed temporal and spatial evolution of the fluid physical quantities, we take a spatial resolution of $128\times 128\times 512$ in $x,~{}y,~{}z$ directions, with an output of plasma density $\rho^{\prime}$, velocity $\boldsymbol{v}^{\prime}$, magnetic field $\boldsymbol{B}^{\prime}$, and thermal pressure $P^{\prime}$ at every grid point in every $0.5t^{\prime}_{0}$ time step. Table 2 lists the conversion between code units and physical units.

RMHD simulations, however, do not include nonthermal particle acceleration and evolution. A detailed description of the nonthermal particle acceleration and cooling in our hadronic model is beyond the scope of this paper. For simplicity, we adopt the following assumptions to approximate the particle evolution. Given that many blazars can still emit radiation during quiescent states, we assume that the emission region has two populations of nonthermal particles for both electrons and protons. One population is a fixed and steady particle distribution $n^{\prime}_{e,b}$ and $n^{\prime}_{p,b}$, which is distributed uniformly in space. The underlying assumption is that this population of particles have reached equilibrium between secondary Fermi acceleration and radiative cooling. The other originates from the shock acceleration of fresh nonthermal particles, $n^{\prime}_{e,s}$ and $n^{\prime}_{p,s}$, which should be related to the time-dependent local shock properties. For both populations, we assume that the particle spectral distributions are of a power-law shape with an exponential cutoff,

[TABLE]

where $n^{\prime}_{e,0}=n^{\prime}_{e,b}+n^{\prime}_{e,s}$ and $n^{\prime}_{p,0}=n^{\prime}_{p,b}+n^{\prime}_{p,s}$ are normalized by the local physical quantities. Previous simulations have shown that the light curves and polarization signatures are mostly affected by the nonthermal particle density rather than the spectral shape (Zhang et al., 2016a, 2018), hence we fix the spectral index throughout our simulations. For $n^{\prime}_{e,s}$, given the very high magnetic field in the hadronic model, the electron cooling is very fast compared to the shock propagation or light crossing time scales. Therefore, we inject $n^{\prime}_{e,s}$ that is proportional to the local kinetic energy at the shock front in each RMHD time step, then in the next RMHD output time step, remove them and inject new $n^{\prime}_{e,s}$ based on the refreshed local kinetic energy and the shock front location. In this way, we mimic the shock acceleration of electrons and their fast cooling. The proton injection $n^{\prime}_{p,s}$, however, needs a more careful treatment. Even the high-energy protons that emit $\gamma$-rays ($\gamma^{\prime}_{p}\sim 10^{8}$) have cooling time scales comparable to the light crossing time scale of the emission region. Therefore, the protons cannot be simply refreshed at each RMHD output time step. Given the very long cooling time scale, we assume that the nonthermal proton energy density $n^{\prime}_{p,s}$ is proportional to the local thermal energy density in the emission region. Zhang et al. (2016b) have demonstrated that if the proton cooling time scale is comparable to the light crossing time scale, as in our simulation setup, this approximation can work. However, in the case that the light crossing time scale is even shorter, the detailed proton cooling physics must be carefully evaluated, which is beyond the scope of this paper.

Finally we calculate the time-dependent polarization signatures based on the above magnetic field evolution and particle distributions. The most important physical effect here is the light crossing time scale. We use the 3DPol code to include all light crossing time effects for both electron and proton synchrotron emission (Zhang et al., 2016b). The light crossing time effect is a physical effect that photons from different location in the emission region may arrive to the observer at different time, which leads to various arrival time delays. The 3DPol code takes Cartesian coordinates as in the RMHD simulation, and evaluate Stokes parameters in the comoving frame at each grid point at each time step. Then it applies ray-tracing method to propagate the beam to the plane of sky. There it adds up all the beams that arrive at the same time, and Doppler boosts all radiation and polarization signatures to the observer’s frame. Due to the relativistic aberration, for blazars observed at $\theta_{obs}\sim 1/\Gamma_{obs}$ from the jet axis, in the comoving frame we are observing at $\theta^{\prime}\sim 90^{\circ}$. For all simulations we take the LOS at $\theta^{\prime}=90^{\circ}$, hence the Doppler factor is $\delta_{obs}=\Gamma_{obs}$. We define the polarization angle $PA=0$ when the electric vector is perpendicular the jet axis.

3.2 Shock Interaction with Magnetic Field

When the shock passes through the emission region, we observe that shocked plasma shows higher pressure and density due to compression (Figure 3 third column). In ideal RMHD, since the magnetic field lines are frozen in the plasma, the toroidal magnetic field component is also enhanced in the shocked plasma. Given the helical magnetic field setup in Equation 10, the toroidal component is at its maximum at radius $r^{\prime}\sim 1.2$ from the $z$-axis, thus we observe strengthened toroidal component near this radius (Figure 3 first column). An interesting feature is that the poloidal component is stronger at the central jet (Figure 3 second column). Additionally, one may quickly notice that the shocked plasma is not flat. Both phenomena are due to the magnetic force triggered by the strengthened toroidal magnetic component in the shocked plasma. As in both simulation setups the emission region is sufficiently magnetized, this magnetic force, which is in the radial direction, pinches the shocked plasma, squeezing the central jet ($r^{\prime}<1$) while expanding outwards ($r^{\prime}>1.5$) radially (Figure 3). Due to the frozen-in magnetic field lines, the squeezed central jet plasma has stronger poloidal component in a bullet shape. Meanwhile, the expansion traces a curved shocked plasma far downstream the shock front. The strength of this radial magnetic force is proportional to the local magnetization. From Figure 3, we can clearly see that the poloidal enhancement is stronger and the shape of the shocked plasma changes more in the high magnetization case than the low magnetization case.

Comparing the low and high magnetization cases, several differences are obvious. One is the magnetic field enhancement in the shocked plasma. Clearly, the toroidal component is more increased and the poloidal component is less enhanced in the low magnetization case. Given that our initial magnetic morphology has a higher poloidal component, we expect a stronger magnetic morphology change in the low magnetization case. The other difference is that the volume of the shocked plasma with enhanced thermal energy and higher velocity is much larger in the low magnetization case. Then, based on our particle injection method, one should expect that in the low magnetization environment both protons and electrons are distributed in a larger region where the magnetic field morphology is also different from the initial conditions. This should increase the polarization variability as well.

3.3 MeV Polarized Variability

The general polarized variability follows well with the magnetic evolution during shock propagation. Here we only plot the MeV light curves and polarization signatures. This is because the X-ray band emission should have a nontrivial contribution from the cascading pair synchrotron, which may alter the general polarization patterns. For higher energies, the general behaviors should be almost identical to the MeV band. Before the shock enters the emission region, since the initial magnetic morphology has a stronger poloidal component than the toroidal component, and the LOS is at $90^{\circ}$ from the jet axis in the comoving frame, the polarization angle stays at 0 along the jet axis with $20\%$ polarization degree. When the shock moves in, the toroidal component is strongly compressed. The same region also has more nonthermal particles due to the shock acceleration, leading to a flare in the light curve 4. Therefore, the toroidal contribution quickly takes over the poloidal contribution, and the polarization angle moves toward $90^{\circ}$. When the shock leaves the emission region, since the toroidal component is frozen in the shocked plasma, the downstream toroidal component is stretched and weakened due to the shock propagation. On the other hand, the poloidal component is hardly affected by the shock propagation. Therefore, the polarization reverts to an even more poloidal dominating situation than the initial state. Also due to the weakened toroidal magnetic field, the flux slightly decreases after the shock.

We can see that the optical and MeV light curves and polarization signatures are nearly identical. This is consistent with Zhang et al. (2016b), where they find that if the light crossing time scale is comparable to or longer than the proton cooling time scale, then the corresponding high-energy polarized variability should be similar to optical counterparts. The only difference between the optical and MeV curves in Figure 4 is that in the high magnetization case the MeV polarization angle does not complete a polarization angle swing. This is because a complete $180^{\circ}$ polarization angle swing requires strong shock compression in the toroidal component, but in the high magnetization case the shock compression is not strong enough. In the optical band, since the shock injected nonthermal electrons strictly follow the magnetic compression at the shock front, the toroidal contribution is just enough to push a complete polarizaiton angle swing. However, the nonthermal protons are proportional to the thermal energy distribution, which covers a larger downstream region where the toroidal component is less dominating (Figure 3), the toroidal enhancement during the shock is not sufficient to push a complete polarization angle swing.

Comparing the low and high magnetization cases, we can easily see that the low magnetization case shows stronger flares and larger polarization variations in both optical and MeV bands. This is straightforward, as the shock compression is less efficient in a more magnetized environment. An interesting feature is that after the shock leaves the region, the polarization degree stays at a high level in the low magnetization case, but in the high magnetization case it can revert back to the initial value. This is because the magnetic force tries to restore a force-balance state. After the shock leaves the emission region, the poloidal component is strong in the central spine. This results in an outward radial magnetic force. If the emission region is sufficiently magnetized, this force can quickly expand the central spine and restore a force-free magnetic field topology, similar to the initial state. However, in the low magnetization environment, this force is relatively weak, thus the dominating poloidal component in the central spine cannot be relieved quickly, leading to a high polarization degree after the shock. Therefore, the restoration of the polarization signatures can help to diagnose the magnetization in the emission region.

4 Implication on Observation

Combining multi-wavelength polarimetry and neutrino detection, we can pinpoint the dominating radiation mechanisms in all high-energy bands. In a pure leptonic model such as Model 1, we do not expect any neutrino detection, and the high-energy polarization should be at most half as polarized as the optical bands. If one observes a decrease of the polarization towards higher energies, it infers the transition from SSC to EC processes. However, if one finds that the X-ray is similarly polarized as the optical band but the $\gamma$-ray polarization is much weaker, then it indicates a IC-dominated model with a cascading pair synchrotron component as in Model 2. We also expect to see neutrinos in such blazars. If one detects similar polarization degree in both X-ray and $\gamma$-ray bands at the level of the optical polarization, then the high-energy spectral component is dominated by the PS mechanism. However, PS models with different magnetic field strengths predict similar high-energy spectral polarization and neutrino flux. Therefore, we need time-dependent polarization signatures to distinguish them.

Two caveats, however, need to be mentioned. One is that if the blazar is an FSRQ that has a strong thermal component (the “big blue bump”), then one should use the infrared to red polarization degree as the optical polarization to avoid thermal contamination. The second caveat is that if the blazar is a high-synchrotron-peaked so that the primary electron synchrotron extends to X-ray bands, then only the $\gamma$-ray polarization constraints can be employed.

Based on our MHD-integrated time-dependent polarization simulation, we can observe that the PS model with moderate magnetization factor ($\sigma^{\prime}\sim 0.1$) predicts much stronger polarization variations in both degree and angle then the high magnetization case ($\sigma^{\prime}\sim 1$). This applies to both optical and MeV $\gamma$-ray bands. We expect $\sim 20\%$ polarization degree during the flare and $>40\%$ after the flare for the low magnetization environment. For the high magnetization case, the polarization degree stays below $\sim 30\%$ throughout the simulation, with $<10\%$ during the flare. Given that the TXS 0506+056 event has an observed optical polarization degree about $7\%$ during the flare (IceCube Collaboration, 2018a), it is likely that the emission zone has a high magnetization.

A very interesting phenomenon is that the MeV band can also present polarization angle swings similar to the optical band. Recent observations have found many optical polarization swings in blazars, and mostly they are accompanied by strong $\gamma$-ray flares (e.g., Blinov et al., 2016; Angelakis et al., 2016; Blinov et al., 2018). In general, both random walks due to turbulence and large-scale magnetic field variations due to shocks or magnetic instabilities can explain these events (e.g., Marscher, 2014; Zhang et al., 2016a; Tavecchio et al., 2018). However, in hadronic blazars, the proton cooling time scale is generally much longer than the random walk time scale, which is determined by the magnetic coherent length scale in the turbulence. As a result, the slow proton cooling should smooth out the stochastic behaviors of high-energy polarization signatures. Therefore, detection of high-energy polarization angle swing infers a major variation of the blazar zone magnetic morphology.

5 Summary and Discussions

In this paper, we predict the high-energy spectral polarization degree as well as their corresponding magnetic field strength and radiation mechanisms, based on detailed one-zone steady state spectral modeling. We find that the PS-dominated models generally predict higher polarization degree than the IC-dominated models across the entire high-energy spectral component. Also we find that the high-energy spectral polarization degree together with the neutrino detection can pinpoint the radiation mechanisms and magnetic field strength in the emission region. Additionally, we evaluate the temporal behaviors of multi-wavelength polarization signatures in PS-dominated models based on combined RMHD and polarized radiation transfer simulations. Our results suggest that the high-energy light curves and polarization signatures are very similar to the optical counterpart, except that the polarization signatures are slightly less variable. An interesting phenomenon that our work predicts is that the polarization signature in the high-energy band may show similar behaviour to the optical one, with polarization angle swings symptom of strong changes in the magnetic field morphology. However, this feature can only be observed with dedicated high-energy polarimeters with good temporal resolution. Furthermore, the time-dependent multi-wavelength polarimetry can shed light on the magnetic field evolution in the emission region. Future high-energy polarimeters such as IXPE and AMEGO are promising tools to test these features and put unprecedented constraints on the blazar emission, especially during their active phases.

Our semi-analytical study and numerical simulations aim at identifying key observational features and predictions of blazars with consistent a description of the underlying physics. However, in practice, there exist several important physical quantities that may impact our predictions. The first one is the size of the emission region. As shown in Table 1, the total jet power in the PS-dominated model is considerably higher than the Eddington luminosity,

[TABLE]

In fact, we can easily estimate the magnetic power in our hadronic model in section 3,

[TABLE]

which by itself already exceeds the Eddington luminosity. Although occasionally it is possible that the emission region power can be higher than the Eddington luminosity, in general one expects that the former should be smaller than or comparable to the latter. Then the size of the emission region should be at least one order smaller than our numerical simulation parameters (Table 2). In this case, however, the proton cooling time is much longer than the light crossing time. This means that for any jet energy dissipation and particle acceleration processes that take place in the emission region, the newly accelerated protons can linger over the entire duration of the flare. Zhang et al. (2016b) has argued that these slow-cooling protons may result in nearly steady high-energy polarization signatures during flares; however, they assume that the magnetic field can revert to its initial state immediately after the passage of the dissipation process, which is not the case based on our RMHD simulations. On the other hand, for a large emission region such as in this paper, our results are consistent with Zhang et al. (2016b). Therefore, we suggest that a full description of the time-dependent high-energy polarization requires at least the combination of fluid dynamics, radiative cooling, and polarized radiation transfer altogether.

Magnetic instabilities and/or turbulence can also efficiently dissipate jet energy to accelerate particles through magnetic reconnection events (e.g., Giannios et al., 2009). In particular, recent particle-in-cell (PIC) simulations have demonstrated that reconnection can accelerate protons as well (Guo et al., 2016). Our spectral modeling parameters suggest that the magnetization factor in the PS-dominated models is generally between $0.1-1$. This indicates that hadronic blazar flares may be even preferably triggered by magnetic energy dissipation rather than shocks. Several previous works have studied the optical polarization signatures under kink instability and magnetic reconnection (Zhang et al., 2017; Nalewajko, 2017; Zhang et al., 2018). However in these cases the magnetic field and nonthermal particles are co-evolving, which require careful treatment of the radiative cooling as well.

The last issue is that all our calculations assume a laminar setup. In practice, turbulence may widely exist in the blazar zone. Given that the gyroradius of the high-energy protons in the hadronic model can be much larger than that of the electrons, if the magnetic coherent length scale is smaller than the proton gyroradius, then protons can experience more turbulent field than electrons, leading to lower high-energy polarization degree than the optical bands. However, notice that on average the turbulent magnetic field should also partially cancel out the optical polarization. Therefore, we expect that overall the polarization degree should be similar in optical and high-energy bands for hadronic models, but optical bands may exhibit some fast high polarization spikes due to the small magnetic coherent lengths. Such stochastic behaviors of the optical polarization have been found in some leptonic models such as Marscher (2014); Tavecchio et al. (2018).

We acknowledge helpful discussions with B. Dingus. HZ acknowledges support from Fermi Guest Investigator program Cycle 10, grant number 80NSSC17K0753, and Cycle 11, grant number 80NSSC18K1723. KF acknowledges the support of a NASA Hubble Fellowship Program Einstein fellowship. HL acknowledges the support by the LANL/LDRD program, NASA/ATP program, and the DoE/OFES program. HZ and DG acknowledge support from the NASA ATP NNX17AG21G and the NSF AST-1816136. MB acknowledges support by the South African Research Chairs Initiative (SARChI) of the Department of Science and Technology and the National Research Foundation111Any opinion, finding and conclusion or recommendation expressed in this material is that of the authors and the NRF does not accept any liability in this regard. of South Africa.