Exact Solutions and Accelerating Universe in Modified Brans-Dicke Theories

Exact solutions are studied in the context of modified Brans-Dicke theory. The non-linearity of the modified Brans-Dicke field equations is treated with the Euler-Duarte-Moreira method of integrability of anharmonic oscillator equation. While some solutions show a forever accelerating nature, in some cases there is a signature flip in the evolution of deceleration parameter in recent past. Importance of these latter models are studied in the context of late-time acceleration of the universe. Constraints on the model parameters are obtained from Markov Chain Monte Carlo (MCMC) analysis using the Supernova distance modulus data, observational measurements of Hubble parameter, Baryon acoustic oscillation data, and the CMB Shift parameter data.


Introduction
Introduced back in 1915 [1], General Theory of Relativity (GR) has remained the most successful description of gravity till date. Based on an equivalence between gravitation and inertia, GR ensures that one can write down the geodesics on a spacetime from the metric structure [2]. The foundation principles of the theory has been well tested over the years through famous experiments, for instance, the Eötvös experiment [3], Michelson-Morley-type experiments [4,5], and the gravitational redshift experiments [6]. Moreover, recent observational evidences tend to confirm the theoretia e-mail: pm14ip011@iiserkol.ac.in b e-mail: soumya@cts.iitkgp.ernet.in cally predicted outcomes of GR as well, for example, the existence of black holes and gravity waves [7], promotes GR as a potentially rich store of possibilities even in the current context.
In spite of passing many theoretical and experimental tests successfully, there remains some unresolved issues hinting towards a possible modification of GR. An immediate motivation comes from the recently discovered late time accelerated expansion of the universe. To account for such a non-trivial acceleration, the simplest possible modification is to consider a cosmological constant Λ playing the role of 'dark energy', a fluid responsible for an effective negative pressure. However, this option is problematic. Defining an empty space as a collection of quantum fields and assuming that the zero-point fluctuations of such quantum fields contributes to Λ , the theoretically predicted vacuum energy scale overwhelmingly mismatches with the observed vacuum energy [8]. A very well-studied alternative is to consider a time-dependent scalar field acting as the generator of the non-trivial acceleration, for example, the quintessence models or the phantom scalar fields (for extensive details on scalar field cosmology, see for instance, Ratra and Peebles [9], Brookfield et. al. [10], Overduin and Cooperstock [11], Bento, Bertolami and Sen [12] and references therein). Another possibility is to treat the dark energy component as a geometric quantity, coming out of a modified Einstein-Hilbert action. For example, replacing the Ricci curvature R in the Einstein-Hilbert action by a general function f (R) produces the so-called f (R) theories. For extensive reviews on different modifications of gravity, we refer to the works of Nojiri and Odintsov [13], Capozziello and De Laurentis [14], Faraoni and Capozziello [15], Bamba and Odintsov [16], Nojiri, Odintsov and Oikonomou [17].
Brans-Dicke (BD) theory is a modification of gravity which was formulated in order to incorporate the Mach's Principle in GR [18]. Mach's principle simply demands that the local motion of a particle must be affected by a largescale matter distribution. This can be violated in GR, for example in the de-Sitter universe dominated by a cosmological constant where matter is entirely negligible [19]. In BD theory, a scalar field φ is included in the action which makes the gravitational coupling a function of coordinates, rather than a constant. A dimensionless parameter ω, called the BD parameter, governs the departure of the results obtained under weak field approximation of the theory from that of general relativity. To be consistent with the local astronomical tests, ω must have a very high value (ω > 500). It can be proved that in the limit ω → ∞, φ reaches a constant value ∼ 1 ω , making the field equations of BD theory equivalent to the corresponding GR equations [20]. Despite carrying this great advantage of giving back GR in some limit, it was eventually proved that the infinite ω limit fails for a traceless stress-energy distribution, by Banerjee and Sen [21], Faraoni [22].
However, BD theory remains well regarded as the prototype of a large class of theories of gravity called the scalartensor theories, where a non-minimal coupling exists between a scalar field and the curvature scalar. These theories receive significant attention cosmologically as they can successfully describe an early inflation of the universe [23]. BD theory predicts a period of extended inflation to tackle the graceful exit problem of the universe as described by Mathiazhagan and Johri [24], La and Steinhardt [25]. Moreover, the theory can successfully generate the late time accelerated expansion as well, for suitable values of the parameter ω without the need of any exotic matter field (see for instance, Banerjee and Pavon [26]). Introduction of an additional self-interaction potential of the BD scalar field forms a straightforward modification of the theory and is well-studied in literature, for example by Faraoni and Gunzig [27], Bertolami and Martins [28]. Banerjee and Pavon proved that the theory can produce a non-decelerated expansion in the presence of an additional minimally coupled scalar field [29]. Sen and Sen looked into the possibility of a late time acceleration in BD theory for some specific choice of an additional potential term [30]. Das and Banerjee showed the possibility of a late time accelerated expansion preceded by a decelerated expansion in the domain of the theory, considering a non-minimal coupling of matter sector and the BD scalar field [31]. Cosmological solutions in BD theory were recently discussed by Jarv, Kuusk, Saal and Vilson for both Einstein and Jordan frames [32]. Generalizing the theory by making ω a function of the scalar field is another interesting possibility, proposed by Bergmann [33], Wagoner [34] and Nordtvedt [35]. Indeed, non-minimally coupled scalar-tensor theories (for example, the works of Barker [36] and Schwinger [37]) become equivalent to a generalized BD theory for a suitable choice of ω(φ ). For a brief review on such modifications and their cosmological motivations we refer to the work of Van den Bergh [38]. For recent discussions on BD theory and more generally, scalar tensor theories and their cosmological implications we refer to the works of Fujii and Maeda [39], Sotiriou [40].
In the present work, we focus on exact solutions in modified BD theory. The specific modified BD actions we study has previously been studied in literature ( [29,31,41]). However, exact solutions have not really been considered. Exact solutions are never guaranteed mainly due to the high degree of non-linearity of the BD field equations, however, they remain an interesting avenue of research carrying a pedagogical importance. Taking a different approach from assuming a cosmic expansion history at the outset, we incorporate a purely mathematical property of general second order differential equations with variable coefficients that can be point transformed into an integrable form. The property is derived from the symmetry analysis of classical anharmonic oscillator equations by Duarte, Moreira, Euler and Steeb [42]; generalized by Euler, Steeb and Cyrus [43] and thereafter expanded by Euler [44], Harko, Lobo and Mak [45]. The exact solutions extracted using this property give accelerating cosmological solutions. While some of the solutions give forever accelerating solutions, some examples indeed show a signature flip of deceleration parameter in recent past, hinting that such models can work well to model a late time acceleration. For relevant models we perform a Markov Chain Monte Carlo (MCMC) analysis using the supernova distance modulus data, observational measurements of Hubble parameter, baryon acoustic oscillation data and the CMB Shift parameter data to study the constraints on the model parameters. We reconstruct the cosmological quantities for the best fit values of the model parameters.
We consider two different modifications of the standard BD action. The first modifiction involves a non-minimal coupling of the matter sector with the BD scalar field. A similar action finds it's existence in the low energy limit of string theory or the dilaton gravity. Under such a setup, the BD scalar field behaves as a chameleon scalar field (for an idea and some examples of chameleon fields and their relevance in cosmology, we refer to the works of Khoury and Weltman [46,47], Mota and Barrow [48,49], Das, Corasaniti and Khoury [50], Mota and Shaw [51,52]). Clifton and Barrow [41] studied the cosmological solutions and their relevance for a similar setup of BD theory. Das and Banerjee [31] studied accelerating solutions in this setup assuming cosmic expansion history at the outset and discussed the evolution of deceleration parameter for different epochs. In the present work, instead of using any apriori ansatz over the scale factor or the BD scalar field, we solve the modified BD field equations for a general power-law non-minimal coupling between matter and scalar field. The second modification involves a minimally coupled self-interacting scalar field serving as a quintessence matter within the standard action of BD theory. A quintessence matter is added to enlarge the scope of the theory since this matter field in known to describe accelerated expansion under the scope of standard GR itself, confirmed by the observational data from the luminosity-redshift relation of type Ia supernovae [53]. Amongst many scalar field models considered as quintessence matter in literature, the model proposed by Zlatev, Wang and Steinhardt [54] serves particular importance, where a tracker field slowly rolls down the potential. The role of such a scalar field in BD theory was studied by Banerjee and Pavon [29] for some particular choices of the self-interaction potential. We study exact solutions for different self-interaction potentials, a simple power law and the other being a combination of power law terms. For a self-interaction potential of the form V (φ ) ∼ φ δ 1 + φ δ 2 , the integrability property of the scalar field equation produces some interesting exact solutions hinting at a late-time accelerated expansion.
In section 2, we briefly outline the mathematics of point transformation and direct integration of a general anharmonic oscillator equation. In section 3 we present exact solutions describing a forever accelerationg cosmology for some modifications of standard BD theory. Section 4 contains exact solutions describing a late time acceleration of the universe for a Quintessence plus BD scalar field setup. Solutions present in Section 4 also closely resemble the observational data as checked in Section 5 through parameter estimation by statistical analysis and the analysis of the evolution of cosmological parameters. We conclude the manuscript in section 6.

