Confidence regions for neutrino oscillation parameters from double-Chooz data

In this work, an independent and detailed statistical analysis of the double-Chooz experiment is performed. To have a thorough understanding of the implications of the double-Chooz data on both oscillation parameters $\sin^{2}(2\theta_{13})$ and $\Delta m^2_{31}$, we decided to analyze the data corresponding to the Far detector, with no additional restriction. By doing this, confidence regions and best fit values are obtained for ($\sin^{2}(2\theta_{13}),\Delta m^2_{31}$). This analysis yields an out-of-order $\Delta m^2_{31}$ minimum, which has already been mentioned in previous works, and it is corrected with the inclusion of additional restrictions. With such restrictions it is obtained that $\sin ^{2}(2 \theta _{13})=0{.}084_{-0{.}028}^{+0{.}030}$ and $\Delta m^2_{31}=2.444^{+0.187}_{-0.215} \times 10^{-3}$ eV$^2$/c$^4$. Our analysis allows us to study the effects of the so called"spectral bump"around 5 MeV, it is observed that a variation of this spectral bump may be able to move the $\Delta m^2_{31}$ best fit value, in such a way that $\Delta m^2_{31}$ takes the order of magnitude of the MINOS value. Finally, and with the intention of understanding the effects of the preliminary Near detector data, we performed two different analyses, aiming to eliminate the effects of the energy bump. As a consequence, it is found that unlike the Far Detector analysis, the Near detector data may be able to fully determine both oscillation parameters by itself, resulting in, $\sin^2(2\theta_{13}) = 0.095 \pm 0.053$ and $\Delta m^{2}_{31} = 2.63^{+0.98}_{-1.15} \times 10^{-3} \text{eV}^2 / \text{c}^4$. The later analyses represent an improvement with respect to previous works, where additional constraints for $\Delta m^2_{31}$ were necessary.

In this work, an independent and detailed statistical analysis of the double-Chooz experiment is performed. In order to have a thorough understanding of the implications of the double-Chooz data on both oscillation parameters sin 2 ð2θ 13 Þ and Δm 2 31 , we decided to analyze the data corresponding to the Far detector, with no additional restriction. This differs from previous analyses, which only aim to estimate the mixing angle θ 13 , without mentioning the effects on Δm 2 31 . By doing this, confidence regions and best fit values are obtained for (sin 2 ð2θ 13 Þ; Δm 2 31 ). This analysis yields an out-of-order Δm 2 31 minimum, which has already been mentioned in previous works, and it is corrected with the inclusion of additional restrictions. With such restrictions it is obtained that sin 2 ð2θ 13 Þ ¼ 0.084 þ0.030 −0.028 and Δm 2 31 ¼ 2.444 þ0.187 −0.215 × 10 −3 eV 2 =c 4 . Our analysis allows us to study the effects of the so-called "spectral bump" around 5 MeV; it is observed that a variation of this spectral bump may be able to move the Δm 2 31 best fit value, in such a way that Δm 2 31 takes the order of magnitude of the MINOS value. In other words, if we allow the variation of the spectral bump, then we may be able to determine both oscillation parameters using Far detector data only, with no further restrictions from other experiments. Finally, and with the intention of understanding the effects of the preliminary Near detector data, we performed two different analyses, aiming to eliminate the effects of the energy bump. As a consequence, it is found that unlike the Far detector analysis, the Near detector data may be able to fully determine both oscillation parameters by itself, resulting in sin 2 ð2θ 13 Þ ¼ 0.095 AE 0.053 and Δm 2 31 ¼ 2.63 þ0.98 −1.15 × 10 −3 eV 2 =c 4 . The later analyses represent an improvement with respect to previous works, where additional constraints for Δm 2 31 were necessary. DOI: 10.1103/PhysRevD.97.093005

