Deep learning model for multiwavelength emission from low-luminosity active galactic nuclei

We assumed sub-Eddington LLAGNs are accreting as radiatively inefficient accretion flow (RIAF). We used a semi-analytical approach to treat the RIAF emitted radiation (Nemmen et al., 2014). Normally LLAGN SEDs are radio-loud (Ho, 1999; Sikora et al., 2007). In our model, we also included the synchrotron emission from a relativistic jet, either. Here in this work, we followed the same approach as Nemmen et al. (2014) to generate our SED sample. ¹¹1Fiducial model source code can be found in: https://bitbucket.org/nemmen/adaf-code/src/master/. Henceforth we will call this model as the “fiducial” one.

Our model is a deep neural network composed of neurons with several layers and weights. In the output layer, we have a function called loss function $\mathcal{L}$ which is any error function that gives us how much the prediction and target are different. The learning procedure is given by the backpropagation method that consists of derivatives of $\mathcal{L}$ concerning the weight we want to learn (Kelley, 1960). In other words, the backpropagation uses gradient descent (Ruder, 2016) to find the minimum value in a space given by $\mathcal{L}(w_{ij})$, i.e., find which weights $w_{ij}$ gives the minimum of $\mathcal{L}$.

The learning part starts when we feed the neural network with data from the training set. We can give the data in batches, which is giving data by parts; this corresponds with faster training. Epoch is passing all batches once. The batch size and epoch are called hyperparameters. They are changed to find which combination of hyperparameters gives better results. In our work, we considered the number of neurons and the number of layers as hyperparameters. In this way, our goal was to find the best architecture for our model. We created a GridSearch to find the best values of our hyperparameters. The GridSearch is a hyperparameter tuning method (Bergstra & Bengio, 2012) that consists of train the NN with several combinations of hyperparameters. In our case, our hyperparameter is the number of neurons in each layer. We created a shell script that runs all the combinations possible and saves the chosen metric. The trainings with the better metrics are the one we use to further analyze. Our best model is composed of two architectures, one is to analyze the RIAF and the other to analyze the jet region. The RIAF model is composed of 4 layers with 56, 60, 99, and 99 neurons, respectively. The jet model is composed of 5 layers with 10, 44, 66, 99, and 130 neurons. Each layer has a ReLU activation except for the last layer which is a linear function. The mean absolute error as $\mathcal{L}$ and the optimizer is Adam (Kingma & Ba, 2014).

To summarize, we build our neural network from scratch using the GridSearch technique. Our architecture is composed of hidden layers meaning that it is a deep neural network. We tested other architectures with different hyperparameters, but this is the best model.

SED calculations for LLAGN are time expensive. The fiducial code which we generated our training RIAF SED’s sample took a minute per spectrum²²2Depends on the hardware and parameters set. For our available computers, it took 0.5-2min. The trained NN, AGNNES, calculates the same component in approximately $0.3$ ms. It is $\sim 400000$ times faster than the original RIAF code. The time improvement allowed us to implement more robust statistic methods to LLAGN SED fitting, which is absent in previous works. The original jet calculation is much faster than the RIAF component, took around $\sim 0.2$ seconds each jet SED. AGNNES calculates the same component in approximately $0.4$ ms. The improvement, for the jet, is close to $\sim 500$ times.

AGNNES results are an approximation from the original calculation. We had an extremely high speed-up in calculations, but the NN introduces inherently small errors when compared to the fiducial code. AGNNES, on average, can reproduce very well the validation data – see appendix C. Comparing the AGNNES predictions with the validation data sample, we defined $\Delta{\rm SED}$ as the “distance” between the model and AGNNES prediction:

Considering the validation sample and comparing original data against AGNNES predicted one, for every point (i.e. frequency). We found the averaged value of $\overline{\Delta{\rm SED}}$ = $0.02\pm 0.05$ (see (2)), we plotted in figure 1. For the RIAF component the NN predicted slightly smaller emission values. The jet component presented $\overline{\Delta{\rm SED}}$ = $0.00\pm 0.01$. These $\overline{\Delta{\rm SED}}$ values are AGNNES’s uncertainty.

The original data is not perfect. Sometimes the original code randomly fails while calculates the inverse Compton emission, returning a lower value for the emission. We automatically generated the data, and we identified the code failures and discarded the wrong data automatically, but it was not 100% efficient. The data set have few“bad” SEDs mixed in the whole set; this can affect the NN training and be responsible for the slightly systematic deviation of 0.02 dex in $log_{10}(\nu L_{\nu})$. However, LLAGNs are very difficult to observe, and the observations typically have considerable uncertainties, higher than 0.02 dex. The theoretical models for the LLAGN emission have their uncertainties. The NN averaged deviation from the original model is negligible if compared with the usual error bars in measurements and model uncertainties.