Integrability of Anharmonic Oscillator Equations
An anharmonic oscillator can be written as a second order differential equation having variable coefficients.
Here, f 1 , f 2 and f 3 are variable coefficients which are functions of t and n is a constant. This equation can be written in an integrable form by a pair of point transformations, if and only if n / ∈ {−3, −1, 0, 1}. Additionally, the coefficients must satisfy the differential condition The point transformations are to be defined as where C is a constant.
The scope of this approach in gravitational physics has been studied only very recently, for example, in massive scalar field collapse [55,56] and Scalar-Einstein-Gauss-Bonnet gravity [57], cosmological reconstruction of modified theories of gravity [58,59].

Brans Dicke Scalar Field as a Chameleon (Model I)
In this section we work with a setup where the BD scalar field φ is non-minimally coupled to the matter distribution through an arbitrary function f (φ ). The relevant action is written as where ω is the Brans Dicke parameter. R is the standard Ricci scalar and L m defines the matter distribution which we assume to be pressureless dust. One may note that for f (φ ) = 1, one gets back the usual BD action. We define the metric for a homogeneous and isotropic, spatially flat universe as where a(t) is the scale factor of the universe. Variation of the action with respect to the metric gives the field equations as, 2ä 4 A dot denotes derivative with respect to cosmic time t. ρ m denotes density of the fluid distribution. The condition for conservation of energy-momentum distribution for a dust fluid can be written aṡ Varying the action with respect to the BD scalar field φ one can also writë From equation (9) one arrives at the evolution of matter density as a function of the scale factor, written as where, ρ 0 is a constant of integration. One can note from equation (11) that due to the presence of the non-minimal coupling function f (φ ), there is a departure from usual matter conservation equation. Equation (11) and (10) together produces the following differential equation as We study the equation (12) for a power law coupling function, i.e., f (φ ) ∼ φ δ . Under such an ansatz for f (φ ) equation (12) falls within the large class of equations that can be identified as a classical anharmonic oscillator equation. These equations can be point-transformed into an integrable form. We now consider a power law coupling such that f (φ ) can be written as f (φ ) = φ m . With this the equation (12) can be written as We write n as − m 2 . Therefore, One may note that this equation indeed is a particular case of anharmonic oscillator with the coefficients identified as, Therefore one can transform the above equation into an integrable form and integrate for the BD scalar field. Moreover, the integrability criterion (2) is expected to reduce into a solvable equation of the cosmological scale factor. For the analysis to stand, n can not be −3, −1, 0, 1. This implies some restrictions over the choice of m as m = 6, 2, 0, −2 i.e., f (φ ) = φ 6 , φ 2 , 1, φ −2 . From equation (2), we find an exact solution for the scale factor which goes as, From the expression of the scale factor (18) it is evident that the nature of the exponent governs the evolution. Now using the condition (3) we point transform the scalar field evolution equation (12) and solve for the BD field. The solution can also be written from (12) simply by using the solution for scale factor as in (18) and is given by, where, φ 0 is to be determined from the consistency check of the theorem. Amongst the equations (7), (8), (9) and (10), only three are independent equations as the fourth can always be derived using the Bianchi identity. Therefore an exact solution coming out from (9) and (10) is a consistent solution as long as it satisfies any one of the field equations. From Eq. (8) the consistency criterion can be written as For a numerical analysis, the value of ρ 0 is obtained from the relation of matter density parameter, Ω m0 . We define Ω m0 = 8πGρ 0 3H 0 2 , where we use the value of Ω m0 ≈ 0.3 [60] and H 0 = 72.8 km Mpc −1 sec −1 [61], being consistent with recent observational data. The value of ω considered is nearly 60, 000 for the results to be consistent with the local astronomical tests. We choose a suitable value of n following Eq. (20) such that this domain of ω is strictly maintained. Fig. 1 shows the logarithmic variation of φ (top graph) as a function of z scaled with φ 0 . The graph in the middle gives the variation of φ w.r.t. z at low redshifts. The lowermost graph gives the variation of φ as a function of the Brans Dicke parameter ω.
We note that the deceleration parameter for Model I is a constant for all z, given by 2n+3 n+3 . Therefore, depending on the respective n value we will get either a constant accelerating or decelerating universe. Thus, Model I becomes redundant as far as present observational cosmology is concerned. However, it still serves as a simple toy model based on simple power-law potential assumption, even more so because the scale factors they produce are very simple. Moreover, a power law evolution of scale factor a(t) ∼ t α may have more interesting role in early universe cosmology.