I. INTRODUCTION
The double-Chooz experiment estimated the reactor neutrino flux of the Chooz-B nuclear plant by means of its operating parameters. This flux, when interacting with the detector target, induces a number of inverse β decays (IBD). This experiment was designed to run with two detectors located at L F ≈ 1000 m (Far), and L N ≈ 400 m (Near). But the current collaboration results report only far observations. The double-Chooz Near detector was finished in 2016, and only preliminary data have been published until now.
In particular, the double-Chooz Far detector reports fewer IBDs than those expected. If a neutrino oscillations model explains this deficit, then the oscillation parameters can be obtained from double-Chooz data.
In the simplified two-flavor oscillation model, the survival probability of aν e with energy E ν ðMeVÞ after traveling a distance LðmÞ is given as The main objective of the double-Chooz experiment is the precise measurement of the mixing angle θ 13 [1]. Table I shows some results for sin 2 ð2θ 13 Þ. In particular, the double-Chooz collaboration determined sin 2 ð2θ 13 Þ ¼ 0.090 þ0.032 −0.029 , without showing confidence regions and using the value obtained by MINOS of 2.44 þ0.09 −0.10 ×10 −3 eV 2 =c 4 for Δm 2 31 [2].
In this paper we use the double-Chooz data in χ 2 tests to determine the best fit for sin 2 ð2θ 13 Þ and Δm 2 31 without assuming an a priori value for the last one to get both parameters as well as their confidence regions.
Consequently, double-Chooz data may be used for a unified analysis with other experiments where both oscillation parameters are obtained simultaneously by means of their corresponding data.
The organization of this work is as follows. In Sec. I the rate þ shape (R þ S) analysis is presented, aligned with Ref. [2]. The statistics defined there is used for Far detector analysis only. This section shows how the quantities that define the function χ 2 RþS are obtained. A consistent definition of the expected number of IBD is introduced in Sec. II. Both sections are used in Sec. III to estimate the oscillation parameters sin 2 ð2θ 13 Þ and Δm 2 31 by minimizing the function χ 2 RþS . Also the confidence regions for these parameters are obtained from a diagonal covariance matrix (DCM), and a full covariance matrix (FCM). The values obtained for Δm 2 31 are significantly different than those expected. In order to line up our results with the MINOS experiment, a χ 2 RþSþM test is presented. In Sec. IV a discussion related to the spectral bump of the neutrino spectrum at 5 MeV is included to estimate their effect in the χ 2 RþS function. Section V is devoted to the Far þ Near detector data. With the purpose of determining the oscillation parameters Δm 2 31 , sin 2 ðθ 13 Þ and their confidence regions from double-Chooz data without using a priori the value of Δm 2 31 from another experiment, two convenient χ 2 functions are introduced. These statistics also suppress the spectral bump effects mentioned before. To do this, we use the preliminary data from [8] as input to the formalism presented in the previous sections. The results obtained are promising and can be used when the collaboration releases new data.
Finally, our conclusions are given in Sec. VI.
II. RATE + SHAPE ANALYSIS [2] Neutrinos are detected through the positron kinetic energy, E vis , in the energy range of 0.5 and 20 MeV, which is divided into 40 energy bins, accordingly to Table 15.2 in [9]. The rate þ shape analysis is determined by the function In the first term, each energy bin requires N obs i (N exp i ), which is the observed (expected) number of IBD. A covariance matrix M ij is introduced to include the correlation terms among energy bins.
N obs i were directly obtained from Fig. 21 in [2]. In this analysis we considered the 17 351 IBD candidates, N tot , that occurred during 460.67 days, T on . The first thirty-one readings are consistent with previous double-Chooz collaboration data [10].
N exp i is proportional to the expected number of antineutrinos without oscillations, n exp i , and to the average survival probability in the ith energy bin,P i , closely related to the flavor-oscillation model (1). In the following section the explicit form of these quantities is given. For now we can write N exp i as follows: The diagonal matrix elements M ij contain information on statistical and systematic uncertainties in each energy bin. The bin-to-bin correlations correspond to the off diagonal elements. M ij is discussed in more detail in [2] and is written as The matrix elements M ij have been taken from Fig. 15.3 in [9]. In this figure there are five matrices that define the M ij numerical value. Since the information is presented through a color code, a basic software was necessary to decode it. We verified that the diagonal elements M stat ii , M flux ii , and M eff ii previously published in [10] were found in the matrices of [9].
In addition to the uncertainties involved in the covariance matrix, eight systematic uncertainties are considered in the second and third terms of χ 2 RþS , using ϵ x parameters.
0.109 AE 0.030 a AE 0.025 b Abe et al. [5] 0.097 AE 0.034 a AE 0.034 b Abe et al. [6] 0.102 AE 0.028 a AE 0.033 b Novella [7] 0.102 AE 0.043 c Abe et al. [2] 0.090 þ0.032 Three of them are the coefficients of a polynomial for the visible energy variation, These are explicitly introduced in the second term of Eq. (2) by means of a matrix that contains uncertainties σ a , σ b , σ c and correlations ρ ab , ρ bc , and ρ ac ( Table II). The other sources of systematic uncertainties are considered through five parameters ϵ k , whose standard deviations are given in Table III.
The fourth term in (2) or χ 2 off is the contribution of two reactors off (2-off) data, in which N obs off (N exp off ) is the observed (expected) number of IBD candidates.
According to Ref. [2], N obs off ¼ 7 and N exp off is determined by the residualν e 's (ϵ 4 ), and the background (n bg ), as where n bg ¼ B · T off . B ¼ 1.56 events=day is the total background rate provided by T off ¼ 7.24 days of reactor-off data, [9].
For IBD events with neutrons captured on Gd, P off ðθ 13 ; Δm 2 31 Þ denotes the average survival probability on all the spectrum of antineutrinos with the reactor off, and is written as [9,11], as suggested by Eq. (1).L ¼ 1050 m is the average distance from the Far detector to both reactors. The term in angle brackets results from averaging the survival probability (1) over all the energy range, 0.5 MeV ≤ E vis ≤ 20.0 MeV, is the energy of the incomingν e written in terms of the corrected visible energy or positron kinetic energy, E vis . The last term of (9) results directly from the observation of the IBD, and depends on the positron and nucleons masses, III. EXPECTED NUMBER OF IBD, N exp i As introduced in Eq. (3), the expected number of IBD, N exp i , is closely related to the oscillation model and is defined as In this equation the expected neutrino spectrum without oscillations, n exp i , depends on the operating parameters of the nuclear reactor. These quantities are obtained from Fig. 21 at [2]; a discrepancy between the observed and the expected neutrino spectrum without oscillations between 4 and 6 MeV has been detected. This energy bump will impact the determination of the oscillation parameters. It is discussed in the next section.
The residual neutrinos ϵ 4 are produced by the radioactive elements in the core of the nuclear reactors, even when they are turned off. The term has been included to take into account this contribution to the spectrum. Both types of events are influenced by the same average oscillation probability over each energy bin,P i ðθ 13 ; Δm 2 31 Þ. Additionally, the three background sources, ϵ 1 , ϵ 2 , and ϵ 3 , mentioned in the description of Table III, are taken into account in Eq. (10).  [2]. (2). The pull parameters ϵ 1 ; …; ϵ 5 are corrections to the predicted antineutrino spectrum. These are as follows: ϵ 1 , antineutrino spectrum error due to β decays of 8 He and 9 Li; ϵ 2 , error due to n þ μ; ϵ 3 , accidentals; ϵ 4 , residuals; ϵ 5 , uncertainty of the squared mass differences Δm 2 31 . The last one is removed later, [2].
Because we use the double-Chooz data to determine the best fit for sin 2 ð2θ 13 Þ and Δm 2 31 without assuming an a priori value for Δm 2 31 , the pull associated with the correction to this parameter, ϵ 5 , is not required anymore.
Therefore, the χ 2 statistics (2) is a multiparametric function made of two oscillation parameters and seven pulls, The Δχ 2 RþS function is defined as the difference between the χ 2 RþS function and its absolute minimum.
IV. CONFIDENCE REGIONS OF OSCILLATION PARAMETERS θ 13 AND Δm 2