The current approach for LLAGN SED fitting in literature is an iterative method in which one changes all the parameters individually while the others remain unchanged. The best fit is determined visually after a certain number of attempts with no estimated uncertainties on the results (Nemmen et al., 2014; Almeida et al., 2018). The SED calculation is computationally expensive. The bottleneck of LLAGN SED calculation is in the Comptonization. The fiducial code takes $\sim 2$min to generate a single SED. It is tough to explore a vast range of parameters; there is insufficient time for this. For example, to run $10^{5}$ SEDs, we would need approximately two months without pause. With a faster method to calculate SEDs, we could explore a broader range of parameters and apply more robust fitting methods.

Time performance is the strongest point of AGNNES when we compare it with the iterative method. A well trained NN can do the IC step smoothly without high time spending. The NN is not a perfect copy of the original model. During the training, there is error propagation, and the NN final result is an approximation of the real value; for simplicity, we assumed the original code gives the perfect SED value, and we compared AGNNES with it. We discussed the NN performance in section 2.3.

We implemented the fitting procedure using Markov-chain Monte Carlo method. We performed the chain calculation with the Python package emcee (Foreman-Mackey et al., 2013). We used emcee to estimate the posterior distributions for the parameters – $\delta$, $s$, $\dot{M}_{0}$, $p$, $\dot{M}_{j}$, $\epsilon_{e}$, $\epsilon_{B}$ – that better describe the data. We defined the likelihood function as a pure Gaussian likelihood. The final likelihood is the sum of all seven parameters’ likelihoods. The priors of parameters have initially been the limits of the data set presented in sections 2.1.1 and 2.1.2.

We started the MCMC with a flat distribution ball around random initial values with small radius, considering all parameters as adimensional. Our MCMC chain ran with 300 walkers for the number of steps N: $N>30000$. We neglected the first 20% steps as burn-in. The value of N depends strongly on the used dataset. We estimated the MCMC convergence using the integrated autocorrelation time ($\tau$) implemented in emcee following Goodman & Weare (2010) method; it should be $N/50\gtrsim\tau$.

The elliptical galaxy M87 harbors a supermassive black hole of $6.5\times 10^{9}M_{\odot}$ at 16 Mpc of distance from us (Event Horizon Telescope Collaboration et al., 2019). The source presents strong radio emission due to a prominent relativistic jet. We showed the SED data points and the best fit in figure 2.

For this source, we used as data points the high-resolution observations presented in Prieto et al. (2016) (c.f. Table 4). We treated observations with the lower spatial resolution– presented in the same paper for different wavelengths –as upper limits.

The jet emission is a known feature of M87 (Doeleman et al., 2012). We tested models considering RIAF and Jet emission, we found the best result considering only the jet emission. Our fit suggested jet dominance over the inner RIAF for all wavelengths. It was not possible to produce a good fit for the observational SED using the RIAF model for this source. Our result considered only jet synchrotron emission and did not reproduce the bump around 100GHz and the most energetic data point is consistent with the best fit with $3\sigma$.

A thermal population of electrons in the accretion flow cannot reproduce the SED of M87 with such accretion rate. It is necessary to change the electronic distribution or assume the jet dominates the observed energy output entirely. The results for M87 fitting parameters are in table 1, and we only considered the synchrotron jet in our model.

NGC 4261 is an elliptic galaxy at $\sim 30$Mpc of distance, located at the Virgo cluster. It harbors a SMBH with $5.2\times 10^{8}M_{\odot}$ (Tremaine et al., 2002). The SED data is available in table 3. de Menezes et al. (2020) presented data points and similar modeling for the same galaxy. We showed the SED data points and our best fit in figure 3.

NGC 315 is another elliptical galaxy at $\sim 64$Mpc of distance, located at Pisces constellation. It contains a SMBH of $7.8\times 10^{8}M_{\odot}$ (Woo & Urry, 2002). The SED data is available in table 4. de Menezes et al. (2020) presented data points and similar modeling for the same galaxy. NGC 315’s SED data points and the best fit are in figure 4.

For these sources, we considered the optical data as upper limits. In the optical frequency range, there is a high possibility of stellar contamination. Stellar clusters typically inhabit the center of galaxies, and these populations emit a considerable amount of light in the optical band (Lauer et al., 2005). Since our primary interest here is the AGN emission, we neglected the potentially contaminated points by stars.