Brans Dicke with a Quintessence and Simple Power Law Potential (Model II)
In this section we study the standard BD action alongwith a self-interacting scalar field playing the role of Quintessence matter. The corresponding action is given by, φ defines the BD Scalar Field, ω is the BD parameter and R is the Ricci scalar. L m is the matter distribution and in the present case, is defined by a combination of a perfect fluid and a spatially homogeneous scalar field ψ with a selfinteraction potential.
The field equations, where the metric is given by equation (6) are, ρ m and p m denotes the density and the pressure for matter sector. ρ ψ and p ψ denotes the density and pressure for the quintessence field. They are written as a function of ψ as V = V (ψ) defines the self-interaction potential for the Quintessence scalar field. Varying the action with respect to the scalar field ψ one gets the wave equation corresponding to the scalar field ψ as, Moreover, varying the action with respect to the BD scalar field, one gets the wave equation corresponding to φ as The matter conservation equation for the fluid giveṡ It is not unphysical to assume that in the current era p m can be put to zero, since the universe is dominantly filled with cold matter. Thus the matter conservation equation gives where ρ 0 is a constant of integration. We note here that the evolution equation (25) for ψ is a simple case of classical anharmonic oscillator for suitable choices of the selfinteraction potentials, as is evident from the form of Eq. (1).
Therefore we intend to extract as much information as possible from the above set of equations by virtue of the integrability analysis. Eq. (25) will provide us with the exact solution for the scale factor and the Quintessence field ψ for which the equation will be integrable. These solutions will in turn be used to determine ρ m from Eq. (28). Finally, using the solutions for a(t), ψ(t) and ρ(t) in Eq. (26), we expect to have an idea of the time evolution of the BD scalar field φ . In the following subsections, we study the system of equations for some relevant self-interaction potentials of the quintessence scalar field.
The first example taken is for a self-interaction potential power law in ψ, given by Therefore the evolution equation for the Quintessence Eq. (25) becomes Straightaway, this can be identified as a simple case of anharmonic oscillator equation. Following the detailed method as in section II, the solutions for a(t) and ψ(t) can be written as provided n / ∈ {−3, −1, 0, 1}. ξ can be evaluated from a consistency check, written as From the expression of the scale factor (31) we note that for flat cosmology one must have n > 0, whereas a negative n very close to zero gives open cosmologies. For −1 < n < 0, a late-time acceleration can be realized. Using the exact solutions, and the expression of density (using Eq. (31) in Eq. (28)) one gets a second order differential equation for the BD scalar field φ (t) as, φ + (n + 3) Eq. (34) is solved numerically to study the evolution of φ . The evolution as a function of redshift z is shown in Fig  3, for different values of n. The value of ρ 0 is obtained from the relation of matter density parameter, Ω m0 just as in case for model I. The BD parameter ω is considered to be 60, 000 to meet the requirement of local astronomical tests.
We note that the value of the deceleration parameter for Model II is a constant for all z. For Model II it is 2n n+3 . So, depending on the respective n value we will get either a constant accelerating or decelerating universe. Thus, Model II become redundant as far as present observational cosmology is concerned. However, it can serve as a simple toy model based on simple power-law potential assumptions, even more so because the scale factors they produce are very simple. Moreover, a power law evolution of scale factor a(t) ∼ t α may have more interesting role in early universe cosmology.