WITH FAR DATA ONLY
We report the minimization of the function χ 2 RþS and its level curves considering the FCM M in Fig. 1. Figure 2 takes into account only the diagonal elements of the covariance matrix M (DCM). In both plots several local minimums are shown.
The absolute minimum (black star) in Fig. 1 has coordinates (0.087, 27.043×10 −3 eV 2 =c 4 ) in the sin 2 ð2θ 13 Þ− Δm 2 31 plane and the χ 2 RþS value of 37.17 for the FCM analysis. Another local minimum (white star) is (0.090, 2.512 × 10 −3 eV 2 =c 4 ) and its χ 2 RþS value is 41.83. These points share close values for sin 2 ð2θ 13 Þ. Nevertheless their corresponding values for Δm 2 31 are significantly different. This effect is closely related to the existence of the energy bump of the neutrino spectrum. These points are listed in Table IV.
The absolute minimum obtained implies Δm 2 31 ¼ 27.043 × 10 −3 eV 2 =c 4 . This value is 1 order of magnitude higher than those reported elsewhere. However, given the quasiperiodic nature of the χ 2 RþS as a function of Δm 2 31 , it can be argued that any one of the local minimums may be the right one, and then additional experimental data and/or improved models would be needed to discriminate between minimums. For this reason, in Sec. V we introduce value. Nevertheless, it is remarkable that a local minimum appears at sin 2 2θ 13 ¼ 0.090, Δm 2 31 ¼ 2.512 × 10 −3 eV 2 =c 4 ; this minimum is included in the confidence region and is consistent with the Δm 2 31 value given by MINOS [12] (see Table IV). The plane sin 2 2θ 13 ¼ 0.087 is indicated with a red vertical line. The intersection of this plane with the χ 2 RþS level curves has been plotted in Fig. 4.

