Decay of distance autocorrelation and Lyapunov exponents

TL;DR

This paper investigates how the decay of distance autocorrelation in discrete dynamical systems relates to Lyapunov exponents, revealing different decay laws for various systems and linking decay exponents to Lyapunov exponents.

Contribution

It provides numerical evidence connecting the decay of autocorrelation to Lyapunov exponents across different dynamical systems, highlighting system-specific decay laws.

Findings

Exponential decay for quadratic map

Logarithmic decay for Hénon map

Power-law decay for standard map with decay exponent near Lyapunov exponent

Abstract

This work presents numerical evidences that for discrete dynamical systems with one positive Lyapunov exponent the decay of the distance autocorrelation is always related to the Lyapunov exponent. Distinct decay laws for the distance autocorrelation are observed for different systems, namely exponential decays for the quadratic map, logarithmic for the H\'enon map and power-law for the conservative standard map. In all these cases the decay exponent is close to the positive Lyapunov exponent. For hyperbolic conservative systems, the power-law decay of the distance autocorrelation tends to be guided by the smallest Lyapunov exponent.

Tables

Table 1. Table 1: Table presents the models and corresponding positive LEs λ i ​ ( i = 1 , 2 , 3 , 4 ) subscript 𝜆 𝑖 𝑖 1 2 3 4 \lambda_{i}(i=1,2,3,4) with the adjustment curves and the associated decay exponents β 𝛽 \beta .

Model	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\lambda }_4$	$\beta$	$DA\left(\Delta t\right)$
HM	0.418	-	-	-	0.419	$-\beta \text{ln}\left(\Delta t\right)$
SM (mixed)	0.750	-	-	-	0.748	$0.422\Delta t^{-\beta }$
SM (chaotic)	1.505	-	-	-	1.478	$0.407\Delta t^{-\beta }$
2-CSM (mixed)	0.860	0.727	-	-	0.862	$0.430\Delta t^{-\beta }$
2-CSM (chaotic)	1.563	1.469	-	-	1.280	$0.341\Delta t^{-\beta }$
4-CSM (mixed)	0.922	0.853	0.776	0.624	0.741	$0.379\Delta t^{-\beta }$
4-CSM (chaotic)	1.608	1.558	1.500	1.437	1.396	$0.358\Delta t^{-\beta }$

Equations

$$D{C}_N(\mathbf{X},\mathbf{Y})=\frac{{\sigma }_N(\mathbf{X},\mathbf{Y})}{\sqrt{{\sigma }_N(\mathbf{X}){\sigma }_N(\mathbf{Y})}},$$

$${\sigma }_N(\mathbf{X},\mathbf{Y})=\frac{1}{N}{\left(\sum _{i,j=1}^N{A}_{ij}{B}_{ij}\right)}^{1/2},$$

$${\sigma }_N(\mathbf{X})=\frac{1}{N}{\left(\sum _{i,j=1}^N{A}_{ij}^2\right)}^{1/2},$$

$${\sigma }_N(\mathbf{Y})=\frac{1}{N}{\left(\sum _{i,j=1}^N{B}_{ij}^2\right)}^{1/2}.$$

$$A_{ij}=a_{ij}-{\bar{a}}_{i.}-{\bar{a}}_{.j}+{\bar{a}}_{..},$$

$$B_{ij}=b_{ij}-{\bar{b}}_{i.}-{\bar{b}}_{.j}+{\bar{b}}_{..},$$

$$\mathbf{X}\Rightarrow$$

$$\mathbf{Y}={\mathbf{X}}^{\prime }\Rightarrow$$

$$DA=D{C}_N(\mathbf{X},{\mathbf{X}}^{\prime })=\frac{{\sigma }_N(\mathbf{X},{\mathbf{X}}^{\prime })}{\sqrt{{\sigma }_N(\mathbf{X}){\sigma }_N({\mathbf{X}}^{\prime })}}.$$

$$x_{n+1}=r-x_n^2,$$

$$\begin{array}[]{ll}x_{n+1}=r-x_{n}^{2}+by_{n},\\ y_{n+1}=x_{n}.\end{array}$$

