Late Time Attractors of Some Varying Chaplygin Gas Cosmological Models

: The goal of this paper is to study new cosmological models where the dark energy is a varying Chaplygin gas. This speciﬁc dark energy model with non-linear EoS had been often discussed in modern cosmology. Contrary to previous studies, we consider new forms of nonlinear non-gravitational interaction between dark matter and assumed dark energy models. We applied the phase space analysis allowing understanding the late time behavior of the models. It allows demonstrating that considered non-gravitational interactions can solve the cosmological coincidence problem. On the other hand, we applied Bayesian Machine Learning technique to learn the constraints on the free parameters. In this way, we gained a better understanding of the models providing a hint which of them can be ruled out. Moreover, the learning based on the simulated expansion rate data shows that the models cannot solve the H 0 tension problem.


Introduction
In modern cosmology, there are several key open problems and various approaches to solve them including gravitational particle creation and modification of general relativity [1][2][3][4][5][6][7][8][9][10][11][12][13][14] (and references therein). Recently another one known as the H 0 tension problem has been added to this list [15][16][17][18]. The goal of this paper is (1) to consider various new cosmological models explaining the late time accelerated expansion of the universe [19][20][21][22][23][24][25][26], and (2) to see whether or not the models solve the H 0 tension problem [27][28][29][30][31][32][33][34][35][36][37]. The analysis of the models is based on two approaches. In particular, we use phase space analysis in the first part of the paper, while in the second part of the paper we use Bayesian Machine Learning to learn the constraints on the model parameters (see for instance [32,35,37,38]). The phase space portrait of a model is one of the mechanisms providing a qualitative understanding of the model. It allows to observe all states of the model without solving a system of differential equations, but solving algebraic equations [39][40][41][42][43][44][45][46][47] to mention a few. On the other hand, Bayesian Machine Learning based on the generative process allows one to infer crucial properties of the model directly from the model used in the generative process. This method was recently applied in several studies indicating very interesting departures from traditional approaches used in cosmology. We will come back to this in the second part of the paper when we discuss learned constraints on the parameters.
In general, dark fluids are actively used to explain the accelerated expansion of the universe. Chaplygin gas is one of such dark fluids [48][49][50][51] where A, B and α are positive constants to be constrained from observational data and ρ de it is the energy density of the gas (or dark energy fluid). In literature there are various modifications of this fluid too [52][53][54][55][56][57][58] (to mention a few).
In this paper, we will consider one of them, assuming that the parameter B in Equation (1) is not constant. In general, a varying Chaplygin gas can be constructed assuming both A and B parameters in Equation (1) are not constant. However, in each case, it is important to have a modification either based on some very well motivated physics or obtain the constraints in order to justify the crafted phenomenological modification. In our work, we follow to the second approach and having phenomenological modifications we will learn which one can survive. It should be mentioned that the interest towards Chaplygin gas is its dark energy and dark matter unifying ability. However, some critics on this issue also exist (see for instance [51]).
In our analysis we consider FRW universe with (c = 8πG = 1) where ρ it is the energy density of the effective fluid, while P it is the pressure. It is easy to see from the structure of Equations (2) and (3), that, for example, an assumption about the effective fluid will allow closing the system of differential equations. In general, we can assume that the effective fluid can be multicomponent one that needs to be confirmed the observational data. However, here we follow a simplified scenario assuming that the energy source entering in Equations (2) and (3) is two-component one and that and P = P dm + P de , where ρ de and P de are the energy density and pressure of dark energy, while ρ dm and P dm are the energy density and pressure of dark matter. Now the question is how to be with the dark energy and dark matter. A phenomenological assumption about the content of the universe is a commonly accepted approach in modern cosmology, which will be used in this paper as well. There are various assumptions about the nature and origin of dark energy and each of them had been investigated from various perspectives. On going research towards dark fluid representation of dark energy, for instance, allows to develop class of viscous dark fluids equally applicable to cosmic inflation and the accelerated expansion of the late time universe [59] (and references therein). Moreover, it has been shown that dark fluids can be described by linear and non-linear EoS. On the other hand, it appears that dark fluids can be solutions of algebraic and differential equations [60].
Besides dark fluids we can use, for instance, scalar field to construct viable dark energy models [61] (and references therein). Quintessence and phantom dark energy models are among them. However, our not complete understanding of the universe shows that dark energy independent of its way of interpretation can be not enough. In this regard, various studies allowed to introduce a phenomenological idea of interacting dark energy models demonstrating that it can be equivalently good as the other ideas. Practically there is nothing against to this idea. An interaction between dark components of the universe can smooth some unpleasant aspects in the dynamics of the models making them more attractive. The main mechanism to implement a non-gravitational interaction is based on the energy transfer between dark components in the following waẏ where Q is the notation of non-gravitational interaction [14,39, (Non-gravitational interaction introduced in Equations (6) and (7) indicates energy transfer between dark matter and dark energy). There are different phenomenological parameterizations of the interaction term and some of them can help to solve/alleviate the H 0 tension problem too (see references of this paper for more details about interacting dark energy models). Up to this point we discussed our goal and shortly mentioned the tools we will use. Then, we presented how to describe interacting dark energy models and what is the background dynamics. But yet we did not discuss the varying Chaplygin models we will consider. We will consider two particular models of varying Chaplygin gas (see for instance [52,58]) and where H is the Hubble parameter, a is the scale factor, while n is a constant, to be the dark energy. Eventually, we see from Equations (6) and (7) that another crucial aspect allowing us to study the cosmological models is the interaction term Q. In other words, it should be given in order to close the system of differential equations describing the background dynamics of the universe with interacting dark energy models. In this paper, we consider non-linear interactions Q to be presented in Section 2. It should be mentioned, that we assume dark matter is a pressureless fluid.
To end this section, let us mention that the phase space analysis allows immediately demonstrate, that considered forms of non-gravitational interactions allow solving the cosmological coincidence problem. Besides the phase space analysis we applied Bayesian Machine Learning approach and learned the constraints on the parameters of each model. It gives a hint that within considered models the H 0 tension cannot be solved. Moreover, we found which of the models eventually needs to be ruled out.
The paper is organised as follows: In Section 2, we will discuss how the phase space analysis can be implemented to find late time attractor solutions for suggested new cosmological models. In Section 3, the phase space analysis is performed, late time attractors are found and classified according to their cosmological applicability. Moreover, in Section 4, for each model the constraints on the parameters using Bayesian Machine Learning have been learned. In the same section very briefly the crucial aspects of the approach have been discussed too. Finally, discussion on obtained results are summarised in Section 5.

Interacting Models and Autonomous System
In the literature there is a huge number of works devoted to the phase space analysis of various cosmological models. In order to start the phase space analysis of a model, an appropriate autonomous system should be found first which is a system of algebraic equations to be solved. The critical points are solutions of the autonomous system. They are stable if appropriate Jacobian matrix has a negative trace and positive determinant. This is in case of linear stability. Following [83] for our models we set and where a it is the scale factor. It is not hard to see, that for interacting models the autonomous system reads as and where is the derivative with respect to N, dot is the derivative with respect to cosmic time, Q is interaction term. Explicit forms of Equations (14) and (15) are obtained when the form of Q is given. It is easy to see, that physically reasonable solutions should satisfy 0 ≤ x ≤ 1 (0 ≤ x de ≤ 1) and 0 ≤ z ≤ 1 (0 ≤ z de ≤ 1) constraints. On the other hand, a stable critical point will be an attractor, which we are looking for. At the same time we should remember that x and z according to Equations (2), (10) and (12) should satisfy to the following constraint Moreover, after some algebra we can see that the EoS parameter of the varying Chaplygin gas in terms of x and y reads as while the EoS parameter of the effective fluid reads as In the above equation, we used the fact that dark matter is cold and pressureless. On the other hand, it is not hard to show that the deceleration parameter q reads as Now, in order to obtain the explicit forms of Equations (14) and (15) we need to define the form for Q. As we mentioned earlier we will consider new non-linear non-gravitational interactions. What we consider here can be obtained from the following more general form of interaction Q where ρ could be either the energy density of the effective fluid, or the energy density one of the components of the effective fluid.ρ it is product of the energy densities of the components of the effective fluid i.e., three possibilities could be ρ 2 de , ρ 2 dm and ρ de ρ dm . In this way we can note that the equations discussed in this section are self consistent and allow performing the phase space analysis. In the next section we will discuss our results.

Phase Space Analysis
As we mentioned earlier, in order to have late time attractors for the cosmological models, we need the explicit form of non-gravitational interaction Q to have Equations (14) and (15) determined. In this work we will pay attention to fixed sign non-linear interactions following from Equation (20). Non-linear interactions having similar structure as here, including also models with non-linear sign-changeable interactions already have been considered in [86]. At this stage of the analysis to simplify the discussion, we impose some constraints on the model parameter taken from the literature. In particular, we will impose constraints. On the other hand, the constraints on the deceleration parameter q and on the EoS parameter of the varying Chaplygin gas gave us an option to reduce the phase space size allowing to find late time scaling attractors more efficiently. On the other hand, it is not excluded that the phantom divide about z ≈ 0.2 redshift took place in our universe, therefore constraint Equation (25) is considered [91][92][93]. Having imposed discussed constraints we deal with conditional attractors. They allow obtaining constraints depended late time states of the universe. The goal of this section is to find and discuss stable critical points preparing appropriate initial seeds allowing to initialize Bayesian Machine Learning-based analysis of the models.

Varying Chaplygin Gas
We start the study from the models where the varying Chaplygin gas, Equation (8), interacts with cold dark matter and the the forms of non-linear interactions are defined from Equation (20).
The first cosmological model we study is a model where the interaction between the varying Chaplygin gas, Equation (8), and cold dark matter is given as and real critical points are presented in Table 1. Table 1. Critical points corresponding to interacting varying Chaplygin gas, Equation (8), for the non-gravitational interaction Q given by Equation (26).

S.P. x y Type of Stability
tending to a constant. For this model the decelerated parameter is q = −1 and ω e f f = −1, while the EoS parameter of the varying Chaplygin gas Equation (8) reads as It is not hard to see that considering varying Chaplygin gas, Equation (8), is a phantom dark energy. Discussed behavior is obtained when the parameters satisfy to Equations (21) and (22) and 0 < b < 2/3. The critical point E.1.3 will be stable when 0 < α ≤ 1, 0 ≤ n ≤ 5, b = 0 and A > 1 + α + n/2. This solution describes a matter dominate state of the universe, since x = 0, and in such universe the accelerated expansion is not possible. Therefore, only E.1.1 will be the physically reasonable solution where the universe will reach starting from the state described by E.1.2.

Interaction
In the second model the interaction between varying Chaplygin, gas Equation (8), and cold dark matter is taken to be In this case, physically reasonable two critical points exist and they are presented in Table 2. The study shows that E.2.1 is a physically reasonable solution, moreover, it is a late time scaling attractor when 0 ≤ n ≤ 5, 0 < α ≤ 1, A ≥ 0 and 0 < b < 2/3. This solution represents a state of the universe where varying Chaplygin gas, Equation (8), is a phantom dark energy with This model is free from the cosmological coincidence problem due to Critical point E.2.2 is physically reasonable solution when b = 0 i.e., x = 1 and y = A (when A > 0), and this is a state of the universe where the accelerated expansion is not possible. Moreover, the varying Chaplygin gas, Equation (8), is a usual fluid and completely dominates the dynamics of the universe. As was expected E.2.2 is not a stable critical point. Figure 1 represents the phase space portraits for E.1.2 and E.2.1 critical points. For a symmetry in plots, we considered −1 ≤ x ≤ 1 interval, however, we should remember that x ∈ (0, 1].   (8). The left plot represents the model where the non-gravitational interaction Q is given by Equation (26). The right plot represents the model where the non-gravitational interaction Q is given by Equation (29).
Another cosmological model has been studied in this paper assuming the interaction term Q has the following form In this case, 3 critical points are physically acceptable (among 4 critical points) presented in Table 3. The critical point E.3.1 is a late time scaling solution since and represents the state of the universe, where the varying Chaplygin gas Equation (8) is phantom dark energy with This is a late time attractor when 0 ≤ n ≤ 5, 0 < b ≤ 2/5, 0 < α ≤ 1. On the other hand, the critical points E.3.2 and E.3.3 are physically reasonable when b = 0 and the critical points are unstable. The phase portrait indicates that in considered case the evolution of the universe started from one of the unstable states (either E.3.2 or E.3.3) eventually will evolve to the state described by E.3.1 solution. It is not hard to see that in case of E.3.3 the universe will be in a matter dominating state and the accelerated expansion is not possible. In considered 3 models of this subsection we observed, that appropriate late time scaling attractors describe the state of the universe with phantom varying Chaplygin gas Equation (8). Moreover, the solution of the cosmological coincidence problem exists. Therefore, it is important to learn the constraints on the free parameters to have better understanding of the models. To this we will come in the second part of this work. Table 3. Critical points corresponding to interacting varying Chaplygin gas, Equation (8), for the non-gravitational interaction Q given by Equation (32).

S.P.
x y Type of Stability unstable node or unstable focus

Varying Chaplygin Gas
In Sections 3.1.1-3.1.3 we considered cosmological models where varying Chaplygin gas was given by Equation (8). In the next three subsections we will consider cosmological models, where varying Chaplygin gas is given by Equation (9).
The model with the varying Chaplygin gas, Equation (9), when non-gravitational interaction is given by Equation (26) has 4 critical points (Table 4, where E.4.3 is not physically reasonable solution). The study shows that the considered model has only one conditional late time attractor corresponding to E.4.4 solution. Moreover, the study shows that for imposed constraints on the parameters of the model we will obtain the states of the universe with a quintessence varying Chaplygin gas, Equation (9). However, for some combinations of the values of the free parameters (from initial constraints imposed on the parameters) E.4.4 solution can describe the state of the universe where the varying Chaplygin gas, Equation (9), is a phantom. Observed dual possibility shows an interesting departure of this model from previously considered models. Moreover, we found that the late time attractor E.4.4 is scaling attractor. In particular, to simplify the analysis, let us demonstrate that E.4.4 is a scaling attractor when 0 ≤ b ≤ 1/10 and n = 1. It is easy to see, that in this case, the attractor E.4.4 is scaling because Moreover, considered constraints on the parameters provide the deceleration parameter q to be 0 < q ≤ −1. On the other hand, we can see that the varying Chaplygin gas, Equation (9), is a quintessence dark energy with .
The solution of the cosmological coincidence problem in form Equation (35) is different from the solutions presented in Section 3.1 and clearly demonstrates that the solution can be obtained due to the non-gravitational interaction. In summary, starting the evolution from one of the initial states described by E.4.1 and E.4.2, the universe eventually will end up on the state described by E.4.4. Table 4.
Critical points corresponding to interacting varying Chaplygin gas, Equation (9), for the non-gravitational interaction Q given by Equation (26).

S.P.
x y Type of Stability The study shows that the late time attractor E.5.2 can represent the states of the universe where the varying Chaplygin gas, Equation (9), either is a quintessence dark energy (Table 6) or a a phantom dark energy (Table 7) depends on the values of the parameters of the model (the general constraints on the parameters are the same as in previous cases). Table 6. Constraints on the model parameters for the interacting varying Chaplygin gas, Equation (9), with non-gravitational interaction term Q given by Equation (29). E.5.2 describes the state of the universe where the varying Chaplygin gas is a quintessence dark energy. Table 7. Constraints on the model parameters for the interacting varying Chaplygin gas, Equation (9), with non-gravitational interaction term Q given by Equation (29). E.5.2 describes the state of the universe where the varying Chaplygin gas, Equation (9), is a phantom dark energy.
The phase space portrait of this models shows that for imposed constraints, the evaluation of the universe will start from the state described by E.5.1 (unstable node) and will reach to the state one described by E.5.2 (stable node). On the other hand, E.5.3 is physically not reasonable solution, while E.5.2 is a scaling attractor with

Interaction
The 4 critical points presented in Table 8 describe the model with interacting varying Chaplygin gas, Equation (9), when the non-gravitational interaction is given by Equation (32). The study shows that only the critical point E.6.4 is late time attractor describing the state of the universe with Moreover, it is a scaling attractor and the solution of the cosmological coincidence problem reads as wherer = −3bb 2 (1 + α) 2 + (n − 3(1 + α)) 2 . Similar to the other two models considered in this subsection, E.6.4 solution can describe the state of the universe where interacting varying Chaplygin gas, Equation (9), is either a phantom dark energy or a quintessence dark energy. In particular, the constraints on the model parameters presented in Table 9 describe the universe where the varying Chaplygin gas, Equation (9), is a phantom dark energy. It should be mentioned that E.6.3 is physically not reasonable solution, while E.6.2 is a stable focus and does not describe the accelerated expansion. Definitely the phase space analysis reveals that the models are very interesting. However, as it is mentioned earlier we have to impose some empirical constraints on the free parameters in order to reduce the phase space size. Definitely it simplifies the discussion. However, we should wonder whether or not we assumed reasonable constraints for the free parameters. In order to understand this, in the second part of this paper, we will apply Bayesian Machine Learning to learn the constraints. The detailed analysis is presented in the next section. The general philosophy behind the approach can be found in [32,35,37,38]). Following earlier works we also use PyMC3 to perform the analysis and learning process. Table 8. Critical points corresponding to interacting varying Chaplygin gas, Equation (9), for the nongravitational interaction Q given by Equation (32), wherer = −3bb 2 (1 + α) 2 + (n − 3(1 + α)) 2 .

