Cosmological evolution of quasar radio emission in the view of multifractality

TL;DR

This study investigates the multifractal behavior of quasar radio emissions across different frequencies and frames, revealing intrinsic multifractality and differences in dynamics, with implications for understanding their cosmological evolution.

Contribution

It applies wavelet transform multifractal analysis and ARIMA modeling to quasar radio emissions, demonstrating intrinsic multifractality and frequency-dependent dynamics.

Findings

Strong multifractality observed in all sources.

Similar multifractality at different frequencies suggests common radiation mechanisms.

ARIMA models fit some quasars but are inadequate for others.

Abstract

Variations in scaling behavior in the flux and emissions of distant astronomical sources with respect to their cosmic time are important l phenomena that can provide valuable information about the dynamics within the sources and their cosmological evolution with time. Different studies have been applying linear analysis to understand and model quasars' light curves. Here, we study the multifractal behavior of selected quasars' radio emissions in their observed frame (at 22 and 37 GHz bands) and and their rest frame. To this end, we apply the wavelet transform-based multifractal analysis formalism called wavelet transform modulus maxima. In addition, we verify whether the autoregressive integrated moving average (ARIMA) models fit our data or not. In our work, we observe strong multifractal behavior for all the sources. Additionally, we find that the degree of multifractality is strongly similar for each source and significantly different between sources at 22 and 37 GHz. This similarity implies that the two frequencies have the same radiation region and mechanism, whereas the difference indicates that the sources have intrinsically different dynamics. Furthermore, we show that the degree of multifractality is the same in the observed and rest frames of the quasars, i.e., multifractality is an intrinsic property of radio quasars. Finally, we show that the ARIMA models fit the 3C 345 quasar at 22 GHz and partially fit most of the time series with the exception of the 3C 273 and 3C 279 quasars at 37 GHz, for which the models are found to be inadequate.

Tables

Table 1. Table 1: The time range of observation and number of data points for each source at 22 and 37 GHz.

Source name		Frequency	Time range	Number of data
				points
		22 GHz	1980-2004	938
3C 273		37 GHz	1979-2018	1954
		22 GHz	1980-2004	757
3C 279		37 GHz	1979-2018	1961
		22 GHz	1980-2004	805
3C 345		37 GHz	1980-2018	1660
		22 GHz	1980-2004	760
3C 454.3		37 GHz	1980-2018	2292

Table 2. Table 2: The calculated multifractal spectrum width ( Δ ​ α Δ 𝛼 \Delta\alpha ) for each source at 22 and 37 GHz.

Source name	Red-shift (z)	Observation frame frequency	$\Delta$ $\alpha$
		22 GHz	1.5961
3C 273	0.158	37 GHz	1.5854
		22 GHz	0.9763
3C 279	0.536	37 GHz	0.9653
		22 GHz	0.9745
3C 345	0.595	37 GHz	0.9633
		22 GHz	0.8951
3C 454.3	0.859	37 GHz	0.8432

Equations

$$W(s,a)=\frac{1}{\sqrt{s}}{\int }_0^T\left(\Psi \left(\frac{t-a}{s}\right).X(t)\right)dt$$

$${\int }_{-\infty }^{+\infty }{t}^m\Psi (t)dt=0,\forall m,0\le m\le n$$

$${\Psi }^{(N)}(t)=\frac{d^N}{d{t}^N}{e}^{-\frac{t^2}{2}},$$

$$\Psi (t)=(1-t^2).e^{-\frac{t^2}{2}}$$

$$W_{s,a}=W(s,a)|(s,a\in N)\wedge (s\in [1,s_{max}])\wedge (a\in [1,T]),$$

$$W_{s,a}^{sq}=(W(s,a){)}^2|(s,a\in N)\wedge (s\in [1,s_{max}])\wedge (a\in [1,T]),$$

$$\begin{gathered} LcM{x}_{s,a}=\text{cases}1|\frac{\partial (W(s,a){)}^2}{\partial a}=0, \\ 0|\neg (\frac{\partial (W(s,a){)}^2}{\partial a}=0) \end{gathered}$$

$$LcM{x}_{s,a}=LcMx(s,a)$$

$$\begin{gathered} LcM{x}_{s,a}=\text{cases}1|(\frac{\partial (W(s,a){)}^2}{\partial a}=0)\wedge (s<a<T-s), \\ 0|\neg (\frac{\partial (W(s,a){)}^2}{\partial a}=0)\vee (T-s\le a)\vee (a\le s) \end{gathered}$$

$$Zq(s)=\sum _{a=1}^{T-1}(C(s).WTMM{)}^q|(LcM{x}_{s,a}=1),$$

$$Zq(s)\sim s^{\tau (q)},$$

$$\tau (q)=\underset{s\to 0}{\lim }\frac{\ln (Zq(s))}{ln(a)},$$

$$\alpha =\alpha (q)=\frac{\partial \tau (q)}{\partial q},$$

$$f(\alpha )=q.\alpha -\tau (q),$$

$$y_t^{\prime }=c+{\varphi }_1{y}_{t-1}^{\prime }+...+{\varphi }_p{y}_{t-p}^{\prime }+{\theta }_1{\epsilon }_{t-1}+...+{\theta }_q{\epsilon }_{t-q}+{\epsilon }_t,$$

Full text

methods: statistical – galaxies: active – (galaxies:) quasar: individual: 3C 273, 3C 279, 3C 345, 3C 454.3

1 Introduction

Radio astronomy remains one of the most significant research areas in astronomy and cosmology to date; this is because a wide variety of astronomical objects that are not detectable in the optical wavelengths as well as thermal and nonthermal radiation mechanisms and propagation phenomena can be studied at radio wavelengths. Radio observations have played a significant role in the discovery, subsequent observation, and classification of active galactic nuclei (hereinafter AGNs), radio galaxies and other radio sources (e.g., Singh (2016); Kharb (2016); Tamhane (2015); Saikia & Jamrozy (2009); Brocksopp (2007); Hota & Saikia (2006); Reynolds & Begelman (1997); Greenstein & Schmidt (1964)). Moreover, as an application of radio astronomy, studying the behavior of the radio luminosity function (e.g., Willott (2001); Dunlop & Peacock (1990)) plays a significant role in understanding the formation and evolution of radio galaxies and the cosmological evolution of radio sources, which in turn help to uncover the physics of the early universe.