$$\begin{array}[]{ll}p_{n+1}=p_{n}+K\sin(2\pi\hskip 0.56917ptx_{n})&\hskip 11.38092pt[\mathrm{mod}\hskip 8.5359pt1],\\ x_{n+1}=x_{n}+p_{n+1}&\hskip 11.38092pt[\mathrm{mod}\hskip 8.5359pt1],\\ \end{array}$$

$$\mathbf{M_{i}}\left(\begin{array}[]{c}p_{i}\\ x_{i}\\ \end{array}\right)=\left(\begin{array}[]{llll}p_{i}+K_{i}\sin(2\pi x_{i})&\hskip 2.84544pt[\mathrm{mod}\hskip 5.69046pt1]\\ x_{i}+p_{i}+K_{i}\sin(2\pi x_{i})&\hskip 2.84544pt[\mathrm{mod}\hskip 5.69046pt1]\\ \end{array}\right),$$

$$\mathbf{T_{i}}\left(\begin{array}[]{c}p_{i}\\ x_{i}\\ \end{array}\right)=\left(\begin{array}[]{llll}p_{i}+\sum_{j=1}^{N}\xi_{i,j}\hskip 1.42271pt\sin[2\pi(x_{i}-x_{j})]\\ x_{i}\\ \end{array}\right),$$

$$|X_i-X_j|=\sqrt{\sum _k^N\left[{\left(x_i^{(k)}-x_j^{(k)}\right)}^2+{\left(p_i^{(k)}-p_j^{(k)}\right)}^2\right]}.$$

Full text

I Introduction

The understanding of decay of correlations in dynamical systems is of crucial relevance in nonequilibrium statistical physics. It may furnish the qualitative and quantitative knowledge about, for example, relaxation processes in complex systems Falcioni et al. (1990). In such complex systems a chaotic behaviour is common and a relation between the decay of correlations with positive Lyapunov exponents (LEs) is intuitively expected. From the theoretical point of view such relation was studied before using the standard autocorrelation function between two observables, as for example, in 1D dissipative systems Badii et al. (1988); Slipantschuk et al. (2013a, b), in area preserving maps Crawford and Cary (1983), in 2D maps Badii et al. (1988) with random perturbations Blumenthal et al. (2018), and conjectured that for higher-dimensional systems the smallest positive LE is an upper bound for the decay of correlation Collet and Eckmann (2004).

In this work we study the decay of a more recently Székely et al. (2007) proposed correlation, named the distance correlation or covariance correlation. It was proposed for testing joint independence of random vectors in arbitrary dimensions and has received attention in applied statistics Székely and Rizzo (2009, 2012, 2013, 2014, 2017), in quantum-classical transitions analysis for ratchet systems Beims et al. (2015), in genetic risk problems Kong et al. (2015), in functional brain connectivity Geerligs et al. (2016) and for the parameter identification of nonlinear systems Ayala Solares and Wei (2015). It was shown recently Mendes and Beims (2018) that the distance correlation is able to detect noise induced escape times from regular and chaotic attractors, to describe mixing properties between chaotic trajectories and to detect properties related to the linear stability of orbital points. However, no direct relation to the LE was achieved.

The crucial issue is to extract from the distance correlation relevant properties related to nonlinear dynamics hidden in time series. Inspired by this, the present work aims to give a qualitative and quantitative connection between the distance correlation and LEs. Being more specific, we show that the positive LE is closely related to the decay of the distance correlation calculated between the time series and itself delayed in time, named here as distance autocorrelation $(DA)$. Results are presented for two dissipative maps, the quadratic and the Hénon maps and for the conservative standard map and coupled standard maps. Distinct qualitative decays of $DA$ as a function of the delay time $\Delta t$ are observed for these maps, ranging from exponential, logarithm to power-law. In all cases the decay exponents are close to the LEs. Exceptions only occur for hyperchaotic systems, for which the decay exponents tend to follow the smallest LEs.

The paper is organized as follows. While in Sec. II we review the general definition of the distance correlation, Sec. III describes the $DA$ used in this work. Section IV shows how the $DA$ can be used to determine the periodicity of a time series and Sec. V shows its relation to the LEs from the quadratic map. In Sec. VI we test our findings for higher-dimensional systems, namely the Hénon map and (coupled) standard maps. Section VII summarizes our main results.

II Distance correlation

In this Section we summary the main properties of the distance correlation Székely et al. (2007) as a statistical measure of dependence between random vectors which is based on Euclidean distances. We present here only the computational procedure to determine the distance correlation. For more details about this procedure and original definitions we refer the readers to Székely et al. (2007).