Brans Dicke with a Quintessence and Higgs Potential (Model III)
The Higgs interaction potential is defined as which is basically a combination of two simple power law self-interaction terms, one a quadratic and another a quartic in φ . This serves an additional purpose as it expands the scope of the integrability analysis. A quadratic potential straightaway does not fall within the scope of the theorem as for a simple φ n term in the equation, n / ∈ {−3, −1, 0, 1}. For the Higgs potential, Eq. (25) becomes Analyzing the equation in a similar fashion as in section II, we find the exact solution for the scale factor as We note that the general solution of scale factor for a Λ CDM cosmology goes as a(t) ∼ (Ae αt + Be −αt ) 2 3 where A, B, and α are independent constants. This may somewhat resemble equation Eq. (37) except for the exponent value being 1 2 instead of 2 3 and that we have two independent parameters λ 0 and µ instead of three.
The general solution for ψ(t) in this case can be written in terms of Gauss Hypergeometric function as follows The solution for ψ is non-trivial to work with without any further simplification. As we are interested in a late time cosmological dynaimics, we note that in the limit t → ∞, a hypergeometric function can be written as a simple series expansion. Morepver, for large t the term containing e √ 2µt dominates over the term of e − √ 2µt in the expression of scale factor as in Eq. (37). Thus, we can write, where, . Using the approximate solution for ψ in Eq. (26) one gets a differential equation for the BD Scalar field φ which is solved numerically. The solution of the BD scalar field φ , scaled with φ 0 is plotted as a function of redshift in Fig. 5. The value of BD parameter ω is again taken to be 60, 000. Moreover, the evolution of the Quintessence scalar field ψ with redshift is plotted in Fig. 4. The nonminimally coupled scalar field theories allow for a variation of the strength of the gravitational interaction as well. As 1 φ behaves as the effective Newtonian gravitational constant G [20], one can write, The values of model parameters have been taken from Table 1. The initial conditions for solving the second order differential equation of φ for the Models II and III are chosen obeying the constraint on G according to Eq. (41).
Comparing Fig. 5 and Fig. 4 we comment that at present epoch z ∼ 0 the BD scalar field φ almost becomes a constant, while the quintessence scalar field ψ has a sharp increasing behavior. At an early epoch z >> 0 this behavior simply reverses. One can interpret that within the domain of early universe cosmology, the BD scalar field plays a much more dominating role over the quintessence field. During present epoch, the role reverses and the quintessence field dominates the accelerated expansion. However, detailed comments on such aspects require further analysis and shall be addressed elsewhere.
Using the solutions as in Eq. (39) one can rewrite ψ anḋ ψ as, Thus, the expression for quintessence energy density ρ ψ becomes In BD theory the scalar field is proportional to inverse of the Gravitational constant. Therefore one can write Using equations (42), (44) and (45) in the Friedmann equation (21), we get the expression for reduced Hubble parameter as, where, Figure 5 depicts that in a late time scenario, φ behaves as a constant (φ 0 ), which implies ε can be considered 0 in present epoch. So, the modified expression for the reduced Hubble parameter becomes, Using the constraint equation such that, in the present epoch at z = 0, a = a 0 = 1, h = 1, we are finally left with, Note that, if µ is chosen as unity, equation (48) exactly mimics a ΛCDM model where, the constraint term V 0 3H 2 0 φ 0 behaves as Cosmological Constant.