Quasars, an extremely luminous compact region thought to reside at the center of galaxies, are distant astronomical objects that belong to a subclass of AGNs (e.g., Padovani et al. (2017); Patiño-Álvarez et al. (2017); Fernandes et al. (2017); Courvoisier (1987)). Flat spectrum radio quasars (hereinafter FSRQs), which include both high- and low-polarization quasars (hereinafter HPQs and LPQs, respectively), along with BL Lac objects, form a subset of AGNs known as blazars. The former are known for their high luminosity and strong and broad emission features in their spectra, whereas the latter are known for their low luminosity and very weak or even absent emission features in their spectra (Marchã, 2001). Radio-loud quasars/blazars show some of the most violent high-energy phenomena observed in AGNs to date. They emit radiation across the electromagnetic spectrum, from the radio to X-rays and/or gamma rays (e.g., Padovani et al. (2017); Patiño-Álvarez et al. (2017); Fernandes et al. (2017); Villata et al. (2009); Courvoisier (1987)) and have the characteristics of very high luminosity, nonthermal radiation, strong radio emission, and large flux fluctuations and have been known as highly variable energetic sources across the electromagnetic spectrum though most (if not all) of them are strongly variable in a very short time scale at higher frequencies (Lawson, 1999; Courvoisier, 1988; Wehrle, 1998; Ulrich, 1997; Netzer et al., 1996; Dondi & Ghisellini, 1995; Cotton, 1979); moreover, high and variable polarization and superluminal motion are characteristics of blazars. The low-energy blazar emission is thought to be the result of electron synchrotron radiation, with the peak frequency reflecting the maximum energy at which electrons can be accelerated (Petropoulou & Dimitrakoudis, 2015).

The timescale variability and other properties of AGNs in general, and of quasars in particular, have been extensively studied at several radio frequencies using different time series analyses approaches. Here, we study the multifractal behavior of selected radio-loud quasars such as 3C 273, 3C 279, 3C 345, and 3C 454.3. The flat spectrum radio-loud quasar 3C 273 is a low polarization quasar (León-Tavares, 2011; Valtaoja, 1991) at a distance of z = 0.158 (Arp, 1995). This radio quasar has been monitored by Mets$\ddot{a}$hovi Radio Observatory at 22 and 37 GHz for 24 and 39 years, respectively, and it is among AGNs whose emission has been studied at all wavelengths (Courvoisier et al., 1990; Courvoisier, 1987). Particularly, the radio variability of 3C 273 has been studied by (Türler, 2000, 1999). The HPQ quasar 3C 279 (León-Tavares, 2011) at z = 0.536 (Netzer et al., 1996) has been observed by Mets$\ddot{a}$hovi Radio Observatory at 22 and 37 GHz for 24 and 39 years, respectively. The flat spectrum radio quasar 3C 279 has been known as one of the brightest quasars at all wavelengths, and its multiwavelength behavior and jet structure have been studied by different authors (e.g. Zheng & Yang (2016); Chatterjee (2008); Hartman et al. (2001)). The HPQ quasar 3C 345 (León-Tavares, 2011) at a distance of z = 0.595 (Netzer et al., 1996; Moore & Stockman, 1984; Burbidge, 1965) is one of the most studied radio-loud flat-spectrum quasars whose radio spectrum is dominated by a flat-spectrum core (Kollgaard, 1989). 3C 345 has been monitored by Mets$\ddot{a}$hovi Radio Observatory at 22 and 37 GHz for 24 and 38 years, respectively. The radio flux observations of 3C 345 have shown that variations are common in the radio frequency range (Aller, 1985). Vio (1991) indicated that the light-curve of 3C 345 is nonlinear and stochastic. The quasar 3C 454.3 is an HPQ (León-Tavares, 2011) at z = 0.859 (Netzer et al., 1996) with a flat spectrum that belongs to the blazar class of active galactic nuclei (Villata et al., 2006), which has been studied across the electromagnetic spectrum (Jorstad et al., 2010) and known as variable source at all wavelengths. This quasar has been monitored by Mets$\ddot{a}$hovi Radio Observatory at 22 and 37 GHz for 24 and 38 years, respectively.

Using the structure function, the discrete correlation, and Lomb-Scargle periodogram analysis approach, Hovatta (2007) has studied the variability time scales of a sample of AGN including our candidates at several frequencies between 4.8 and 230 GHz, has shown that the slopes between frequencies were different and did not find a significant difference between most of the sources. Similarly, Valtaoja (1992b) studied variability timescales of a sample of AGNs including our sources at 22 and 37 GHz and found that the spectral and variability characteristics between most of the sources were very similar, except for the difference observed between HPQs and LPQs. Additionally, Hovatta (2008) has studied the variability behavior of 80 AGNs including our candidates at 22, 37 and 90 GHz using wavelet analysis and has shown that the timescales at 22 and 37 GHz did not differ significantly from source to source. Using the structure function analysis, Wang & Yang (2010) studied the variability characteristics of 3C 273 and 3C 345 at 22 and 37 GHz and found that each source has similar slopes at these two frequencies, indicating that for both sources, the two frequencies have the same region and mechanism. Moreover, Valtaoja (1992a) studied the flux variations of extragalactic radio sources at 37, 22, 14.5, 8, and 4.8 GHz and suggested that all outbursts in AGN have similar evolutionary tracks, defined by the motion of the turnover peak of the shock spectrum in time, and most of the observed differences result from variations of the frequency at which the outburst reaches its maximum development. Similarly, Valtaoja (1999) studied the total flux density variations of a group of AGNs at 22 and 37 GHz and has shown that the variations were adequately decomposed into flares having an exponential rise.

See Zhang (2018); Berton (2018); Gupta et al. (2017); Jorstad et al. (2013); Bachev (2011); Dong (2010); Schinzel (2010); Villata et al. (2009); Chatterjee (2008); Raiteri et al. (2008); Hovatta (2007); Lindfors (2006, 2005); Lawson (1999); Valtaoja (1999); Wehrle (1998); Stevens (1994); Valtaoja (1992a, 1991); Courvoisier et al. (1990); Courvoisier (1987)) for more studies on the variability timescale and other physical properties of our candidates across the electromagnetic spectrum.