Assume two random samples $(\boldsymbol{X},\boldsymbol{Y})=\{(X_{k},Y_{k}):k=1,\dots,N\}$ for $N\geq 2$ and $X\in\mathbb{R}^{s},Y\in\mathbb{R}^{t}$ with $s$ and $t$ integers. The empirical distance correlation is given by the expression

[TABLE]

where the empirical distance covariance $\sigma_{N}(\boldsymbol{X},\boldsymbol{Y})$ for a joint random sample $(\boldsymbol{X},\boldsymbol{Y})$ is defined by the expression

[TABLE]

where $A$ and $B$ are matrices. The empirical distance variance for a random sample $\boldsymbol{X}$ is given by

[TABLE]

and for a random sample $\boldsymbol{Y}$

[TABLE]

For $X\in\mathbb{R}^{s}$ with $i=1,\ldots,N$ and $j=1,\ldots,N$, the matrix $A$ is obtained from

[TABLE]

where $a_{ij}=|X_{i}-X_{j}|$ is the Euclidean norm of the distance between the elements of the sample, $\bar{a}_{i.}=\frac{1}{N}\sum_{j=1}^{N}a_{ij}$ and $\bar{a}_{.j}=\frac{1}{N}\sum_{i=1}^{N}a_{ij}$ are the arithmetic mean of the rows and columns, respectively, and $\bar{a}_{..}=\frac{1}{N^{2}}\sum_{i,j=1}^{N}a_{ij}$ is the general mean. Similarly to $Y\in\mathbb{R}^{t}$ defined by $i=1,\ldots,N$ and $j=1,\ldots,N$, we also define the matrix

[TABLE]

where the terms $b_{ij}$, $\bar{b}_{i.}$, $\bar{b}_{.j}$ and $\bar{b}_{..}$ are similar to those for the matrix $A_{ij}$.

The $DC_{N}(\boldsymbol{X},\boldsymbol{Y})$ is defined inside the interval $[0,1]$ and its main characteristic is that it will be zero if and only if the random vectors are independent Székely et al. (2007); Székely and Rizzo (2009, 2012, 2013, 2014). In addition, it is easy to check that $DC_{N}(\boldsymbol{X},\boldsymbol{Y})$ is scale independent. This means that samples $\boldsymbol{X}$ and $\boldsymbol{Y}$ can be multiplied by $\alpha$ and $\beta$, with $\alpha,\beta\in\mathbb{R}$, and $DC_{N}(\boldsymbol{X},\boldsymbol{Y})$ remains unaltered. The $DC_{N}(\boldsymbol{X},\boldsymbol{Y})$ only requires a series of data in order to be analyzed. This can be crucial when analyzing experimental data.

III Distance autocorrelation

Defining matrix $B$ to be equal to matrix A with a delay $\Delta t$, the matrix elements are related by $B_{i,j}=A_{i+\Delta t,j+\Delta t}$. In the present work we always consider time series, so that $\Delta t$ is a time delay. Using increasing values of $\Delta t$ this is a discrete convolution between $A$ and $A_{i+\Delta t,j+\Delta t}$ and we expect to obtain some relevant information about the dynamics of the original time series as a function of $\Delta t$. In other words, for the calculation of $DA$, we use here only one time series of states of a particular initial condition (IC) to compose the sample $\{\boldsymbol{X}\}$. From this time series, we obtain the data set $\{\boldsymbol{Y}\}$ by only shifting the states of $\{\boldsymbol{X}\}$ on time. Thus, we must consider the following representation for the joint sequence $(\boldsymbol{X},\boldsymbol{Y})$

[TABLE]

Here $n$ is an integer referring the times for which the variables $(\boldsymbol{X,Y})$ were determined. For maps, $1\leq\Delta t\leq N$ is a positive integer. Thus, distance autocorrelation function is then defined by

[TABLE]

IV Detecting periodicity

It is easy to realize that the $DA$ calculated between samples $\boldsymbol{X}$ and $\boldsymbol{X^{\prime}}$ will be $DA=1$ when $\Delta t=m$, with $m$ being the periodicity of the time series. To check this obvious property we present numerically results only for one time discrete dynamical system, the quadratic map.