Brans Dicke with a Quintessence and a general combination of Power Law Potentials (Model IV )
We also present an example where the self-interaction potential of the quintessence scalar field is a simple combination of power functions ∼ ψ 2 + ψ δ , not restricting δ to 4 only. Similar to the higgs potential, this allows us to observe the role of a quadratic term in the potential. Moreover, this allows one more parameter in the solutions, which is the exponent, allowing wider variety of solutions and possibilities. The potential can be written as using which one can write the scalar field evolution Eq. (25) as On a comparison, this gives a class of the general oscillator equations (1) for different n. Using the aforementioned method of integrability we arrive at the evolution equation for scale factor a(t) and Quintessence field ψ(t) written as Carrying out a similar treatent as in Section 3.1 one may find that the qualitative behavior of the Quintessence scalar field ψ and the BD scalar field φ are similar to that of Model III, apart from some scaling.

Parameter Estimation
The Models in section 3, i.e., models in III and IV , involving a combination of two potential functions, gives the deceleration parameter as a function of redshift, which might show a possibility for signature flip in the expression of q. Therefore, it is prescribed to estimate the model parameters and study the evolution of the different cosmological quantities extensively for Model III and IV . In this section, we estimate the model parameters and study the confidence contours on the parameter space, the marginalized likelihood function using four different data sets; the Supernova distance modulus data (SNe), observational measurements of Hubble parameter (OHD), Baryon Acoustic Oscillation (BAO) data and the Cosmic Microwave Background (CMB) data.
We first write the cosmological quantities in a dimensionless way from the expression of scale factor a(t) and it's time derivatives, using the redshift as the argument instead of cosmic time t. Redshift z is a dimensionless observational quantity defined as, a 0 being the present value of the scale factor. The Hubble parameter H(t) =ȧ a can also be written as a function of the redshift z as We can estimate the Hubble parameter from equation (37) and (51) Table 1 The parameter values and the associated 1σ uncertainty of the parameters of Model III, obtained from the analysis with different combinations of the data sets.
We note that the contribution of the Higgs field towards the total energy density of the Universe resembles a radiation dominated Universe in the early epoch. For parameter estimation, we use the distance modulus measurements of type Ia supernova from the Joint Light-Curve Analysis following the work of Betoule et. al. [62], who studied cosmological constraints from the SN-Ia observations of SDSS-II and SNLS collaborations. We incorporate the estimation of the Hubble parameter as a function of redshift [63], alongwith the measurement of Hubble parameter from Lyman-α forest at redshift z = 2.34 by Delubac et. al. [64] and measurement of H 0 from Planck [60]. For BAO data, three independent measurements of r s (z d ) D v (z BAO ) are used. r s (z d ) gives the sound horizon at photon drag epoch (z d ) and D v is the dilation scale at the redshift of BAO measurement. Three measurements are for three different values of redshift, for instance, from 6dF Galaxy Survey at z = 0.106 [65], from Baryon Oscillation Spectroscopic Survey (BOSS) at z = 0.32 (BOSS LOWZ) and at z = 0.57 (BOSS CMASS) [66]. The BAO measurements have been scaled by the acoustic scale (l A ) estimated from Planck [60]. The CMB shift parameter is related to the position of the first acoustic peak in power spectrum of the temperature anisotropy of the CMB radiation. Value of the parameter is estimated from the CMB data along with some assumption about the model of background cosmology, as estimated from Planck data [60].
Uncertainty of the parameters are estimated by 'Markov Chain Monte Carlo' (MCMC) method with the assumption of a uniform prior distribution. In the present analysis, we have adopted a Python implementation of the ensemble sampler for MCMC, the 'emcee', introduced by Foreman-Mackey et al. [67]. Fig. 6 shows the confidence contours on the parameter space and the marginalized likelihood function of Model III and IV obtained from the combined analysis with different datasets. It clearly depicts that, the parameters h 0 and λ 0 for Model III have a positive correlation between themselves. The best-fit values of the model parameters have been estimated for both cases, and the associated 1σ uncertainty, obtained for different combinations of the data sets in case of Model III and IV are represented in Table 1 and 2 respectively.