FIG. 2. Behavior of the χ 2
RþS statistics for DCM analysis. The confidence region up to 90% of C.L. for (sin 2 2θ 13 , Δm 2 31 ) shows two disjoint regions. The best fit is found at sin 2 2θ 13 ¼ 0.091, Δm 2 31 ¼ 27.043 × 10 −3 eV 2 =c 4 , again, an inconsistent Δm 2 31 value. This time, the local minimum appears at sin 2 2θ 13 ¼ 0.085, Δm 2 31 ¼ 2.422 × 10 −3 eV 2 =c 4 and still belongs to the confidence region. This local minimum remains consistent with the Δm 2 31 value given by MINOS [12] (see Table IV). A discrimination criterion is needed. The plane sin 2 2θ 13 ¼ 0.091 is indicated with a red vertical line. The intersection of this plane with the χ 2 RþS level curves has been plotted in Fig. 4. data-data analyses that allow us the elimination of the energy bump effect. To do this, Near detector data are required.
The simplest discrimination criterion is that which allows us to establish the value Δm 2 31 , closest to that obtained in another experiment.
Consistently with [4], we introduce to the function (2) the additional term where Δm 2 MINOS ¼ 2.44 × 10 −3 eV 2 =c 4 , and σ MINOS is the average of the Δm 2 31 uncertainties reported by MINOS, Table IV. So, by minimizing  Table V. The confidence regions generated from the χ 2 RþSþM statistics are consistent with those published in [9] and shown in Fig. 3. Figure 4 shows Δχ 2 RþS as a function of Δm 2 13 where sin 2 ð2θ 13 Þ has been fixed at 0.087 and 0.091. These values correspond to the sin 2 ð2θ 13 Þ coordinate of the absolute minimum obtained from the FCM and DCM analyses, respectively (Table IV). A succession of Δχ 2 RþS local minimums appear and are denoted as

V. SPECTRAL BUMP EFFECTS
We can note that the jth minimum has the Δm 2 31 j j coordinate; then In the DCM analysis, the separation between two consecutive minimums of Fig. 4 is given as  (14). As a consequence of the addition of χ 2