Quadratic map (QM). The quadratic map is a one-dimensional dynamical system given by the expression Alligood et al. (1996); Gilmore and Lefranc (2011)

[TABLE]

where $r$ is a control parameter that belongs to the range $[-0.25,2]$, $n=0,1,2,3,\ldots$ the integer which represents the number of iterations and $x_{n}\in[-2,2]$ is the state of the system at time $n$. Results for $DA$ as a function of $\Delta t$ are presented in Fig. 1. $DA$ is shown for two periodic windows located at the parameters: $r=1.47786$ (period-6) and $r=1.77043$ (period-3). We notice that in both cases $DA=1$ when $\Delta t=6$ and $\Delta t=3$, respectively. This means that $DA$ furnishes a measure of how much the time series repeats itself. Since this is an evident property we will not show it in other dynamical systems.

V Relation to Lyapunov exponents

In this Section we show that the $DA$ calculated in Sec. III is related to the maximal LE. We mention that while the LE gives the exponential divergence of nearby trajectories in the linear approximation, the $DA$ may contain in itself nonlinear extensions of this property and is not restricted to furnishes solely the usual LE.

To start the discussion we have firstly to guarantee that the $DA$ converges to a reasonable value for a given $\Delta t$. This is presented in Fig.2(a) for the QM with $r=1.6$, $1.7$, $1.9$ and $2.0$. In this case, $\Delta t=1$. This shows that for each value of $r$ the $DA$ converges asymptotically, as a function of $n$, to distinct finite values. In other words, even for a chaotic trajectory, $DA$ does not converges to zero. Such situation changes drastically when $DC_{N}(\boldsymbol{X},\boldsymbol{Y})$ is calculated between two distinct ICs inside the chaotic regime. This leads to mixing properties for which $DC_{N}(\boldsymbol{X},\boldsymbol{Y})$ tends to zero with $n^{-1/2}$, as shown in Mendes and Beims (2018).

The next step is to show how such asymptotic values of $DA$ change with $\Delta t$. This is presented in Fig. 2(b), which displays $DA$ as a function of $\Delta t$ for parameters inside the chaotic regimes of the QM, namely $r=1.6$, $1.7$, $1.9$ and $2.0$. No values $DA=1$ are observed since the chaotic trajectory is not periodic anymore. We note that all $DA$ curves converge to small, but finite values when $\Delta t=20$.

From Fig. 2(a) we observe that when increasing $r$ the corresponding $DA$ decrease. Even though we expect that for time series with larger LEs the $DA$ is smaller, we wonder if there is an analytical relation between both. To find such relation we firstly determine $DA$ for many values of parameters inside the interval $1.25\leq r\leq 2.0$ of the QM. This is shown in Fig. 3(a). In this case we considered just two time delays, $\Delta t=1$ (blue curve) and $\Delta t=2$ (orange curve) from a total iteration time of $n=10^{4}$. A direct comparison between $DA$ and the LE ($\lambda$) from the QM, shown in Fig.3(b), allows us to realize that, besides the regions with periodic motion, $DA$ decreases while $\lambda$ increases with the increment of $r$. It becomes clear that $DA$ is able to identify the periodic windows that occur at specific values of $r$.

Now we combine data for $DA$ from Fig. 3(a) with those for $\lambda$ from Fig. 3(b). This is plotted in Fig. 4, and shows the relation between the LE and $DA$ for $\Delta t=1,2$. Clearly almost all positive LEs belong to decreasing lines when increasing $DA$. Linear adjusted lines are given in the inset of this figure. Since the above analysis depends on the specific $\Delta t$ used, it is not appropriate to extract general nonlinear properties from $DA$. Thus, we move forward and argue that there should exist some relation between the LEs and the time decay of $DA$ as a function of $\Delta t$. This argument is reasonable due to the fact that the $DA$ calculated for the time series with itself for distinct values of $\Delta t$ is related to the convolution between both series.