Most (if not all) of the approaches used in the abovementioned works were linear analyses. Here, we apply a nonlinear analysis approach to reveal the nonlinear characteristics about the sources considered, and this is the main motivation for investigating the new methods and emphasizes the importance of the present work. In this work, we apply a wavelet transform-based multifractality analysis approach called wavelet transform modulus maxima (hereinafter WTMM) and study the multiscaling or multifractal behavior of the quasars 3C 273, 3C 279, 3C 345, and 3C 454.3 radio observations at 22 and 37 GHz. Why multifractality analysis? Most astrophysical objects are possibly associated with continuous nonlinear stochastic systems due to their complexity in nature. A fractal behavior can be observed in the time series of complex systems (Maruyama, 2016). It has been shown that quasars, in general, are among complex systems that have nonlinear time series characterized by fractal behavior (Vio, 1991) and also by sudden bursts of very large amplitude (Barbieri, 1990; Kidger, 1989), which implies that the dynamical evolution of quasars is nonlinear (i.e., described by nonlinear stochastic differential equations) (Vio, 1991). Additionally, there is a suggestion that extra-galactic radio sources are intermittent on timescales of $\sim$ $10^{4}-10^{5}$ yr (Reynolds & Begelman, 1997). Traditionally, different classical approaches, such as power spectrum function (PS), structure function (SF), PS-periodogram, covariance analyses, and others, which are suitable only for addressing the signals characteristic of linear systems (Vio, 1992), have been in use to analyze and study nonlinear signals due to the unavailability of better approaches that are necessary to gain detailed information about the dynamics of complex systems. The multifractal formalism was introduced in the mid-1980s to provide a statistical description of the fluctuations of regularity of singular measures found in chaotic dynamical systems (Halsey, 1986). Currently, multifractal analysis is being used in several fields of science to characterize nonlinearity or detect singularity in various types of signals from complex systems and study correlations between different physical parameters (Trevino & Dal Negro, 2012). Currently, the multifractal analysis approach has been applied to a large number of empirical as well as theoretical studies of a wide variety of problems. As an example, using multifractality analysis, it has been shown that the X-ray light curve of the BL Lac object PKS 2155-304 is monofractal, and the optical light curve of the quasar 3C 345 has a multifractal nature - nonlinear behavior (Vio, 1992). In addition, using multifractality analysis, it has been shown that the distribution of The Infrared Astronomical Satellite (IRAS) galaxies is homogeneous at the large scale (Pan & Coles, 2000), which is in agreement with the cosmological principle. Additionally, Bewketu Belete et al (2018) has studied the fractal nature of the light curves of 3C 273 at specific frequencies across its electromagnetic spectrum and shown that most of the light curves have presented a multifractal signature, confirming the nonlinear and intermittent nature of the source. For more studies on multifractality analysis by different authors in different science cases, see Maruyama (2017, 2016); Kasde (2016); Agarwal (2016); de Freitas (2016); Aliouane & Ouadfeul (2013); Jagtap (2012); Ouahabi & Femmam (2011); Nurujjaman (2009); Lin & Sharif (2007); Degaudenzi & Arizmendi (1998).

The aims of our work are as follows: (i) to study the multifractal behavior (if any) of selected radio-loud quasars located at different redshifts, (ii) to verify how fractal signatures of each source, and between sources, behave at 22 and 37 GHz (any similarity or difference? If yes/no, what can we learn about them?), and (iii) to determine any possible correlation between multifractal behavior and redshift at those two frequencies, i.e., how the fractal signature of these radiations change with redshift (the cosmological evolution of quasars’ radio emissions from the view of multifractality analysis). The first aim helps to draw a conclusion regarding whether radio-loud quasars are multifractal/intermittent and nonlinear systems. From the second aim, we can roughly understand whether the selected radio frequencies have the same radiation region and mechanism or not. Furthermore, knowledge of the correlation between multifractal behavior and redshift is essential not only for robustly understanding the cosmological evolution of quasars’ multifractal signature at specific radio frequencies and for interpretation of the behavior of the relativistic plasma and nonthermal radiation associated with the jet out flowing from the black hole/accretion disk systems but also for supporting the claim that quasars’ redshift is cosmological in nature. In all our discussions, a flat $\Lambda$CDM cosmology with $\Omega_{\lambda}$ = 0.70, $\Omega_{m}$ = 0.3, and $H_{o}=70kms^{-1}$ Mp$c^{-1}$ is used, unless otherwise specified. Our work is structured as follows: in section 2, we present the light curves of the data used, the method and procedures. The results obtained and a discussion are presented in section 3, and the summary and conclusions are included in section 4.

2 Data Collection, Method and Procedures

2.1 Data Collection: Light-curves

We have collected the radio flux data of the sources 3C 273, 3C 279, 3C 345, and 3C 454.3, at 22 and 37 GHz frequencies from Aalto University Mets$\ddot{a}$hovi Radio Observatory in Finland. We have full light curves for 37 GHz and only up to 2004 for 22 GHz. Some of the most recent data have not been published yet, and the light curves up to 2004 are published in Teräsranta (2005). Our selection of these two frequencies takes advantage of the very long history of radio flux observations of those quasars at these bands. The light curves of the sources in the observation and rest frames are given in Figs. 1 and 2, respectively. We have corrected our data as $f_{rest}=f_{obs}*(1+z)$ , $t_{rest}=t_{obs}/(1+z)$, and $S_{rest}=S_{obs}/(1+Z)^{(1+\alpha)}$, where $f$ is frequency, $t$ is time, $S$ is flux, $z$ is redshift, and $\alpha$ is the spectral index. In the radio domain, we assume that $\alpha=0$ as the spectrum is more or less flat.

For all our candidates, we have an unequal length of data streams at both frequencies (Table 1). In WTMM-based multifractaltiy analysis, though the choice of the scale parameter is dependent on the length of a time series, it mainly depends on the absence and presence of local maxima lines at the scale considered, i.e., whether the scale at which the calculated wavelet coefficients hold maxima lines or not. At different scales, we can have different numbers of local maxima lines that carry local information about any singularity contained in that part of the time series. In general, the data length determines the choice of scale parameter. It is one of the parameters we use to calculate wavelet coefficients, from which we calculate the local maxima lines that we use in the multifractality analysis part. Obviously, the scale parameters at which we obtain wavelet coefficients that hold maxima lines for the current light curves change when the length of the time series changes, which in turn affects the multifractality strength to be calculated. For our case, we have chosen the most informative scale parameter, i.e., the scale parameter at which we have better maxima lines, at each frequency for all the time series, and therefore, our discussion of results takes this into consideration.

2.2 Method and Procedures