Evolution of Cosmological Parameters
Depending on the best fit choice of model parameters, the functional form of the cosmological quantities can be plotted as a function of redshift. Figs. 7 and 8 show that the plots Table 2 The parameter values and the associated 1σ uncertainty of the parameters of Model IV , obtained from the analysis with different combinations of the data sets. The kinematic quantities related to the expansion of our universe play a vital role in cosmology. They basically are connected to the second and third order derivatives of the scale factor. For instance, the deceleration parameter is written as, The jerk parameter can be written as The effective equation of state parameter w e f f , represented by a ratio of the total pressure contribution to the total  Fig. 9 shows that for Model III, there is a transition in the signature of q(z) from a decelerated phase to an accelerated phase of expansion. Moreover, the transition redshift z t < 1 is consistent with direct observational results [68,69]. Thus, it can be said that the Higgs' interaction potential model is consistent with the observed evolution of q(z). But in case of Model IV no such transition behavior could be seen. So, Model IV is inconsistent with the observed evolution of q, and hence can be ruled out. This situation probably arises because of the fact that in Model IV the coefficients for combined potential functions is considered as unity, which hints towards a uncompensated reduction in the necessary model parameters.
The deceleration parameter is expected to show a nontrivial evolution with respect to redshift as recent observational cosmology suggests. Therefore it is very important to investigate the next order derivative of the scale factor or the jerk parameter, whose value is unity for Λ CDM model. Model III shows a non-trivial evolution of jerk parameter as a function of redshift, shown in Fig. 10. The present value of jerk parameter obtained from Model III remain in between best fit 1 2 Fig. 9 The deceleration parameter q(z) plot for Model III using the best fit values and the associated 1σ , 2σ confidence regions from the combined analysis.   [60] 0.9 to 1.3 at 1σ level, well in agreement with the requirement of observational evidences. This also confirms the recent prediction that jerk parameter is in general expected to be in close proximity with the corresponding Λ CDM value [70].
The graph on the top of Fig. 11 shows the evolution of w e f f for Model III. At the present epoch, z = 0, w e f f is around −1 but slightly greater. While this may indicate a 'phantom' nature, it also strongly indicates a Λ CDM behavior depending on different model parameter values. With increase in z, w e f f gradually rises and becomes a constant at ∼ 0.3 which resembles a radiation dominated universe at high redshift (bottom graph of Fig. 11).
In Table 3, we compare different values cosmographical parameters for the Brans-Dicke + Higgs (Model III) setup with corresponding values of the parameters from Λ CDM and observations.