To check this, we obtained the samples $\{\boldsymbol{X}=x_{n+1}:x_{1},x_{2},\dots,x_{N}\}$ and $\{\boldsymbol{X^{\prime}}=x_{n+1}:x_{1+\Delta t},x_{2+\Delta t},\dots,x_{N+\Delta t}\}$ by iterating the QM (8) $n=10^{4}$ times after to discard the first $5\times 10^{3}$ transitory iterations. Thus, to determine $DA$ we calculate an average over $5\times 10^{3}$ ICs randomly chosen with equal probability inside the interval $[-2,2]$ to reach a satisfactory level of accuracy in the asymptotic $DA$ values. We write the mean distance autocorrelation as $\langle DA\rangle$. Figure 5 shows $\langle DA\rangle$ (red circles) as function of $\Delta t$ in semi-log scale for three different values of the parameter $r$ in the chaotic regime. We consider the displacement only in the interval $1\leq\Delta t\leq 10$. In Fig. 5(a) we display the case $r=1.7$, which has a LE $\lambda=0.438$. In Figs. 5(b) and 5(c) the cases $r=1.9$ and $r=2.0$ are shown, for which $\lambda=0.548$ and $\lambda=0.693$ are found, respectively. Straight lines in these curves can be fitted by exponential decays of the type $\langle DA\rangle\propto\alpha\exp{(-\beta\Delta t)}$, where $\alpha$ and $\beta$ are adjustment parameters. Remarkably, the parameters $\beta=0.439$, $\beta=0.541$ and $\beta=0.692$ provides similar values as the corresponding LEs $\lambda=0.438$, $\lambda=0.548$ and $\lambda=0.693$, as can be observed from the adjusted curves (dashed black line) in Fig. 5. We have checked the relation $\beta\approx\lambda$ for other parameters of the QM. This is displayed in Fig. 6(a), which shows the LE $\lambda$ (continuous black line) and $\beta$ (circles) as a function of $r$. From this figure, we can conclude that the relation $\beta\approx\lambda$ is valid for the values of $r$ that lead to chaotic dynamics. Figure 6(b) combines data for $\lambda$ and $\beta$ in the plane $\lambda\times\beta$. These results show that the decay of $DA$ as a function of $\Delta t$ is closely related to the LE of the dynamical system. We have also tried to increase the number of samples, *i. e. * iteration $n$ and ICs, but no relevant changes were observed.

VI Extension to higher-dimensional systems

In the following we analyze the relation of the decay of $DA$ with the LE in nonlinear systems with more dimensions. We start with the dissipative Hénon map and then consider the conservative (coupled) standard map.

Hénon map (HM). The Hénon map is a generalization of the quadratic map for two dimensional systems, and it is given by Hénon (1976)

[TABLE]

The additional parameter $b$ determines the dissipation of the systems. The HM is dissipative for $|b|<1$. This two-dimensional map has two LEs. It is well known that the HM has a chaotic attractor for $r=1.4$ and $b=0.3$ Hénon (1976) with positive LE $\lambda=0.418$, as shown in Fig. 7(a). Since we are considering two dimensions, the quantities for $DA$ are obtained from $|X_{i}-X_{j}|=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}$. Results are displayed in Fig. 7(b) and show the corresponding behavior of $DA$. When we refer to the adjusted curve we use the notation $DA(\Delta t)$. As seen in this figure, the decay of $DA$ is well adjusted by the function $DA(\Delta t)\sim 0.956-0.419\hskip 2.84544pt\mbox{ln}(\Delta t)$. Thus the time decay coefficient $0.419$ is close to the LE $\lambda=0.418$.

Standard map (SM). Here we consider the two-dimensional conservative SM given by Chirikov (1979)

[TABLE]

where $x_{n}$ is the position at the iteration $n=0,1,2,\ldots$, and $p_{n}$ its conjugated momentum. $K$ is the nonlinear positive parameter, and its value determines the topology of the phase space. Some values of $K$ lead to a mixed phase space compound by Kolmogorov-Arnold-Moser (KAM) tori and a stochastic region Lichtenberg and Lieberman (1992). The chaotic trajectory can be trapped for long, but finite times close to the tori, given origin to the stickiness effect Meiss and Ott (1985). This scenario can be obtained using $K=0.57$, for which $\lambda=0.750$, as shown in Fig. 7(c). A totally chaotic phase space takes place if $K=1.43$, with $\lambda=1.505$ as displayed in Fig. 7(e). In both cases the appropriate adjustment was $DA(\Delta t)\propto\alpha\Delta t^{-\beta}$. Also here the quantities for $DA$ are obtained from $|X_{i}-X_{j}|=\sqrt{(x_{i}-x_{j})^{2}+(p_{i}-p_{j})^{2}}$. For the mixed case we have $\alpha=0.422$ and $\beta=0.748$, which is very close to the associated LE [see Fig. 7(d)], and for the chaotic case we found $\alpha=0.407$ and $\beta=1.478$ [Fig. 7(f)], which is also close to the respective LE $\lambda=1.505$.