Signals can be efficiently represented by decomposing them in different frequencies using the Fourier analysis method. However, the most valuable information in a complex system signal is contained by its irregular structures and transient phenomena called singularities, and particularly in physics, it is important to analyze irregular structures in a signal to deduce properties about the underlying physical phenomena (Arneodo, 1988; Mandelbrot & Whitrow, 1983), which is beyond the capability of Fourier analysis because it only decomposes a signal into its frequency domain. Therefore, the Fourier transform is not powerful and preferable for multifractality analysis, which requires a special technique of decomposing a signal into time and frequency domains. The continuous wavelet transform is an excellent tool for mapping the changing properties of nonstationary signals. Because of its capability of decomposing a signal into small fractions that are well localized in time and frequency and of detecting local regularities of a signal (areas on the signal where a particular derivative is not continuous) such as nonstationarity, oscillatory behavior, breakdown, discontinuity in higher derivatives, the presence of long-range dependence, and other trends, wavelet analysis remains one of the most preferable signal analysis techniques to date (Maruyama, 2016; Andrejs Puckovs & Andrejs Matvejevs, 2012). These strengths of wavelet transform make it preferable to other traditional singularity analysis techniques, and there is a claim that it is suitable for multifractal analysis and allows for reliable multifractal analysis to be performed Muzy & Arneodo (1991). Therefore, due to these and other reasons not mentioned here, we have chosen to apply a multifractality analysis approach that requires the continuous wavelet transform known as the wavelet transform modulus maxima. The WTMM was originally introduced by (Muzy & Arneodo, 1991). Basically, WTMM-based multifractality analysis consists of two statistically connected parts: the wavelet transform part and multifractality analysis part. Each part is discussed below.

**I. Wavelet transform formalism

**A. Continuous wavelet transform:

The direct continuous wavelet transform of a given signal $X$($t$) can be represented by:

[TABLE]

where $W$ are the wavelet coefficients; $\Psi$($s$,$a$,$t$) - the mother wavelet function; $s$ - the scaling parameter; $a$ - the shift parameter; $X$ - the signal; $t$ - the time at which the signal is recorded; and $T$ - maximal time value or signal length. The analyzing wavelet $\Psi$($t$) is generally chosen to be well localized in space and frequency. Usually, $\Psi$($t$) is only required to be of zero mean, but in addition to these requirements, for the particular purpose of multifractal analysis, $\Psi$($t$) is also required to be orthogonal to some low-order polynomials, up to the degree $n$-1 (i.e., to have $n$ vanishing moments) Enescu (2006):

[TABLE]

A class of commonly used real-valued analyzing wavelets, which satisfies the condition given by eq. 2, is given by the successive derivatives of the Gaussian function Enescu (2006).

[TABLE]

for which $n$ = $N$. Our analyzing wavelet is the Mexican Hat (MHAT) wavelet (second-order Gaussian wavelet), which is one of the wavelets that has been applied for WTMM-based analysis, and represented by the relation:

[TABLE]

where $\Psi$($t$) - the mother wavelet function; $t$ - the time at which the signal is recorded. At lower scales $s\approx 0$, the number of local maxima lines (hereinafter LcMx) tends to infinity. Since it is the maxima line that points toward each regularity or carries information about any singularity or nonlinearity in a signal (Muzy, 1994; Mallat & Hwang, 1992; Mallat & Zhong, 1992), it is unnecessary to calculate wavelet coefficients that do not contain maxima line(s). Though there is a suggestion that the scaling parameter $s$ used in the WTMM approach is limited to $s\leq$ [128], it should be in the interval [1, $\frac{T}{2}]$ and can also be in the interval [1, $\frac{T}{4}$], which is still informative, mainly to reduce computation time (Andrejs Puckovs & Andrejs Matvejevs, 2012). The shifting parameter $a$ cannot be greater than the signal length $T$, and therefore, $a$ $\leq$ T.

The calculated wavelet coefficients $W_{s},{a}$ can be written in a matrix form given by (Andrejs Puckovs & Andrejs Matvejevs, 2012):

[TABLE]

where $W_{s,a}$ - the wavelet coefficients; $s_{max}$- the maximal scaling parameter; $s$ - the scaling parameter; $a$ - the shifting parameter; and $T$ - the signal length. Additionally, one can calculate the absolute wavelet coefficients in matrix form as:

[TABLE]

where $W^{sq}$ - the squared wavelet coefficients matrix. Other parameters are as explained above. In the wavelet transform output plot, wavelet coefficients are colored by their absolute values.

B. Skeleton function construction: The skeleton function is nothing but a collection of maxima lines at each scale of the calculated wavelet coefficient matrix, i.e., it is a scope of all local maxima lines that exist on each scale $s$. In other words, the skeleton matrix construction is a technique of excluding coefficients in the absolute wavelet coefficients matrix that are not maximal. As a result, in the skeleton matrix, only absolute wavelet coefficients that belong to local maxima lines exist. The need to collect all the maxima lines at each scale together in matrix form, the skeleton function, is from the fact that it is the maxima lines that carry valuable information about the signals, i.e, maxima lines point toward regularity in the signal. We construct the skeleton function as follows (Andrejs Puckovs & Andrejs Matvejevs, 2012):

[TABLE]

under the conditions ($s$,$a\in N$) $\wedge$ ($s\in[1,s_{max}$]) $\wedge$ ($a\in[1,T]$), where $LcMx$ - the wavelet skeleton function; $W$($s$, $a$) - the wavelet coefficeints; $s$ - the scaling parameter; $s_{max}$ - the maximal scaling parameter; $a$ - the shift parameter; and $T$ - the signal length. The wavelet skeleton function can also be calculated from squared wavelet coefficients matrix and expressed in a matrix form as follows:

[TABLE]

under the same conditions as eq. 7. Since wavelet coefficients on corners provide little or no information, and consequently, local maxima lines (LML) on corners also provide no significant information about the singularity in the signal, we therefore take the edge effect into consideration by removing the local maxima lines on corners using the following formula (Andrejs Puckovs & Andrejs Matvejevs, 2012):

[TABLE]

the skeleton function is a logical function that has only two variables, 0 and 1; one is used if the skeleton matrix element is a local maximum, and zero is used otherwise.

There is one problem in WTMM-based multifractality analysis: the skeleton function does not contain local maxima lines as required for fractal analysis but has disconnected broken lines, gaps, and single points in the LcMx matrix. Therefore, it is mandatory to apply some technique to fix these limitations. To that end, we applied an algorithm called the supremum algorithm, which consists of the following steps.(Andrejs Puckovs & Andrejs Matvejevs, 2012):