S.P.
x y Type of Stability stable node or stable focus Table 9. Constraints on the model parameters for the interacting varying Chaplygin gas, Equation (9), with the non-gravitational interaction term Q given by Equation (32). E.6.4 describes the state of the universe where the varying Chaplygin gas, Equation (9), is a phantom dark energy.

Constraints from Bayesian Machine Learning
In the first part of this paper we discussed the phase space analysis of the models. It appears possible to find all critical points analytically allowing to discuss conditional critical points and attractors. The last means that for imposed constraints some of the obtained solutions can serve as attractors and have specific type of stability, which can be changed when another set of constraints on the model free parameters will be imposed. A necessity somehow to constrain the possible parameter space and decrease the uncertainty we applied the Bayesian Machine Learning. We will omit all technical details behind this approach referring the readers to several recent papers [32,35,37,38]) for more details. The crucial aspect of this approach to be mentioned is the option connecting posterior and prior without a need to evaluate the likelihood. This allows us to depart from the real observational data concept and use simulated one. It is obvious that the method is good for forecasting purposes too. However, the validation of the learned result, as in the case of any other Machine Learning approach is required. In this early stage of analysis we used the background dynamics of each model to generate only the expansion rate data used for the learning process. Therefore available expansion rate data presented in Table 10 has been used during the result validating process. It should be mentioned that our interest to the expansion rate data is due to the estimation of H(z) from the cosmic chronometers. Since they are model independent estimations then there is a possibility to learn features of the models not depending on the features of the other models used to extract H(z) data. However, the analogues of other observations can be crafted/generated and used to learn the constraints, too. This and other interesting questions for this moment have not been studied. They have been left to be tackled in forthcoming papers. In the next two subsections we will discuss learned constraints and see whether or not the models can be used to solve the H 0 problem (see [15][16][17][18] for more details). We start with the case where the dark energy model is given by Equation (8). We took into account Equations (3), (6) and (7) to construct the background dynamics used in the generative process. We considered three forms of non-gravitational interaction given by Equations (26), (29) and (32) respectively to complete the field equations. Assuming that we have cold dark matter with P dm = 0 we learned the constraints on 7 free parameters in each case using generated expansion rate data. The learned constrained can be found in Table 11. During the learning process we imposed H 0 ∈ [64.0, 80 and Ω dm ∈ [0.1, 0.45] flat priors on free parameters. The results presented below are from the analysis based on 10 chains and in each chain, 10,000 "observational" data-sets from the models have been simulated. After a very detailed analysis of our results we found that: • When the interaction is given by Equation (26)  and Ω dm = 0.247 ± 0.014. The contour map is given in Figure 2, in red colour.
As we see the models cannot solve the H 0 tension problem. Moreover, we found a hint that the model with interaction given by Equation (32) should be ruled out due to very low Ω dm . Moreover, all of them are in huge tension with high redshift expansion rate data and non of them can explain the BOSS experiment results [94]. Another interesting result to be mentioned is about the learned constraints for free parameter n. We see that the learned best fit value of it can be negative for the model with Equation (26) interaction term. In general, during our analysis we learned that the expansion rate data will not put tight constraints on it. On the other hand, this gives a hint that considered constraints on the free parameters during the phase space analysis just allows partially cover the phase space of the models. However, we obtained all critical points analytically, therefore with future improved new tight constraints on the free parameters we can explore phase space regions more carefully. To summarize the results of the learning process we can say that there is a hint that all three models should be ruled out. One thing is clear that still additional analysis is needed for the final conclusion. Table 11. Best fit values and 1σ errors estimated for the Model given by Equation (8), when z ∈ [0, 2.5]. The results have been obtained from a Bayesian Machine Learning approach, where the interaction Q is given by Equation (26) (first row), Equation (29) (second row) and Equation (32) (third row), respectively. Moreover, when the interaction is given by Equation (26) we obtained Ω dm = 0.268 ± 0.009, while when the interaction is given by Equation (29) Ω dm = 0.263 ± 0.008 has been obtained. Finally, Ω dm = 0.247 ± 0.014 has been obtained when the interaction has been given by Equation (32). We  and Ω dm ∈ [0.1, 0.45] flat priors have been imposed during the generative process used to generate the "observational" data. The analysis is based on 10 chains and, in each chain, 10,000 "observational" data-sets from the model have been generated.

Varying Chaplygin Gas
Finally, we analyzed the model where the dark energy is given by Equation (9). Again, we took into account Equations (3), (6) and (7) to construct the background dynamics used in the generative process. In this case, also we considered three forms of non-gravitational interaction given by Equations (26), (29) and (32) respectively. Assuming that we have cold dark matter with P dm = 0 we learned the constraints on 7 free parameters in each case using generated expansion rate data. Similar to the first case, during the learning process we and Ω dm ∈ [0.1, 0.45] flat priors on free parameters. The results presented below are from the analysis based on 10 chains and in each chain, 10,000 "observational" data-sets from the models have been simulated. The learned constrained can be found in Table 12 and to simplify our discussion we present them below. In particular, using Bayesian Machine Learning based on generated expansion rate data we found that: • When the interaction is given by Equation (26)  Again we see the models cannot solve the H 0 tension problem. Moreover, we found a hint that the models should be ruled out due to very low Ω dm . Moreover, all of them are in huge tension with high redshift expansion rate data and non of them can explain the BOSS experiment results. In other words, the models considered in this subsection should be rejected. We would like to mention that the learning based on generated other data sets can induce a significant clarification on this issue, therefore an additional analysis is still required before the final rejection of the models. It has been left to be discussed in a forthcoming paper.  (9), when z ∈ [0, 2.5]. The results have been obtained from a Bayesian Machine Learning approach, where the interaction Q is given by Equation (26) (first row), Equation (29) (second row) and Equation (32) (third row), respectively. Moreover, when the interaction is given by Equation (26) we obtained Ω dm = 0.233 ± 0.012, while when the interaction is given by Equation (29) Ω dm = 0.229 ± 0.012 has been obtained. Finally, Ω dm = 0.215 ± 0.017 has been obtained when the interaction has been given by Equation (32).    and Ω dm ∈ [0.1, 0.45] flat priors have been imposed during the generative process used to generate the "observational" data. The analysis is based on 10 chains and, in each chain, 10,000 "observational" data-sets from the model have been generated.

Conclusions
In this paper, we considered six different cosmological scenarios. Three different forms of the interaction term Q have been used to model non-gravitational interaction supposedly existing between dark energy and cold dark matter. The Chaplyging gas (and its various modifications) is still among very actively considered dark energy gas/fluid models. Actually two modifications of it has been considered. In particular, with specific phenomenological assumptions about its free parameter it is possible to construct new interesting modifications. Following to this approach, in this work we considered two models where two parametrizations of B are taken into account. We performed the phase space analysis of all models. In all cases, late time scaling attractors have been found analytically and by imposing empirical constraints their nature have been explored. In particular, it has been found that for the first type of cosmological models the late time scaling attractors describe the state of the universe where varying Chaplygin gas has only a phantom nature. On the other hand, the study of the second class of models shows that for some values of the model free parameters the varying Chaplygin gas will be a phantom dark energy, while for some cases it will be a quintessence dark energy. An interesting result that has not been observed during the study of the first model. Moreover, in the second part of the paper we applied Bayesian Machine Learning approach to learn the constraints on the model parameters. It is a good way to understand the models and extract additional information justifying the reason specific phenomenological modifications should exist. On the other hand, it appears to be very useful for future studies pointing out the shape of the parameter space. In our analysis the background dynamics of each model has been used to generate the expansion rate data to be used in the learning process. We have a very brief discussion on this, referring the reader to a series of papers demonstrating how the Bayesian Machine Learning can be used to achieve very interesting and unique results. The results from the Bayesian Machine Learning based on generated expansion rate data revealed interesting results. In particular, the results of our analysis showed that it is very hard to say the models eventually should be rejected or not. However, we found a hint that they should be rejected since (1) they cannot solve the H 0 tension problem, and (2) there is a huge tension between learned and observational results at high redshifts. However, we need to take into account that the learned constraints indicate-the expansion rate data is not good to learn the constraints on the model free parameters. It can be that the learning process based on other generated data sets can do the job better and eventually it would be possible to learn tight constraints on the parameters. This among other problems has been left to be tackled in forthcoming papers. In this stage of the analysis the Bayesian Machine Learning indicates several interesting aspects of considered varying Chaplygin gas models. One of them is related to the forms of non-gravitational interaction. In particular, a recent analysis of [37] demonstrated that a deviation form cold dark matter paradigms can easily be confused with non-gravitational interaction between dark energy and dark matter. In this regard, we need to mention that a similar situation could be observed here too. In this case, if it is true, then we can significantly reduce the phenomenology and craft new cosmological models with varying Chaplygin gas, where the model rejection question will be more transparent and easily understandable. Our initial attempts to study the Chaplygin gas with Machine Learning approaches are very promising and should be extended, revealing whether or not it can be really used to unify dark energy and dark matter. Other interesting questions related to this issue we left to be tackled in the future papers. In general with future research we still need to understand whether or not the expansion rate data can capture the features of non-linear non-gravitational interaction between dark energy and dark matter properly. In addition, if it cannot, then how to be with the bias that can be introduced in the analysis where the expansion rate data has been used with other observational data sets and the features of the non-gravitational interaction have been captured?
In summary, the Bayesian Machine Learning reveals that something is not smooth with varying Chaplygin gas models we considered here. However, a final conclusion requires additional work which can provide a fresh look at (varying) Chaplygin gas cosmology and significantly reduce so far existing phenomenology.