For NGC 4261 and NGC 315, our models were more focused on the RIAF component, especially for higher frequencies. In our model, the innermost regions of the RIAF produce most of the X-ray emission. We used the jet component to explain mostly the radio region of the SED. We did not consider the jet emission to explain the observed high energy emission (optical, UV, or X-ray bands). If we try to explain the whole SED only with a synchrotron jet, we will get different results and higher values of $\chi^{2}_{\rm red}$.

For both NGC 4261 and NGC 315,– assuming the RIAF dominates high energy emission –there are no strong constraints in the high-frequency jet emission ($\nu>10^{12}\text{Hz}$) –, this assumption makes the value of $p$ poorly constrained. The jet model fits the radio points, and this gives us information about the accretion rate $\dot{M}_{j}$ and the distribution of energy, represented by $\epsilon_{e}$ and $\epsilon_{B}$. In our model, we restrained $\dot{M}_{0}<10^{-2}\dot{M}_{\rm Edd}$. We trained AGNNES to calculate SEDs with accretion rates up to $10^{-1}\dot{M}_{\rm Edd}$, but we expect even lower accretion rates for LLAGNs. Our assumption of the maximum value for $\dot{M}_{0}$ resulted in the truncated distributions for $\dot{M}_{0}$ in NGC 315 and NGC 4261, respectively figures 8 and 7.

If the sources have higher accretion rates, the RIAF assumption is not valid. This truncation can be a hint that the emission could not come from an inner hot accretion disc. The X-ray emission for NGC 4261 could be better constrained if we allow an accretion rate of $\dot{M}_{0}=3.2\times 10^{-2}\dot{M}_{\rm Edd}$.

For NGC 4261 and NGC 315, the fitting was not able to fully explain the UV emission. Both galaxies presented much lower observational value in UV than our RIAF model prediction.

Astronomers have studied M87 extensively over the past sixty years (Felten, 1968; Schreier et al., 1982; Biretta et al., 1999; Doeleman et al., 2012; Event Horizon Telescope Collaboration et al., 2019). Recently, Bandyopadhyay et al. (2019) modeled this source with a combination of RIAF and synchrotron jet. They used a population of nonthermal electrons, which is absent in our work, and they manage to fit the bump at $\sim 100$GHz combining RIAF and jet emission. If they considered only thermal electron RIAF, the disc emission is not strong enough compared to the jet emission.

For NGC 4261 and NGC 315, de Menezes et al. (2020) did the previous SED modeling for the same data set with our original code with iterative procedures. In this work, we reported an improved version of these fits with AGNNES. Our results suggested that $\dot{M}_{0}\gtrsim 10^{-2}\dot{M}_{\rm Edd}$, which is not consistent with the RIAF scenario – notwithstanding our model fits well the data. de Menezes et al. (2020) Gamma-ray observations indicate the emission comes from a jet emitting synchrotron (low-frequency) and synchrotron self-Compton (high-frequency).

The original model used to calculate the SEDs has its limitations, and AGNNES inherited all of them. The RIAF SED followed a straightforward model, with some strong approximations. First of all, it is a 1-dimensional calculation. It was assumed an axisymmetric accretion flow, and we integrated over the $z$-component (considering a cylindrical geometry). Furthermore, the code works on stationary accretion flow, hence it cannot account time variability.

Our jet SED took account of synchrotron emission only. We did not calculate the Inverse Compton effect in the jet photons. Some observations support a scenario with synchrotron self-Compton emission in the jet dominating the LLAGN emission (Nagar et al., 2005; Finke et al., 2008; Takami, 2011).

The SED data of every object was not simultaneously measured. Many authors observed each data point at a different time, with years of separation. Moreover, we did not take into account any variability. We assumed all data comes from the quiescent state of the source. Such assumptions can be a problem for galaxies with few data points and not known variability patterns. M87 has more reliable data because it is a source continuously observed.

AGNNES works properly inside its training parameters range (see section 2.1). For fainter sources with very low luminosity, like SgrA*, AGNNES can not fit the SED reliably. For the SgrA* case, the accretion is too small to be inside our training specifications. The fiducial code can calculate SEDs for objects with such low accretion rates. However, numerical errors increase considerably as the accretion rate decreases. To generate data for low accretion rates are harder than for mildly higher accretion rates. The difficulty in creating SEDs with $\dot{M}_{0}\gtrsim 10^{-5}$ impacted on AGNNES results. Below this $\dot{M}_{0}$ value, our NN can not calculate a correct SED. AGNNES best reproduces LLAGN SEDs with $(\nu L_{\nu})_{\rm peak}>10^{36}$ erg/s.