Conclusion
It can not be denied that Brans-Dicke theory might have renounced it's original appeal a little bit as long as it can not reproduce GR at some limit of the BD parameter ω. However, this does not diminish the stature of the theory as a prototype of scalar-tensor theories of gravity. Moreover, the theory has the elegance of describing both the inflationary universe and the present accelerated expansion of the universe without any need of dissipative processes or an exotic fluid. One often considers generalizations or modifications of the theory in a hope of tackling the shortcomings and making it a 'better' theory of gravity.
The present work deals with accelerating solutions in modifed BD theory. Keeping in mind the non-linearity of the system of field equations, a mathematical method of treating an anharmonic oscillator equation system is incorporated. The BD scalar field evolution equation is treated with an Euler-Duarte-Moreira method of integration. The advantage of this method is that, it helps one to solve the system of equations without any apriori assumption on the scale factor or the scalar field from the cosmic history, but rather, the restriction over the choice of the functions come from a purely mathematical property.
The solutions for three different models are discussed. For model III, a parameter estimation is carried out and the confidence contours on the parameter space are studied using four data sets, namely, SNe, OHD, BAO and the CMB data. The cosmological quantities, for instance the effective equation of state parameter ω e f f , the deceleration and the jerk parameters are plotted as a function of redshift for the best fit choices of the model parameters. For Model III the evolution of ω e f f closely resembles a Λ CDM behavior around z ∼ 0 and a matter dominated universe at high redshift. Model III also shows a transition in the signature of q(z) and the transition redshift is consistent with direct observational results. Models I and II has zero value of the jerk parameter whereas Model III shows a non-trivial evolution of j(z) with the present value of jerk parameter remaining in between 0.9 to 1.3 in the present epoch. We can therefore note that the Brans Dicke with quintessence or a 'BDQ' setup with a Higgs interaction potential case is well consistent with the observed evolution of cosmological quantities, as compared to the other significant options. A general combination of power functions of the quintessence field is also considered as potential. In such a case, the scale factor is seen to describe an accelerated expansion. For some choices of the potential, the hubble parameter also follows closely the observational data-points. However, not all combinations produce a result consistent with observational prediction as they fail to describe the signature flip in the evolution of deceleration parameter.
The simplification of the BD field equations and a subsequent extraction of an exact solution is extremely nontrivial. This difficulty is often by-passed by studying simplified systems who do not fail to describe the physics involved. In that sense, the present work also presents some special cases, as the solutions come under a specific assumption over the integrabilty of one particular equation only. However, the cosmological solutions found through the present method seem to describe the accelerated expansion of the universe quite well, at least for some specific cases. The solutions are simple and easy to work with for further allied investigations. Moreover, the simplicity provides one with the opportunity to study the evolution of the BD scalar field in a general manner.
It must also be mentioned that the issue regarding the value of BD parameter ω remains to be solved such that the theory can solve cosmological requirements along with satisfying the lower limit on the BD parameter, imposed by the solar system experiments. To conclude we note that different single or multiple scalar field models remain extremely popular even in the present context candidates to fill in for the fluid responsible for the late-time acceleration. However, it is always better to have a complete mathematical theory providing a unified description of the whole expansion history of the universe, from an early inflationary epoch to a late time cosmic acceleration, and beyond. It was quite extensively discussed by Elizalde et. al. [72] that given a certain universe expansion history, one can reconstruct a wide class of minimally or non-minimally coupled scalar field theories presenting a number of explicit examples which show a unified description of the inflationary era and a late-time cosmic acceleration epoch. While the models discussed in the present work do not describe a unified cosmic history in that sense, it shows potential in producing interesting exact solutions describing atleast some patches of our universe, namely the late-time era. Some simple power law solutions also have the potential to describe early inflation of the universe. Some modification of the starting action of our models may solve these problems, for instance, considering ω a function of φ which gives a Nordtvedt type theory [35]. We note therefore, in conclusion that the application of the anharmonic oscillator treatment for a varying ω theory is perhaps the next prescribed step. This may produce novel solutions under the scope of the theory bridging an early inflation, a decelerating radiation and a late time accelerated expansion without violating the astronomical constraints.