MINOS
to the statistics, the absolute minimum is discarded. In this way, the best fit is found at sin 2 2θ 13 ¼ 0.092, Δm 2 31 ¼ 2.444 × 10 −3 eV 2 =c 4 for FCM analysis, and sin 2 2θ 13 ¼ 0.084, Δm 2 31 ¼ 2.444 × 10 −3 eV 2 =c 4 for DCM analysis. The wider region corresponds to the full analysis, and therefore, this one has greater uncertainties (see Table V). For comparison purposes we introduce the Daya Bay data for the parameters sin 2 2θ 13 and Δm 2 31 up to 95.45% of C.L. [13]. All the analyses are consistent to each other.
Besides, we can identify the absolute minimum with χ 2 m 3 ¼ 0 at Δm 2 31 ¼ 27.043 × 10 −3 eV 2 =c 4 and χ 2 m 1 ¼ 3.27, with Δm 2 31 ¼ 2.422 × 10 −3 eV 2 =c 4 , which is closer to the currently accepted Δm 2 31 value. A weighted average of the neutrino energy can be defined asĒ where ω i is the percentage of the observed IBD in each energy bin. Substituting this value into the term sin 2 ð1.27Δm 2 31L =Ē ν Þ, we found that it vanishes when Δm 2 31 ¼ 0.012 eV 2 =c 4 . This value is approximately equal to the average of λ 1 , λ 2 , and λ 3 . The small variation of these values may be attributed to the complex dependence of the χ 2 RþS on the squared sine function and to the average value ofĒ ν used.
The Far detector results alone can be used to discuss how the spectral bump around 5 MeV in the neutrino spectrum affects the Δm 2 31 fit and how the distribution of χ 2 m j might change. RþS profiles when a hypothetical source of rector neutrinos η is added to the prediction in the energy bump, when the unknown source is 5% (black solid line), 10% (red dashed line), and 20% (blue short dashed line) of the total prediction and their respective sin 2 ð2θ 13 Þ best fit. The oscillatory behavior of the χ 2 functions remains, but the Δm 2 31 best fit adopts different values: 2.73 × 10 −2 eV 2 =c 4 at ξ ¼ 5%, 1.79 × 10 −1 eV 2 =c 4 at ξ ¼ 10%, and 8.20 × 10 −3 eV 2 =c 4 at ξ ¼ 20%. The Δm 2 31 best fit is sensitive to the energy bump changes.
The origin of the energy bump is still undetermined, as discussed by [2,14]. The distortion seems to be hardly correlated to the reactor flux. This hypothesis was tested by the double-Chooz collaboration finding that the number of reactors on has influence on the distortion rate.
The pattern of minimums χ 2 m j is sensitive to the energy bump changes. As an example, we introduce a hypothetical source of rector neutrinos, given by η i ≡ ξn exp i added to the prediction in the energy bump. As a consequence, the oscillatory behavior of the χ 2 functions remains, but the Δm 2 31 can diverge 2 order of magnitude higher than expected or even fall into the order set by MINOS Fig. 5. Note that when ξ ¼ 5%, the Δm 2 31 is still 1 order of magnitude higher than expected, and diverges 2 orders of magnitude when ξ ¼ 10%, but when ξ ¼ 20% the difference of squared masses falls into 8.2 × 10 −3 eV 2 =c 4 .
Thus, the effect of the spectrum distortion is relevant but its source is unknown. If we want to obtain both parameters simultaneously, it is necessary to change our point of view, suppressing the energy bump as is discussed in the next section. This is encouraging to perform a unified analysis with other experiments.

VI. CONFIDENCE REGIONS OF OSCILLATION
PARAMETERS θ 13 AND Δm 2 31 WITH FAR + NEAR DATA Although the Near detector was built in 2016, only preliminary results have been published [8]. These preliminary double-Chooz two-detector results can be used as input to the formalism presented in Secs. II and III.
A direct comparison between two sets of data (a data-data analysis) has been considered to cancel the spectrum distortion for the determination of the oscillation parameters.
In order to perform a data-data analysis we are restricted to compare only the Far II to the Near data from [8]. In particular, we propose a χ 2 ð1Þ statistics defined as where N obs Far=Near are the observed number of IBD candidates at the Far/Near detector in the bin with energy E vis i ; Ω i is a weight factor, andP i ðθ 13 ; Δm 2 31 ; L Far=Near Þ is the averaged survival probability over each energy bin at the Far/Near detector.
This statistics suppresses the use of the prediction of the unoscillated reactor neutrino signal spectrum n exp Near=Far at the Near/Far detector.
Another way to define a data-data analysis is  Table VI.
By means of the data-data analyses, the influence of the spectral distortion for the Δm 2 31 determination is highly suppressed. Figure 6 shows the 68.27%, 90%, and 95.45% C.L. regions. Three main points are remarkable.
(i) Data-data analyses no longer show two disjoint regions as the χ 2 RþS in Sec. III, and best fit. Through data-data analyses the spectral bump effect in the Δm 2 31 determination is highly suppressed. Even when the oscillatory behavior is still present, the Δm 2 31 is fully defined now; this is shown in Fig. 7. (ii) Δm 2 31 does not differ by 1 or more orders of magnitude with respect to MINOS, (iii) χ 2 ð1Þ and χ 2 ð2Þ do not depend on external information, as χ 2 RþSþM . Even when both χ 2 ð1Þ and χ 2 ð2Þ statistics still have the oscillatory behavior in Δm 2 31 , there is a well-defined difference between the absolute minimum and the local minimums, as can be seen in Fig. 7.
It is important to recall that in this section, preliminary data from [8] were used. In fact, the formalism described in Secs. II and III, in addition to the data-data statistics of this section, can be used to analyze the double-Chooz twodetector data to determine both sin 2 ð2θ 13 Þ and Δm 2 31 without any restrictions from other experiments.
This work represents a useful tool to build a unified analysis of double-Chooz, Daya Bay, and RENO, as suggested in [15] and [16], even without solving the spectrum bump problem.