1. Define matches (relations between single local maxima points);
2. Define match conflicting (one cell to more) and nonconflicting cases;
3. Create chains from pairs;
4. Chain interpolation (to fill missing wavelet coefficients on WTMM line);
5. Add points to LcMc map (on line gaps);
6. Add single points to LcMc map; and
7. Change variables (at the end).

**II. Multifractality analysis part

**C. Fractal partition function estimation: At this time, the multifractality analysis starts. Using the wavelet modulus maxima coefficients obtained in the wavelet transform part, we calculate the thermodynamic partition function, a function that connects the wavelet transform and multifractality analysis part. Muzy & Arneodo (1991) defined the thermodynamic partition function $Zq$($s$) as the sum of the $q$-th powers of the local maxima of the modulus to avoid division by zero:

[TABLE]

where $Zq(s)$ - the thermodynamic partition function; $WTMM$ - the wavelet modulus maxima coefficients; $C$($s$) - the constant depending on the scaling parameter $s$; $s$ - the scaling parameter; $q$ - the moment, which takes any interval with zero mean, for our case $q$ $\in$ [-5, 5]; and $LcMx$ - the wavelet skeleton function (aggregate of local maxima lines in matrix form). The thermodynamic partition function is a function of two arguments - the scaling parameter $s$ and the power argument $q$. The moment $q$ discovers different regions of the singularity measurement in the signal, i.e., it indicates the presence of wavelet modulus maxima coefficients of different values. The condition $LcMx_{s,a}$ = 1 is to inform that only modulus maxima coefficients are used. The thermodynamic partition function is finite if the wavelet modulus maxima coefficients are not equal to zero ($WTMM_{s,a}\neq 0$). To satisfy the condition, all zero coefficients should be neglected in the wavelet modulus maxima matrix WTMM. The origin of zero coefficients in wavelet modulus maxima matrix is in the LcMx wavelet skeleton function - all elements in the skeleton matrix that are not local maxima are zero-valued elements. What is important here is the relationship between the $Zq$($s$) and $s$, which determines the scalability of the signal under consideration. We investigate the changes of $Zq$($s$) in the time series at a different scales $s$ for each $q$. A plot of the logarithm of $Zq$($s$) against the logarithm of the time scale $s$ was created. This plot shows how $Zq$($s$) scales with $s$ and reveals the strength and nature of local fluctuations at each scale in the signal. In the WTMM approach, the wavelet transform maxima are used to define a partition function whose power-law behavior is used for an estimation of the local exponents. At small scales, the following relation is expected:

[TABLE]

where $\tau(q)$ - the scaling exponent function, which is the slope of the linear fitted line on the log-log plot of $Zq$($s$) and $s$ for each $q$.

D. The scaling exponent function $\tau(q)$: This is a function of one argument $q$ that is determined from the slope of the linear fitted line on the log-log plot of $Zq$($s$) against the logarithm of the time scale $s$ for each $q$, which means that the behavior of the scaling function $\tau$($q$) is completely dependent on the nature of the thermodynamic partition function, or in other words, the behavior of $\tau$($q$) is dependent on the scaling relationship between $Zq$($s$) and $s$. The mathematical relation of calculating $\tau$($q$) is given by the relation:

[TABLE]

where $Zq$($s$)- the thermodynamic partition function; $\tau$ - the local scaling exponent; $s$ - the scaling parameter; and $q$ - the moment. The condition $\tau$($q$ = 0) $+$ 1 = 0 is important for the multifractal spectrum calculation(Andrejs Puckovs & Andrejs Matvejevs, 2012). We define monofractal and multifractal as follows: the time series is said to be monofractal if $\tau$($q$) is linear with respect to $q$, and if $\tau$($q$) is nonlinear with respect to $q$, then the time series considered is classified as multifractal (Frish and Parisi, 1985).

E. The multifractal spectrum function, $f$($\alpha$): Once we determine the scaling function $\tau$($q$), it is necessary to estimate the multifractal spectrum $f$($\alpha$) to be able to fully draw conclusions about the multifractal behavior of the signal considered. Using the calculated scaling function, we estimate the multifractal function via Legendre transformation as (Halsey, 1986):

[TABLE]

where $\alpha$ is the singularity exponent or Holder exponent.

[TABLE]

where $f$($\alpha$) is the multifractal spectrum function. We extract two important items of information from $f$($\alpha$) against the $\alpha$ plot: the width ($\Delta\alpha$ = $\alpha_{max}$ - $\alpha_{min}$) and the symmetry in the shape of $\alpha$ defined as $A=(\alpha_{max}-\alpha_{0})/(\alpha_{0}-\alpha_{min})$ , where $\alpha_{0}$ is the value of $\alpha$ when $f(\alpha)$ assumes its maximum value. Ashkenazy (2003) and Shimizu (2002) have proposed that the width of a multifractal spectrum is the measure of the degree of multifractality. Smaller values of $\Delta\alpha$ (i.e., $\Delta\alpha$ becomes close to zero) indicate the monofractal limit, whereas larger values indicate the strength of the multifractal behavior in the signal Telesca (2004). For the symmetry in the shape of $\alpha$, the asymmetry presents three shapes: asymmetry to the right-truncated ($A>1$), left-truncated ($0<A<1$), or symmetric ($A=1$). Ihlen (2012) presented that the symmetric spectrum originated from the leveling of the $q^{th}$-order generalized Hurst exponent for both positive and negative $q$ values. The leveling of $q^{th}$-order Hurst exponent reflects that the $q^{th}$-order fluctuation is insensitive to the magnitude of the local fluctuation. When the multifractal structure is sensitive to the small-scale fluctuation with large magnitudes, the spectrum will be found with a right truncation, whereas the multifractal spectrum will be found with left-side truncation when the time series has a multifractal structure that is sensitive to the local fluctuations with small magnitudes. Therefore, the width and shape of a multifractal spectrum are able to classify small and large magnitude (intermittency) fluctuations and determine the degree of the multifractality signature in a given signal.

2.2.1 Autoregressive Integrated Moving Average (ARIMA) Model