Coupled Standard Maps (CSMs). In this model we consider a coupling between identical conservative SMs. For our numerical investigation we used the $2$-dimensional SM Altmann and Kantz (2007); da Silva et al. (2015):

[TABLE]

and for the coupling

[TABLE]

with $\xi_{i,j}=\xi_{j,i}=\frac{\xi}{\sqrt{N-1}}$ (all-to-all coupling). This constitutes a $2N$-dimensional Hamiltonian system and the total map is a composition of $T$ and $M$. We considered two cases: (i) $N=2$ and (ii) $N=4$, for which $DA$ is obtained from

[TABLE]

Case (i). Here we have two positive LEs ($\lambda_{1}=0.860$ and $\lambda_{2}=0.727$) for the mixed dynamics obtained using $K_{1}=0.57$ and $K_{2}=0.59$. Calculating $DA(\Delta t)$ for this case, we obtain $\beta=0.862$ [see Fig. 8(a)]. Considering the chaotic case with $K_{1}=1.43$ and $K_{2}=1.45$, we have LEs $\lambda_{1}=1.563$ and $\lambda_{2}=1.469$ and the parameter that fits the decay of $DA(\Delta t)$ is $\beta=1.280$, as seen in Fig. 8(b). In either case the function used to the adjustment was $DA(\Delta t)\sim\alpha\Delta t^{-\beta}$.

Case (ii). Here we have four positive LEs ($\lambda_{1}=0.922,\lambda_{2}=0.853,\lambda_{3}=0.776$ and $\lambda_{4}=0.624$) for the mixed dynamics obtained with $K_{1}=0.57$, $K_{2}=0.58$, $K_{3}=0.59$ and $K_{4}=0.60$. In this case, again we found for the adjustment the function $DA(\Delta t)=\alpha\,\Delta t^{-\beta}$, with $\beta=0.741$, as can be seen in the Fig. 8(c). For the chaotic case, obtained with $K_{1}=1.42$, $K_{2}=1.43$, $K_{3}=1.44$ and $K_{4}=1.45$, we have LEs $\lambda_{1}=1.608,\lambda_{2}=1.558,\lambda_{3}=1.500$ and $\lambda_{4}=1.437$. Results are shown in Fig. 8(d) and we obtain $\beta=1.396$.

Results for higher-dimensional systems are summarized in Table I. The decay of $DA$ in higher-dimensional systems with more than one positive LE does not give precise information about the LEs. In such cases $\beta$ is close to the smallest LE, a feature also observed for the autocorrelation function Collet and Eckmann (2004).

VII Conclusions

The present work analyzes numerically the relation between the decay of the distance autocorrelation and LEs. Results are shown for the dissipative quadratic and Hénon map and for the conservative standard maps and coupled standard maps. For all conservative cases we found a decay for $DA$ proportional to $\Delta t^{-\beta}$, where $\beta$ is very close to the LE $\lambda$ in case the system has just one positive LE. For hyperchaotic systems, the parameters $\beta$ tends to follow the smallest LE (with one exception). For the 1D dissipative system (quadratic map) the observed decay obeys $e^{-\beta\Delta t}$ with $\beta\sim\lambda$. For the 2D dissipative system (Hénon map), the decay follows $-\beta\,\mbox{ln}(\Delta t)$, and $\beta\sim\lambda$. Thus, for systems with one positive LE the decay of the distance autocorrelation is nicely related to the Lyapunov exponent.

Further investigations will study in more details the decay of $DA$ in hyperchaotic and in time continuous systems. Moreover, it should be possible to extract from $DA$ more general aspects of the dynamics, since $DA$ certainly contains information which go beyond the linear stability analysis of the exponential divergence of trajectories, quantified by the LE.

acknowledgments

CFOM thanks FAPEAM, RMS thanks CAPES and MWB thanks CNPq for financial support. The authors also acknowledge computational support from Professor Carlos M. de Carvalho at LFTC-DFis-UFPR.