VII. CONCLUSIONS
The proposed χ 2 RþS statistical analysis yields consistent results to those published by the double-Chooz collaboration using Far data. The approach followed allows us to generate the confidence regions for the oscillation parameters Δm 2 31 and sin 2 ð2θ 13 Þ shown in Figs. 1 and 2. The effect of the nondiagonal elements of the covariance matrix on the oscillation parameters can be compared with DCM analysis, which is in some sense a zero statistical error case based on the double-Chooz Far detector analysis only.
It is observed that in the FCM analysis the confidence regions are wider in sin 2 ð2θ 31 Þ, and therefore have greater uncertainties.
Each one of the FCM and DCM analyses reports the existence of a χ 2 RþS absolute minimum, corresponding to a Δm 2 31 value, which is inconsistent with the MINOS Δm 2 31 value. In fact, the double-Chooz Far data do not provide by themselves enough evidence to perform a squared mass difference Δm 2 31 estimation with the data available before 2016. However the first local minimum agrees with MINOS Δm 2 31 value, as shown in Fig. 4. In order to force the first local minimum to become the absolute minimum we introduced the additional term (13), in χ 2 RþS . Hence by minimizing the χ 2 RþSþM we got the best fit parameters, sin 2 ð2θ 13 Þ ¼ 0.084 þ0.030 −0.028 and Δm 2 31 ¼ 2.444 þ0.187 −0.215 × 10 −3 eV 2 =c 4 , as can be seen in Table V. In Fig. 3 we have established the confidence regions for neutrino oscillation parameters θ 13 and Δm 2 31 from double-Chooz Far data.
In Sec. IV we have introduced a hypothetical source of rector neutrinos to show how the spectrum distortion affects the oscillatory behavior of the χ 2 functions and the Δm 2 31 value. We found that a correction of 20% in the expected spectrum distortion, independently of its source, corrects the order of magnitude of Δm 2 31 , as indicated in Fig. 5. To cancel the spectrum distortion in the determination of the oscillation parameters, we performed two data-data analyses, using preliminary two-detector data. In both cases, the Δm 2 31 values obtained are not so different than those currently accepted by the community as shown in Figs. 6 and 7 and Table VI. Data-data analyses no longer show two disjoint regions as the χ 2 RþS in Sec. III. Also the value of Δm 2 31 found does not differ by 1 or more orders of magnitude with respect to MINOS and the whole analysis is independent of external information. to their respective sin 2 ð2θ 13 Þ best fit. The spectral bump effect in the Δm 2 31 determination is highly suppressed by means of data-data analyses. The oscillatory behavior is still present, but the Δm 2 31 is fully defined now and agrees with the currently accepted value for this parameter. The horizontal lines at Δχ 2 RþS ¼ 2.3, 4.61, and 6.18 represent the 68.27%, 90.0%, and 95.45% C.L., respectively. Notice that only the absolute minimum falls into these regions and it is near the expected one.
The formalism described in Secs. II and III in addition to the data-data statistics from Sec. V can be used to analyze the two-detector double-Chooz data to determine both sin 2 ð2θ 13 Þ and Δm 2 31 without any restrictions from other experiments. In this way, this work extends the facilities of the double-Chooz experiment by allowing us to measure two oscillation parameters, Δm 2 31 and sin 2 ð2θ 13 Þ.
This work might contain elements of a future unified analysis with other experiments, such as Daya Bay and RENO even with the spectrum bump problem. ACKNOWLEDGMENTS B. V. P. acknowledges the Escuela Superior de Física y Matemáticas, Instituto Politécnico Nacional, for the hospitality during his PhD studies in sciences. We also thank the kind referee for the positive and invaluable suggestions that have improved the manuscript greatly. Special thanks go to Karla Rosita Téllez Girón Flores for her suggestions. This work was partially supported by COFAA-IPN, Grants No. SIP20180062 and No. SIP20170031 IPN and the Consejo Nacional de Ciencia y Tecnilogía through the SNI-México.