Quasars are known to be extremely variable in a short time scale from hours to months due to the rapid change in their accretion rate (e.g., Wagner (1995); Quirrenbach (1993)). Thus, a time series model that is suitable for a nonstationary time series with short memory may be preferable to model a quasar time series. The autoregressive integrated moving average (hereinafter ARIMA) model is the extension of autoregressive moving average model (ARMA), which in turn is a combination of autoregressive (AR) and moving average (MA) linear time-series models. The AR, MA, and ARMA models are only applied for stationary time series, whereas ARIMA includes the case of nonstationarity as well. In addition, ARIMA models are capable of describing short-memory autocorrelations that exist in a time series. The ARIMA model is one of the autoregressive time-series models rarely used to understand the nature of astronomical times series (Feigelson, 2018). Here, we apply the ARIMA model for the light curves in Fig. 1. The aim is only to check whether ARIMA fits our data or not, from which we can learn the type of correlations that exist in quasars’ time series and their behavior. Since seasonality is not common in quasar time series but trends are, we use the nonseasonal ARIMA(p,d,q) model discussed in Hyndman and Athanasopoulos (2014):

[TABLE]

where $y^{\prime}_{t}$ is the differenced time series; p is the order of the autoregressive part- AR(p); q is the order of the moving average part- MA(q); and $\varepsilon_{t}$ is white noise. The ARIMA(p,d,q) model contains two predictors: a linear combination of lagged values of the variable- AR(p) and the linear combination of lagged errors- MA(q) model. The reference used for this model can be used to understand the details.

In fitting the model, we follow the procedures discussed in Hyndman and Athanasopoulos (2014): (1) We check for stationarity - if nonstationary, we apply differencing. Usually the autocorrelation function (hereinafter ACF) of a time series decays exponentially to zero if it is stationary and slowly to zero if it is nonstationary; (2) we plot the ACF and partial autocorrelation function (hereinafter PACF) of the already differenced time series; (3) we first estimate the parameters $p$ and $q$ based on ACF and PACF in (2); (4) we fit the model ARIMA(p,d,q) using the $p$ and $q$ estimated in (3); and (5) we calculate and plot the corresponding residuals. We accept the model only if its residual is white noise. If the residual of the fitted model is not white noise, we repeat the steps from 3 to 5. In estimating the parameters $p$ and $q$, we use the function auto.arima() that provides us with a better model, though not always.

3 Results and Discussion

In our work, we apply WTMM based-multifractality analysis and study the multifractal or nonlinear behavior of radio emissions of selected radio-loud quasars 3C 273, 3C 279, 3C 345, and 3C 454.3 in the observation frame (at 22 and 37GHz ) and in the rest frames ($f_{rest}=f_{obs}*(1+z)$). The aim is to search for the presence of a multifractal signature in the light curves of the sources at each band and verify whether there is any similarity or difference in the degree of multifractaltiy between the bands and between the sources both in their observed and rest frames.

3.1 Analysis of light curves in the observation frame:

Here, we analyze the multifractal behavior of the light curves in the observed frame at 22 and 37 GHz, Fig. 1. The scenario emerging for each light curve is discussed separately as follows. In Figs. 3 and 6, we present the continuous wavelet transform maps and the constructed skeleton functions for all sources at 22 and 37 GHz, respectively. The wavelet maps present the calculated wavelet coefficients in their absolute values and colored as dark and light in accordance with the color map given. The dark and light colors show lower and higher absolute wavelet coefficients, respectively. From the wavelet coefficient matrix, we select local maxima lines, or maxima points at each scale. The collection of all local maxima lines, or maxima points, at each scale forms a function called the skeleton function. The need to construct the skeleton function is to consider only local maxima lines at each scale and simplify the multifractality analysis. In Figs. 4 and 10, we present the 3D plots of the thermodynamic partition functions, a function of two arguments - scaling parameter $s$ and moment $q$ - at 22 and 37 GHz for all the light curves to show how the thermodynamic function $Zq(s)$, the moment $q$, and the scaling parameter $s$ behave together. It is the thermodynamic partition function that connects the wavelet formalism to the multifractal formalism. To clearly see the scaling behavior between $Zq$($s$) and $s$, the scalability of the signals, we have created the log-log plots of $Zq$($s$) and $s$ using eq. 12 for all the light curves and represented the calculated slope by the scaling exponent function $\tau$($q$) plots as shown in Figs. 5 and 8. The observed nonlinearity between the scaling exponent $\tau(q)$ vesus the moment $q$ at both bands reveals the presence of nonlinear scaling behavior between the thermodynamic partition function $Zq$($s$) and the scale $s$ for all the light curves. Furthermore, the observed nonlinearity between $\tau$($q$) and $q$, which is the slope of log($Zq$(s)) against log($s$) plots, clearly indicates the presence of multifractal behavior in all the light curves though the degree of nonlinearity that varies between sources at both bands.

Based on the scaling exponent function $\tau$($q$), Figs. 5 and 8, we can see the similarity and difference in the degree of nonlinearity at the two bands for each source, Fig. 9, and between the sources, Figs. 5 and 8. After observing the presence of the multifractal signature based on the nonlinear scaling exponent function, the degree of multifractality is determined (i.e., how strong is the observed multifractal signature) at each band for all the sources. To this end, we estimate the multifractal spectrum function $f$($\alpha$) and calculate the width ($\Delta\alpha$ = $\alpha_{max}$ - $\alpha_{min}$) using eqs. 13 and 14. The broader the spectrum (the larger the value of $\Delta\alpha$), the richer the multifractality is. The width value tells us how strong is the observed multifractal signature, which is an additional parameter one can use to see the difference in the degree of nonlinearity, or multifractality between signals. The calculated width $\Delta\alpha$ values for 3C 273, 3C 279, 3C 345, and 3C 454.3 at 22 and 37 GHz are 1.5961/0.9763/0.9745/0.8951 and 1.5854/0.9653/0.9633/8432, respectively (Table2). Therefore, for all the sources considered, the width $\Delta\alpha$ value is not near zero at both frequencies, indicating the multifractal and intermittent nature of the sources at both bands, which is in agreement with the conclusion based on the nonlinear scaling exponent function in Figs. 5 and 8. Additionally, the similarity and difference in the degree of the nonlinearity observed at 22 and 37 GHz in the scaling exponent function for each source and between the sources are further confirmed by the calculated width values of the multifractal spectrum function given in the same figures.

We have tested the stability of $\Delta\alpha$ by applying different fluctuations in the range of errors recorded for each observation to the original time series and found that our $\Delta\alpha$’s are stable with uncertainty $<$ 10 %. The observed multifractality could be due to different physical mechanisms. It is the variation in the flux that results in the multifractal structure in a time series, and therefore, the observed multifractal signature could be due to turbulence in the radio radiation region (the blobs that propagate in the jet) since turbulent dynamics create multifractal and intermittent structures (Leonardis, 2013; Yordanova, 2004). Additionally, it has been indicated that a higher magnetic field strength will increase both Compton and synchrotron losses in blazars, which could result in an increase in variability at millimeter and longer wavelengths. Additionally, it has been shown that the change in Doppler factor resulting from the change in shock orientation could result in rapid flux variation (Stevens, 1994), which may culminate in a multifractal structure.

The similarity and difference observed in the degree of nonlinearity at 22 and 37 GHz for each source and between the sources, respectively, provide us with valuable information about the radiation region and mechanisms of the sources considered. The similarity in the slope at those frequencies tells us that the signals of each source at 22 and 37 GHz fluctuate in a similar fashion. The difference in the behavior of the light curves between the sources at both frequencies is clearly visible, and therefore, our finding of different degrees of nonlinearity or multifractality between the sources is somewhat expected.

Indeed, the scale parameter for a time series of given length changes if the length of the time series changes, i.e., the choice of the scale parameter is dependent on the length of the time series. As seen from the calculated width values for each light curve at both bands, we found that the multifractality (nonlinearity) strength at 22 and 37 GHz is strongly similar though not the same for each source and differs from source to source. The similarity in the slope indicates that the sources have the same radiation region and radiation mechanism at the two frequencies as indicated by Wang & Yang (2010) for different light curves of 3C 273 and 3C 345. However, the data streams currently at hand are not of identical length, and the similarity observed at 22 and 37 GHz for each source might be affected by the difference in the data length at these two frequencies. Importantly, when we have full light curves (large data points) at 22 GHz, the scale parameters that we used to obtain better maxima lines for the present data may change based on the change in the singularity behavior of the light curves in time. Consequently, we may have different local maxima lines that affect the current relationship between the scaling exponent function $\tau$($q$) and the moment $q$. This in turn affects the behavior of the multifractality spectrum function $f$($\alpha$). Of course, still, there is difference in nonlinearity at 22 and 37 GHz for all the sources, though not strong, and this could be due to the difference in the number of data points. For the difference in nonlinearity between the sources at both bands, for example, we analyzed the light curves of the sources 3C 273 and 3C 279 at 37 GHz, where there are only 7 data points difference between them using the same scale parameter, and found strongly different results. For the time being, since we are considering scale parameters that give us better maxima lines, informative scale parameters, it is somewhat logical to perform a comparison between the sources in terms of the degree of multifractality at the two frequencies based on what we have on hand. In general, despite the difference in the number of data points (time series length) between the sources, all their light curves at 22 GHz span from 1980 - 2004 and from 1979/1980 -2018 at 37GHz, and therefore, the results obtained using the chosen informative scale parameters could inform us at least about how the multifractal (nonlinear) behavior of the light curves at the two radio bands behave throughout these observation periods.

3.2 Analysis of light curves in the rest frame:

Following the same procedures as applied for light curves in the observation frame, we have repeated the same multifractality analysis for the corresponding light curves in the rest frame, as shown in Fig. 2. We obtained the same results as shown in Figs. 5 and 8 for all the light curves. The similarity in the degree of multifractality (nonlinearity), $\Delta\alpha_{obs}=\Delta\alpha_{rest}$, of the light curves in the observation and rest frames implies that redshift-correction does not affect the multifractal behavior of quasars radio emissions, indicating that multifractality is an intrinsic behavior of quasars radio emissions. The redshift versus multifractality strength is given in Fig. 10.

3.3 Analysis of the ARIMA(p,d,q) models

In this section, the light curves shown in Fig. 1 are analyzed using the ARIMA model. The scenario regarding the ARIMA model analysis of each light curve is discussed hereafter. The first step is verifying whether the time series are stationary or not. A visual inspection of the light curves show that they are not stationary, but one cannot be sure whether a time series is stationary or not only through visual inspection. There are different techniques to test for stationarity such as the augmented Dickey-Fuller unit root test. Though our light curves are shown to be nonstationary in the first subsection of section 3, we have checked the nonstationarity of each light curve based on their ACF plots, which are not included here. The ACF plot of a nonstationary time series decays slowly to zero, whereas stationary time series decay exponentially to zero. As a first step of our analysis, we have plotted ACF for all the time series considered here and found that all of them decay very slowly to zero, implying the presence of nonstationarity. Therefore, we proceed to our analysis directly by differencing, first differencing (d=1), the time series. The differenced time series and corresponding ACF and PACF plots are given in Figs. 11, 12, 13, and 14 for 3C 273, 3C 279, 3C 345, and 3C 454.3, respectively.

As seen from the differenced time series plots, the differencing operator reduces the long-memory autocorrelations, although it left a few significant lags behind in most of the light curves except for 3C 273, where it does not reduce the autocorrelations that much. The next step is model selection and fitting. As a first attempt of best model selection, we estimate the parameters $p$ and $q$ based on the ACF and PACF of the differenced time series. Using estimated values of $p$ and $q$, we fit the time series to the ARIMA(p,d,q) model. The value of d is already known, d=1, since we differenced the time series only once. Additionally, we use the auto.arima() function in our best model selection. We select a model with the lowest log likelihood and the Akaike information criteria (AIC). If the models obtained based on the ACF and PACF plots and using auto.arima() do not fit our data, we try another model by changing the values of the parameters $p$ and $q$ and fit them to the ARIMA(p,d,q) model until we obtain the best one. The last step, in our case, is checking the diagnostic test to determine whether the residuals of the fitted models are white noise or not. We accept or reject the model selected based on two criteria: we accept the model if $p>0.05$ for the Ljung-Box test and the residual is white noise (randomly distributed), indicating no significant lags outside of the 95$\%$ limit, and we reject if these conditions are not satisfied. For the second condition, the residual should be white noise, and if it is not satisfied for all possible combinations of $p$ and $q$, we take the one that nearly provides a white-noise residual and $p>0.05$. Our interest here is not to forecast but only to determine whether ARIMA(p,d,q) models fit our data or not. Except for 3C 273 and 3C 279 at 37 GHz, where the fit to ARIMA models does not reduce the autocorrelations left behind more than 2 significant autocorrelations outside of the 95$\%$ limit, the ARIMA models fit sufficiently to reduce most of the autocorrelations for the rest of the time series. However, the models require a different number of coefficients, and as a result, the residuals of the ARIMA model fits are nearly white noise, leaving behind no more than two significant lags outside of the 95$\%$ limit, as shown from Figs. 11, 12, 13, and 14 for 3C 273, 3C 279, 3C 345, and 3C 454.3, respectively. Except for 3C 345 at 22 GHz, the ARIMA(4,1,4) model fit reduces all significant autocorrelations leaving behind a white noise residual, indicating no single significant autocorrelation outside of the 95$\%$ limit as shown in Fig. 13. Note that the ARIMA(p,d,q) models capture only short-memory components, and the observed failure could be due to the presence of long-memory components and/or nonlinear signatures in the time series that cannot be captured by the ARIMA models. The parametric autoregressive model, called the autoregressive fractionally integrated moving average (ARFIMA) model, where the differencing operator $d$ is a real number, is preferable to model time series with long-memory components. Of course, multifractality analyses as applied in this work are also capable of detecting both short and long memory trends.

4 Summary and Conclusions

In this study, we analyze the multiscaling signatures in the radio emissions of selected radio-loud quasars 3C 273, 3C 279, 3C 345, and 3C 454.3 both in the observed frame (at 22 and 37 GHz passbands) and the rest frame($f_{rest}=f_{obs}*(1+z)$) using a WTMM-based multifractality analysis approach. In addition, we fit the light curves of the sources to the ARIMA models. In our work, we first calculate the wavelet coefficients using a continuous wavelet transform and form a matrix of maxima lines (construct the skeleton function) by aggregating the absolute wavelet coefficients that only hold maxima lines. Second, using the collected local maxima lines and the constructed skeleton function, we determine the thermodynamics partition function. Third, by creating the log-log plots of the thermodynamics function $Zq(s)$ and the scale $s$, we estimate the slope using the least squares fitting method. The behavior of the estimated slopes are presented by the scaling exponent $\tau$($q$) versus $q$ plots in Fig. 3 and Fig. 5 for each source in the observed frame at 22 and 37 GHz, respectively. Finally, we estimate the multifractal spectrum functions at each band for all the light curves and calculate the multifractality strength from the width $\Delta\alpha$ of the spectrum. Additionally, we analyze the corresponding light curves in the rest of the frame following the same procedures. Finally, we fit the time series in the observed frame to the ARIMA models. Our main conclusions are as follows. In this work, we have shown that (i) the scaling nature of quasars’ radio emissions at 22 and 37 GHz is strongly multifractal and intermittent, (ii) the degree of nonlinearity or multifractality is similar for each source and strongly different between sources at both frequencies, (iii) the redshift correction does not affect the nonlinear or multifractal behavior of quasars’ radio emission, and (iv) the ARIMA models fit most of the time series partially, except for 3C 345 at 22 GHz where the ARIMA(4,1,4) model fit obtains a white noise residual and for 3C 273 and 3C 279 at 37 GHz where the models are shown to be inadequate. The strong multifractal signature observed in quasars’ radio emission further supports what has been previously posited: that quasars are intrinsically multifractal and complex systems that have nonlinear time series characterized by fractal behavior (Bewketu Belete et al, 2018; Vio, 1991). It is the physical mechanism that causes variation in the flux of radio emissions that possibly presents the detected multifractal signatures. It has been explained that the low-frequency emissions in general and radio emissions from radio-loud quasars/blazars in particular are due to synchrotron emissions from nonthermal electrons in a relativistic jet (Singal, 2011; Türler, 2000; Robson, 1993; Courvoisier, 1988). This means any physical process that causes a change in the dynamics of relativistic jets and presents variations in the flux of quasars’ radio emissions is likely responsible for the existence of multifractal signatures in the radio observations of quasars. Therefore, our results play a significant role in providing valuable information for those working to develop models to better understand emissions in terms of synchrotron radiation from shocks propagating along relativistic jets and the dynamics of relativistic jets coming out from the center of radio-loud quasars, which in turn helps to constrain some physical properties of quasars in relation to the dynamics of their relativistic jets.

The observed similarity in the slope (degree of nonlinearity) or multifractality strength at 22 and 37 GHz for all the sources considered further supports the claim that the radiations of the sources at 22 and 37 GHz have the same emission region and mechanism, at least for 3C 273 and 3C 345 (Wang & Yang, 2010). Despite that all our sources are flat-spectrum radio-loud quasars, the difference in the degree of multifractality (nonlinearity) between them at those radio bands provides very useful information about the physics of the sources of relativistic jets. This difference could be due to the different nature of turbulence in the accretion rate, internal shocks in the relativistic jets of the sources, fluctuations in the local magnetic fields and particle density, variation in the activity of the central engine, the difference in their black hole mass since it determines the mass accretion rate, which in turn affects the radio emissions, or due to the difference in any other variability mechanism not mentioned here. The other result we have found is that the degree of multifractality (nonlinearity) is the same both in the observed and rest frames of the sources, providing physically important information that multifractality is an intrinsic behavior of quasars’ radio emissions. This finding further supports the conclusion that most of the long-term variations of quasars, in general, are intrinsic to the quasars themselves (de Vries, 2005). In our recent work Bewketu Belete et al (2019), we have shown that extrinsic variations in relation to gravitational lensing, mainly microlensing effects, affect the multifractal behavior of quasars, and the microlensing effect increases the degree of multifractality. Therefore, the results obtained in this work, the sameness in the degree of multiractality in the observation and rest frames, possibly indicate the absence of extrinsic variations, mainly due to microlensing, in the light curves of the sources. A possible reason why the ARIMA(p,d,q) models do not fit most of the time series well could be due to the presence of memories and/or trends that cannot be easily captured by the models. This work provides valuable information mainly to model relativistic jet dynamics in particular and to understand the interior of quasars in general.

Acknowledgments

This publication makes use of data obtained at the Metsähovi Radio Observatory, operated by Aalto University in Finland. Research activities of the Observational Astronomy Board of the Federal University of Rio Grande do Norte (UFRN) are supported by continuous grants from CNPq and FAPERN Brazilian agencies. We also acknowledge financial support from INCT INEspaço/CNPq/MCT. AB acknowledges a CAPES PhD fellowship. ICL acknowledges a CNPq/PDE fellowship. We warmly thank the anonymous referee for the fruitful comments and suggestions that greatly improved this work.