1 Introduction

General theory of relativity is widely regarded to be the best theory of gravity as over the years it has passed many observational tests with flying colors. Significant observational aspects such as the perihelion precission of Mercury’s orbit, gravitational lensing, redshift in the light spectrum from extragalactic objects are well documented in the framework of general relativity. However, general relativity can pose considerable intrigue in certain aspects such as the possibility of an extremely strong gravitational field imploding without limit and ending up in a region where the density of matter and the strength of the gravitational field can in principle, become infinite. Such a region is called a spacetime singularity. Aspects of gravitational collapse and the formation of a spacetime singularity forms an integral part of gravitational physics today.

In general, a continual gravitational collapse occurs whenever a massive astronomical body, upon exhausting it’s nucleur fuel supply, fails to support itself against the force of gravity. For a simple enough configuration of collapsing matter, a horizon generally develops prior to the formation of the singularity, thereby enveloping the singularity producing a black hole end-state. The first analytic model of such an unhindered contraction of an idealized star ending up in a black hole was given by Oppenheimer and Snyder [1] and independently by Dutt [2, 3] and these serve as the paradigm of gravitational collapse today. However, whether or not every sufficiently massive star undergoing gravitational collapse ultimately ends up in a black hole, remains an intriguing question even in the current context. In this regard, Penrose proposed [4] the cosmic censorship hypothesis (CCH), which roughly states that gravitational collapse of physically reasonable matters with generic regular initial data will always end up in a covered singularity.

The CCH, however, remains one of the most thought-provoking problems in gravitational physics today. There is no general proof of the CCH as yet which applies under all conditions. Moreover, there are many counterexamples in general relativity, where it is shown that the singularities formed in a collapse of reasonable distribution of matter can stay exposed, giving rise to the concept of a naked, i.e., an observable singularity (for a brief summary of such examples, we refer to [5] and references therein). This much is realized that the nature of the final outcome of a collapse depends on the initial configuration from which the collapse evolves, and the allowed dynamical evolutions in the spacetime, as permitted by the non-linear field equations of gravity. For recent works and different aspects regarding the dynamics of a continued gravitational collapse we refer to the summary by Joshi [6, 7].

The singularity or a spacetime region with infinitely high curvature can be realized only during the final stage of gravitational collapse, where the collapsing body almost reaches a zero proper volume and known laws of physics are expected to break down. A possibility of naked singularity implies that one may have a chance to observe indication of quantum effects of gravity. Amongst the many quantum theories of gravity proposed, superstring/M-theory is a promising candidate, which motivates the presence of a higher-dimensional spacetime [8,9,10,11,12,13,14,15,16,17]. During the final stages of the stellar evolution, where the curvature of the central high-density region is very high, the effects of extra dimensions can perhaps play a crucial role. From such a perspective, higher-dimensional gravitational collapse models are studied in general relativity (we refer to the works of Banerjee, Debnath and Chakraborty [18], Patil [19], Goswami and Joshi [20, 21] in this regard). Given that the entire aspects of superstring/M-theory are not understood completely so far, taking their effects perturbatively into classical gravity is one possible approach to study higher curvature effects. The Gauss-Bonnet (GB) term, \(G = R^2 - 4R_{\mu \nu }R^{\mu \nu } + R_{\mu \nu \alpha \beta }R^{\mu \nu \alpha \beta }\) in the standard Einstein-Hilbert Lagrangian is the higher curvature correction to general relativity which finds it’s motivation from heterotic superstring theory [22, 23]. Such a theory is called the Einstein-Gauss-Bonnet gravity.

The additional elegance of including the GB term is that this is a Lovelock scalar and if included linearly in the action, the field equations include no higher than second order partial derivatives of the metric tensor (unlike f(R) gravity whose field equations are fourth order in metric components). In a four dimensional spacetime the Gauss-Bonnet term does not modify the field equations. However, if the non-minimal coupling of a scalar field with the GB term is considered, the dynamical equations are quite different from the standard field equations and the influence of the GB term in a four dimensional universe is effective.

Recently there is an increasing interest in Gauss-Bonnet theory with a non-minimally coupled scalar field to suffice for the possible candidature of the late time acceleration of the universe [24,25,26,27,28]. From such a perspective, spherically symmetric solutions has been studied by Boulware and Deser[29, 30] and Gurses [31]. It has also been discussed that the effective action including correction terms of higher order in the curvature can perhaps play a significant role in the dynamics of early universe by Zwiebach[32], Zumino[33], Boulware and Deser[29, 30]. Questions regarding gravitational instability, cosmological perturbation were also considered by Kawai, Sakagami and Soda[34, 35], Kawai and Soda[36]. Observational restrictions over different cosmological aspects of the scalar field coupled Einstein Gauss Bonnet gravity were investigated by Guo and Schwarz[37], Koh et. al.[38]. Spherically symmetric collapsing solutions of this theory have also gained some interest quite recently. For instance, Maeda presented a model of \(n \ge 5\) dimensional spherically symmetric gravitational collapse of a null dust fluid in Einstein-Gauss-Bonnet gravity [39] and illustrated the possibility of a formation of massive naked singularity in higher dimensions. He comparitively analyzed the results with the general relativistic cases as well [40], which serves as a higher order generalization of the Misner-Sharp formalism of the four-dimensional spherically symmetric spacetime with a perfect fluid in general relativity. Hamiltonian Formulation of spherically symmetric scalar field collapse in Gauss-Bonnet gravity was studied in detail by Taves, Leonard, Kunstatter and Mann [41]. They also proved that such a formulation can readily be generalized to other matter fields minimally coupled to gravity. Apart from their role in cosmology, the role of scalar fields in a gravitational collapse is worthy of attention. It is indeed important to investigate if the CCH necessarily holds or violates in the collapse of fundamental matter fields. Moreover, a scalar field alongwith an interaction potential is known to mimic different kind of reasonable distribution of matter as discussed by Goncalves and Moss [42]. Numerical simulations of the spacetime evolution of a massless scalar field minimally coupled to gravitational field studied by Choptuik [43], Brady [44] and Gundlach [45, 46] hint interesting possibilities, such as, the critical behavior observed around the threshold of black hole formation. There is a self-similar solution that sits at the black hole threshold, dubbed a critical solution, and also a very interesting mass scaling law for the formed black hole end-state. Motivated by this, self similar collapsing scenario of a massive scalar field was analytically studied by Banerjee and Chakrabarti [47] very recently. Deppe, Taves, Leonard, Kunstatter and Mann presented a numerical analysis in generalized flat slice co-ordinates of self-gravitating massless scalar field collapse in five and six dimensional Einstein Gauss Bonnet gravity near the threshold of black hole formation [48]. Effects of higher order curvature corrections to Einstein’s Gravity on the critical phenomenon near the black hole threshold, were investigated by Golod and Piran [49]. In general relativity also, possibilities and scope of scalar field collapse have been analytically studied extensively [50,51,52,53]. Finally, the fact that the distribution of the dark energy component or the fluid still remains unknown, naturally inspires a continuing study of scalar field collapse under increasingly generalized setup (preferably alongwith a fluid distribution) towards a better understanding of the possible clustering of dark energy.

In this work, we study aspects of a scalar field collapse in a Scalar-Einstein-Gauss-Bonnet gravity, where the self-interacting scalar field \(\phi \) is non-minimally coupled to the GB term. Very recently Banerjee and Paul [54] studied such a scalar field collapse where the coupling term was proportional to \(e^{2\phi }\). In the present case the coupling term is proportional to \(\phi ^2\), i.e., the coupling is quadratic in \(\phi \). Quite recently, Doneva and Yazadjiev showed that in a very similar setup of Scalar Einstein Gauss Bonnet theory (the conditons they imposed on the coupling function \(f(\phi )\) are \(f'(\phi = 0) = 0\) and \(b^2 = f''(\phi = 0) > 0\)), there exists new black hole solutions which are formed by spontaneous scalarization of the Schwarzaschild black holes in the extreme curvature regime and below a certain mass, the Schwarzschild solution becomes unstable and new branch of solutions with nontrivial scalar field bifurcate from the Schwarzschild one [55]. They also proved the existence of neutron stars in a class of extended scalar-tensor Gauss-Bonnet theories for which the neutron star solutions are formed via spontaneous scalarization of the general relativistic neutron stars [56]. Very recently, the spontaneous scalarization of black holes and compact stars from such a Gauss-Bonnet coupling has been investigated and dubbed as the Quadratic Scalar-Gauss-Bonnet gravity by Silva et. al. [57]. In this context, existence of regular black-hole solutions with scalar hair in the Einstein-scalar-Gauss-Bonnet theory was investigated by Antoniou, Bakopoulos and Kanti [58, 59], with a general coupling function between the scalar field and the quadratic Gauss-Bonnet term which highlighted the limitations of the existing no-hair theorems. In a recent study Kanti, Gannouji and Dadhich have addressed the importance of such a coupling from a cosmological purview and proved by some simple analytical calculation that a quadratic coupling function, although a special choice, allows for inflationary, de Sitter-type solutions to emerge [60].

The inclusion of the Gauss-Bonnet terms make the dynamical field equations even more non-linear. We study a spatially homogeneous model where the energy momentum tensor is contributed by the self-interacting scalar field as well as a perfect fluid. To track down the system of equatons, we use the method of invertible point transformations and integrability of anharmonic oscillator equations; an approach which has been quite useful recently, in the study of minimally coupled massive scalar field collapse by Chakrabarti and Banerjee [47, 61].

The paper is organized as follows. In Sect. 2, we introduce the action and basic field equations. Section 3 contains the method of finding the exact solution and in Sect. 4 we study the evolution of the scale factor for different initial data. The time evolution of the scalar field is studied in Sect. 5. Evolution of the curvature scalars and the strong energy condition is studied in Sect. 6. The physical nature of the singularity is addressed in 7. We complete the model by matching the solution with an exterior Vaidya solution in Sect. 8 and conclude in Sect.  9.

2 Action and basic equations

The relevant action for a four dimensional action containing the Einstein-Hilbert part, massive scalar field and the Gauss-Bonnet term coupled with the scalar field. The corresponding action is given by

$$\begin{aligned}&S = \int d^4x \sqrt{-g} \bigg [R/(2\kappa ^2)-(1/2)g^{\mu \nu }\partial _{\mu }\phi \partial _{\nu }\phi -V(\phi )\nonumber \\&\quad - \xi (\phi )G \bigg ], \end{aligned}$$
(1)

where R is the Ricci scalar, \(1/(2\kappa ^2) = M_{p}^2\) is the four dimensional squared Planck scale and \(G = R^2 - 4R_{\mu \nu }R^{\mu \nu } + R_{\mu \nu \alpha \beta }R^{\mu \nu \alpha \beta }\) is the GB term. \(\phi \) and \(V(\phi )\) denote the scalar field and the self interaction potential respectively. \(\xi (\phi )\) defines the coupling between scalar field and GB term. Variation of the action with respect to metric and scalar field leads to the field equations as follows,

$$\begin{aligned}&\frac{1}{\kappa ^2}[-R^{\mu \nu }+(1/2)g^{\mu \nu }R] + (1/2)\partial ^{\mu }\phi \partial ^{\nu }\phi \nonumber \\&\quad -\, (1/4)g^{\mu \nu }\partial _{\rho }\phi \partial ^{\rho }\phi + (1/2)g^{\mu \nu }\big [-V(\phi )+\xi (\phi )G\big ] \nonumber \\&\quad -\, 2\xi (\phi )RR^{\mu \nu } - 4\xi (\phi )R^{\mu }_{\rho }R^{\nu \rho } - 2\xi (\phi )R^{\mu \rho \sigma \tau }R^{\nu }_{\rho \sigma \tau } \nonumber \\&\quad +\, 4\xi (\phi )R^{\mu \rho \nu \sigma }R_{\rho \sigma } + 2[\nabla ^{\mu }\nabla ^{\nu }\xi (\phi )]R-2g^{\mu \nu }[\nabla ^2\xi (\phi )]R \nonumber \\&\quad -\, 4[\nabla _{\rho }\nabla ^{\mu }\xi (\phi )]R^{\nu \rho } - 4[\nabla _{\rho }\nabla ^{\nu }\xi (\phi )]R^{\mu \rho } + 4[\nabla ^2\xi (\phi )]R^{\mu \nu }\nonumber \\&\quad +\,4g^{\mu \nu }[\nabla _{\rho }\nabla _{\sigma }\xi (\phi )]R^{\rho \sigma } + 4[\nabla _{\rho }\nabla _{\sigma }\xi (\phi )]R^{\mu \rho \nu \sigma } = 0, \end{aligned}$$
(2)

and

$$\begin{aligned} g^{\mu \nu }[\nabla _{\mu }\nabla _{\nu }\phi ] - V'(\phi ) - \xi '(\phi )G = 0, \end{aligned}$$
(3)

where a prime denotes the derivative with respect to \(\phi \).

The metric for the interior is assumed to be a spatially flat Friedmann-Robertson-Walker metric and is given by

$$\begin{aligned} ds^2 = -dt^2 + a^2(t)\big [dr^2 + r^2d\theta ^2 + r^2\sin ^2{\theta }d\varphi ^2\big ], \end{aligned}$$
(4)

where the scale factor a(t) governs the time evolution of the interior spacetime. For such a metric, the expression of Ricci scalar R and GB term G take the following form,

$$\begin{aligned} R= & {} 6[2H^2 + \dot{H}],\\ G= & {} 24H^2[H^2 + \dot{H}], \end{aligned}$$

where \(H = \dot{a}/a\) and dot denotes the derivative with respect to time (t).

Here, we have assumed the interior of the collapsing star to be consisting of the scalar field as well as a perfect fluid. Therefore, using the metric in eqn. (4), the field equations can be written as,

$$\begin{aligned}&-(3/\kappa ^2)H^2 + (1/2)\dot{\phi }^2 + V(\phi ) + 24H^3\dot{\xi } + \rho = 0, \end{aligned}$$
(5)
$$\begin{aligned}&\quad \frac{1}{\kappa ^2}\big [2\dot{H}+3H^2\big ] + (1/2)\dot{\phi }^2 - V(\phi ) - 8H^2\ddot{\xi } \nonumber \\&\qquad - 16H\dot{H}\dot{\xi } - 16H^3\dot{\xi } + p = 0, \end{aligned}$$
(6)

and

$$\begin{aligned} \ddot{\phi } +3H\dot{\phi } + V'(\phi ) + 24\xi '(\phi )(H^4+H^2\dot{H}) = 0. \end{aligned}$$
(7)

\(\rho \) and p signifies the density and pressure of the constituent fluid inside the collapsing star.

3 Exact solution

It can be easily noted that, due to the inclusion of the additional fluid component (two additional unknown functions \(\rho \) and p), the system of equations becomes even more difficult to deal with analytically, even more so without assuming any particular equation of state. However, we try to do away with the difficulty by incorporating a strategy, where the scalar field evolution Eq. (7) is identified as an anharmonic oscillator equation and integrated straightaway without any apriori assumptions regarding the equation of state or the scale factor. In this way, the other field equations can be used to study the evolution of the fluid/scalar field distribution in a general manner. However, the only assumption being implemented here is that the Eq. (7) is integrable. The criterion for such an integrability can be defined in terms of an invertible point transformation, first worked out by Duarte et. al. [62], Euler et al. [63]. Although the main motivation over this assumption of integrability is extracting information from the field equations at any cost without resorting to any assumption of equation of state, this assumption is by no means unphysical. In recent past, this approach has proved to be very useful to produce interesting solutions depicting dynamics of minimally coupled massive scalar field [61], self-similar solutions [47] and also showed promise towards being a tool for reconstructing modified gravity lagrangians. Using this approach, Chakrabarti, Said and Farrugia have, quite recently studied a reconstruction method for teleparrallel gravity [64]. This work therefore, carries a subtle motivation of testing the scope of this approach in the domain of Scalar-Einstein-Gauss-Bonnet gravity.

While a simple linear harmonic oscillator has a straightforward sinusoidal solution, an anharmonic oscillator has more contributing terms and can represent more physical features of a dynamical system. This takes the form of a nonlinear second order differential equation with variable coefficients as

$$\begin{aligned} \ddot{\phi }+f_1(t)\dot{\phi }+ f_2(t)\phi +f_3(t)\phi ^n=0, \end{aligned}$$
(8)

where \(f_i\) are functions of t and \(n \in \mathcal{Q}\) is a constant. Overhead dot represents differentiation with respect to cosmic time, t. Using Euler and Duartes [62, 63, 65] work on the integrability of the general anharmonic oscillator equation and the more applicable reproduction by Harko et al. [66], this equation can be integrated under certain conditions. The essence of the integrability criterion is that, an equation of the form Eq. (8) can be point transformed into an integrable form iff \(n\notin \left\{ -3,-1,0,1\right\} \), provided the coefficients of Eq. (8) satisfy the differential condition

$$\begin{aligned}&\frac{1}{n+3}\frac{1}{f_{3}(t)}\frac{d^{2}f_{3}}{dt^{2}} - \frac{n+4}{\left( n+3\right) ^{2}}\left[ \frac{1}{f_{3}(t)}\frac{df_{3}}{dt}\right] ^{2} \nonumber \\&\quad + \frac{n-1}{\left( n+3\right) ^{2}}\left[ \frac{1}{f_{3}(t)}\frac{df_{3}}{dt}\right] f_{1}\left( t\right) \end{aligned}$$
(9)
$$\begin{aligned}&\quad + \frac{2}{n+3}\frac{df_{1}}{dt}+\frac{2\left( n+1\right) }{\left( n+3\right) ^{2}}f_{1}^{2}\left( t\right) =f_{2}(t). \end{aligned}$$
(10)

Introducing a pair of new variables \(\Phi \) and T given by

$$\begin{aligned} \Phi \left( T\right)= & {} C\phi \left( t\right) f_{3}^{\frac{1}{n+3}}\left( t\right) e^{\frac{2}{n+3}\int ^{t}f_{1}\left( x \right) dx }, \end{aligned}$$
(11)
$$\begin{aligned} T\left( \phi ,t\right)= & {} C^{\frac{1-n}{2}}\int ^{t}f_{3}^{\frac{2}{n+3}}\left( \xi \right) e^{\left( \frac{1-n}{n+3}\right) \int ^{\xi }f_{1}\left( x \right) dx }d\xi , \end{aligned}$$
(12)

where C is a constant, Eq.(8) can then be written in an integrable form as

$$\begin{aligned} \frac{d^{2}\Phi }{dT^{2}}+\Phi ^{n}\left( T\right) = 0. \end{aligned}$$
(13)

We focus our investigation on a particular case of polynomial coupling, i.e., where \(\xi (\phi ) = \xi _{0} \frac{\phi ^2}{2}\). We also take the self-interaction potential \(V(\phi ) = V_{0} \frac{\phi ^{(n+1)}}{(n+1)}\). Both positive powers and inverse powers of \(\phi \) are very useful in the cosmological setting, in particular, inverse power law models are extremely useful as quintessence fields among other interesting properties. With this assumption, the scalar field evolution equation becomes

$$\begin{aligned} \ddot{\phi } +3H\dot{\phi } + 24 \xi _{0} \phi (H^4+H^2\dot{H}) + V_{0} \phi ^n = 0. \end{aligned}$$
(14)

One can easily identify \(f_{1}(t) = 3H\), \(f_{2}(t) = 24 \xi _{0} (H^4 + H^2 \dot{H})\) and \(f_{3} = V_{0}\).

Provided \(n\notin \left\{ -3,-1,0,1\right\} \), the integrability criterion produces a differential equation for \(H = \frac{\dot{a}}{a}\) given by

$$\begin{aligned} 24 \xi _{0} H^4 - 18 \frac{(n+1)}{(n+3)^2} H^2 + \dot{H} \left( 24 \xi _{0} H^2 - \frac{6}{(n+3)}\right) = 0, \end{aligned}$$
(15)

which we rewrite in the form

$$\begin{aligned} \dot{H} = \frac{18\frac{(n+1)}{(n+3)^2} - 24 \xi _{0} H^2}{24 \xi _{0} - \frac{6}{(n+3) H^2}}. \end{aligned}$$
(16)

During the final stages of the evolution of the collapsing scalar field, it is expected that the proper volume is very small and sharply decreasing in nature. Moreover, the rate of collapse must also be increasing rapidly. With that in mind one can say, \(H = \frac{\dot{a}}{a}>> 0\) and is a sharply increasing function of time. Therefore, \(\frac{1}{H^2}\) can be neglected compared to other term on the denominator on the RHS. With that simplification and the fact that \(\dot{a} < 0\) for a collapsing model, we can write a solution for the scale factor, defining the time evolution of the collapsing scalar field as

$$\begin{aligned} a(t)= & {} \frac{1}{6(1+n)} e^{\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)} \nonumber \\&-\, 2 \xi _{0} a_{0} (3+n)^2 e^{-\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)}. \end{aligned}$$
(17)

Here both \(t_{0}\) and \(a_{0}\) are constants of integration and serves as the initial condition of the model alongwith the choice of self-interaction potential (value of n). In order to have a real time evolution one must enforce the restrictions \(n > -1\) and \(\xi _{0} > 0\).

The time of reaching the zero proper volume can be calculated from Eq. (17) as

$$\begin{aligned} t_{s} = t_{0} - \frac{2(3+n)\sqrt{\xi _{0}}}{\sqrt{3(1+n)}} ln [2\sqrt{3 a_{0} \xi _{0} (1+n)} (3+n)]. \end{aligned}$$
(18)

4 Evolution of the scale factor for different initial conditions

4.1 Evolution of the scale factor for \(a_{0} > 0\)

We present different possible outcomes of the collapse graphically and study the evolution varying different initial conditions.

Fig. 1
figure 1

Evolution of the scale factor with time for \(V(\phi ) = \frac{\phi ^4}{4}\) and for different positive values of \(a_{0}\). (Colour Code : \(Blue \rightarrow a_{0} = 0.6\), \(Yellow \rightarrow a_{0} = 0.5\), \(Green \rightarrow a_{0} = 0.4\) and \(Red \rightarrow a_{0} = 0.3\)). For all values of the initial parameter \(a_{0}\), the collapse reaches a zero proper volume, however the rate of collapse depends on the choice of the parameter

Full size image

In Fig. 1, the scale factor is plotted as a function of time for a fixed value of \(t_{0}\), for different positive values of \(a_{0}\) and \(n = 3\), i.e., \(V(\phi ) = \frac{\phi ^4}{4}\). The evolution shows a rapid collapsing behavior For all values of \(a_{0}\), and an ultimate zero proper volume end-state, however, the rapidity of the collapse and the time of formation of zero proper volume depends on the choice of the parameter.

An important note to make here is that, for the particular case of \(n = 3\), i.e., \(V(\phi ) = \frac{\phi ^4}{4}\), \(t_{0} = 20\) and as long as a positive choice of \(a_{0}\) is made, \(\xi _{0}\) must be taken in between 0 and 1. For \(\xi _{0} < 0\), there is no real evolution and for \(\xi _{0} > 1\), the evolution becomes negative. Therefore a condition of \(0 < \xi _{0} \le 1\) must be enforced upon the strength of the coupling. However, this range can differ for a different set of initial conditions, i.e., for a different set of n and \(t_{0}\). We have presented a particular case with the note that for any such model, the strength of the coupling is very important and therefore the restriction over the allowed domain of \(\xi _{0}\) must be accounted for.

Fig. 2
figure 2

Evolution of the scale factor with time for different self-interaction potential defined by \(V(\phi ) = \frac{\phi ^{(n+1)}}{(n+1)}\) and for fixed \(a_{0}\) and \(\xi _{0}\). \((Colour Code : Blue \rightarrow n = 3\), \(Yellow \rightarrow n = 1.5\), \(Green \rightarrow n = 0.5\) and \(Red \rightarrow n = 0.001)\). The two different plots are for different ranges

Full size image

In Fig. 2, we plot the collapsing behavior for \(\xi _{0} = 0.4\) and for different choices of n, i.e., for different choice of self interaction potential defined by \(n = 3, 1.5, 0.5\) and 0.001. For all the cases, the evolution starts at a finite value of the scale factor (as shown in the first figure) and a zero proper volme is reached at a finite future. The time of formation of singularity changes depending on the choice of n (as shown by the figure below).

Fig. 3
figure 3

Evolution of the scale factor with time for \(V(\phi ) = \frac{\phi ^4}{4}\), fixed \(a_{0}\) and different values of \(\xi _{0}\)

Full size image

In Fig. 3, we show the evolution of the scale factor with time for different choice of the coupling \(\xi _{0}\), fixing other initial parameters. For different value of \(\xi _{0}\) the collapse ends up in a zero proper volume afterall, but for different \(t_{s}\), which is also evident from Eq. (18).

Fig. 4
figure 4

Evolution of the scale factor with time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(\xi _{0} = 0.4\) and for different negative values of \(a_{0}\). \((Colour Code : Blue \rightarrow a_{0} = -10\), \(Yellow \rightarrow a_{0} = -1\), \(Green \rightarrow a_{0} = -0.1\) and \(Red \rightarrow a_{0} = -0.00001)\). The evolution suggests that the scalar field experiences a collapse initially, only until a critical point after which the collapse changes into an expanding phase

Full size image

4.2 Evolution of the scale factor for \(a_{0} < 0\)

However, depending on the initial conditions, the collapse may not always lead to a zero proper volme. As shown in Fig. 4, for all negative values of \(a_{0}\), the system experiences a collapse initially, but only until a critical point (a non-zero minimum radius) after which it can not shrink further and experiences rapid expansion. The plot here is for \(V(\phi ) = \frac{\phi ^4}{4}\) and \(\xi _{0} = 0.4\). This behavior is valid for all negative values of \(a_{0}\).

Fig. 5
figure 5

Evolution of the scale factor with time for different n, i.e., for different choices of \(V(\phi ) = \frac{\phi ^{(n+1)}}{n+1}\), \(\xi _{0} = 0.4\) and for a particular negative value of \(a_{0} = -1\). \((Colour Code : Blue \rightarrow n = 3\), \(Yellow \rightarrow n = 1.5\), \(Green \rightarrow n = 0.5\) and \(Red \rightarrow n = 0.001)\)

Full size image

We present the evolution graphically for different choices of the self-interaction potential (depending on the choice of n) and for a particular choice of negative \(a_{0}\), and fixed \(\xi _{0}\). The evolution shown in Fig. 5 suggests that the scalar field experiences a collapse initially, only until a critical point after which it experiences a bounce. The overall qualitative behavior remains the same for different n, however, there is an indication that, depending on n, the bouncing behaviour may be scaled.

Fig. 6
figure 6

Evolution of the scale factor with time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(a_{0} = -1\) and for a different values of \(\xi _{0}\). \((Colour Code : Blue \rightarrow \xi = 0.4\), \(Yellow \rightarrow \xi = 1\), \(Green \rightarrow \xi = 2\) and \(Red \rightarrow \xi = 3)\)

Full size image

In Fig. 6, we plot the scale factor as a function of time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(a_{0} = -1\) and for a different values of the coupling parameter \(\xi _{0}\). The evolution suggests that the overall qualitative behavior may change depending on the value of \(\xi _{0}\), in the sense that, the phase of an initial contraction of the scale factor depends on the choice of \(\xi _{0}\). For instance, for \(\xi _{0} = 0.4\), (shown by the blue curve) there is an initial collapsing phase before reaching an eventual non zero minimum cutoff, and the expanding phase begins thereafter. However, if one increases the value of \(\xi _{0}\) gradually, it can be seen (Yellow for \(\xi = 1\), Green for \(\xi = 2\) and Red for \(\xi = 3\)) that the initial collapsing phase seems to become negligible and the entire solution turns into an expanding solution.

5 Evolution of the scalar field

Using the defined point transformations (11) and (12), a general solution of the scalar field can be constructed from the transformed integrable form of the scalar field evolution Eq. (13). The solution is given in the form

$$\begin{aligned} \frac{dT}{d\Phi } = \frac{1}{\sqrt{2 \Big (C_{0} - \frac{\Phi ^{(n+1)}}{(n+1)} \Big )}}. \end{aligned}$$
(19)

Calculating the coefficients \(f_{1}(t) = 3 H\), \(f_{2}(t) = 24 \xi _{0} (H^4 + H^2 \dot{H})\) and \(f_{3} = V_{0}\) using the solution for scale factor (17), one appears at a differential equation governing the behaviour of the scalar field itself given by

$$\begin{aligned} \Big (\dot{\phi } + \frac{2\phi }{(n+3)} \Big )^2 = \frac{C^{-(n+1)} a^{-6\frac{(n+1)}{(n+3)}}}{2C_{0} + \frac{2}{(n+1)}\Big [C \phi a^{\frac{6}{(n+3)}}\Big ]^{(n+1)}}. \end{aligned}$$
(20)

It is extremely non-trivial to solve this equation analytically without any choice of n. We present a particular case here. Using \(C = 1\), \(C_{0} = 0\) and investigating the equation numerically for \(n = 3\) (i.e., for \(V(\phi ) = \frac{\phi ^4}{4}\)), we numerically solve the equation and the solution is written as

$$\begin{aligned} \phi (t) = e^{-\frac{1}{3}} \Bigg [C_{1} + 3 \int {\frac{g_{1}(t)}{(g_{2}(t) - \delta )^{4}}} \Bigg ]. \end{aligned}$$
(21)

Here, \(g_{1}(t)\) and \(g_{2}(t)\) are functions of time which we have written in this form for the sake of brevity.

ALthough it would have been better if a neat and closed for of the evolution could be written, however, the numerical integration of the Eq. (21) produces a very large expression for the scalar field (about 145 terms). We present in brief the functional form so as to give an idea regarding how the scalar field can, in principle, evolve.

$$\begin{aligned}&\phi (t) = \psi (t)^{\frac{1}{3}}, \\&\psi (t) = s_{0} e^{-t} + \frac{g(t)}{j(t)}, \\&j(t) = s_{1} \left( 1728 b e^{\frac{\sqrt{\frac{1}{m}}t}{\sqrt{3}}} m-e^{\frac{d \sqrt{\frac{1}{m}}}{\sqrt{3}}}\right) ^{3}, \\&g(t) \sim \Bigg [a_{1}e^{a_{2}-t}\Bigg (a_{i} + a_{j}e^{a_{k}+a_{l}t}\Bigg )\Bigg ] + b_{i}e^{b_{j} + b_{k}t}\\&{_2}F_1\left( 1,2+\sqrt{3\xi _{0}};3+\sqrt{3\xi _{0}};1728 \xi _{0} a_{0} e^{\frac{1-t_{0}}{\sqrt{3\xi _{0}}}}\right) \\&\quad + c_{i}e^{c_{j} + c_{k}t} {_2}F_1\left( 1,2+\sqrt{3\xi _{0}};3+\sqrt{3\xi _{0}};1728 \xi _{0} a_{0} e^{\frac{t-t_{0}}{\sqrt{3\xi _{0}}}}\right) . \end{aligned}$$

All the parameteres \(a_{1}\), \(a_{2}\), \(a_{i}-\)s, \(a_{j}-\)s etc, are infact defined in terms of the parameters of the theory, i.e., \(\xi _{0}\), n and \(a_{0}\). The function g(t) consists of a total of \((80+33+31)\) terms where there are 80 terms of the form of \(\Bigg (a_{i} + a_{j}e^{a_{k}+a_{l}t}\Bigg )\), 33 terms of the form of \(b_{i}e^{b_{j} + b_{k}t}\) and 31 terms of the form of \(c_{i}e^{c_{j} + c_{k}t}\). It is quite obvious that a numerical examination is absolutely necessary over the time evolution of the scalar field. We note that, the parameters \(\xi _{0}\) (i.e., the strength of the coupling of scalar field with the Gauss-Bonnet term) and the constant of integration \(C_{1}\) play an important part in determining the behavior of the scalar field afterall. In the next subsection we present the numericla results of the evolution of the scalar field as a function of time for different parameters.

Fig. 7
figure 7

Scalar field as a function of time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(\xi _{0} = 0.4\) and for different values of \(C_{1}\) (\(Color Code : Blue \Rightarrow C_{1} = 0.001, Yellow \Rightarrow C_{1} = 0.01, Green \Rightarrow C_{1} = 0.1\) and \(Red \Rightarrow 1.0\))

Full size image

5.1 Evolution of the scalar field for \(a_{0} > 0\), i.e., for collapsing scale factor

In Fig. 7, the scalar field is plotted as a function of t for \(V(\phi ) = \frac{\phi ^4}{4}\), \(\xi _{0} = 0.4\) and for different values of \(C_{1}\) (\(C_{1} = 0.001, C_{1} = 0.01, C_{1} = 0.1 \mathrm{and} 1.0\)). The scalar field diverges around the time of formation of singularity. However, depending on the value of \(C_{1}\), the nature of the evolution before reaching singularity may be a little different but the qualitative behavior is this; the time evolution of the scalar field starts at a finite value and then decreases more-or-less steadily, before reaching a minimum critical value. Thereafter, the scalar field increases with time rapidly and diverges.

Fig. 8
figure 8

Evolution of the scalar field as a function of time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(C_{1} = 0.01\) and for different choices of coupling parameter \(\xi _{0} = 0.4 (Blue)\) and \(\xi _{0} = 0.5 (Yellow)\)

Full size image

Evolution of the scalar field as a function of time is plotted in Fig. 8 for \(V(\phi ) = \frac{\phi ^4}{4}\), \(C_{1} = 0.01\) and for different choices of \(\xi _{0} = 0.4\) and 0.5. For different values of the coupling parameter \(\xi _{0}\), the scalar field diverges at different time. This is quite consistent because the time of formation of singularty depends on \(\xi _{0}\). The same can be verified from Eq. (18) as well.

5.2 Evolution of the scalar field for \(a_{0} < 0\)

As discussed by Figs. 3, 4 and 5, the evolution is not always collapsing forever, depending on the choice of the parameter \(a_{0}\). For \(a_{0} < 0\), the scale factor experiences a transition from a state of contraction into a rapid expansion. Here we show the behavior of the scalar field with time for such cases, i.e., for \(a_{0} < 0\).

Fig. 9
figure 9

Evolution of the scalar field as a function of time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(C_{1} = 10\), \(\xi _{0} = 0.4\) and \(a_{0} = -1\)

Full size image

In Fig. 9, we have plotted \(\phi (t)\) vs t for \(V(\phi ) = \frac{\phi ^4}{4}\), \(C_{1} = 10\), \(\xi _{0} = 0.4\) and \(a_{0} = -1\). The scalar field starts at some finite value and gradually decreases, exhibiting some periodic behavior with time. Eventually it decays into a negligibly small positive value as shown in the figure.

Fig. 10
figure 10

Evolution of the scalar field as a function of time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(C_{1} = 0.01\), \(a_{0} = -1\) and for different choices of \(\xi _{0} = 0.4 (Blue)\), \(\xi _{0} = 1 (Yellow)\) and \(\xi _{0} = 2 (Green)\)

Full size image

However, the periodic time evolution of the scalar field is sensitive over the choice of \(\xi _{0}\) as shown in figure 10. The periodic nature seems to be absent as one gradually increases the value of \(\xi _{0}\), (here, plotted for \(\xi _{0} = 0.4, 1\) and 2).

6 Evolution of the strong energy condition

A collapsing perfect fluid is physically reasonable if it obeys the strong energy condition which is satisfied if for any timelike unit vector \(w^{\alpha }\) and the following inequality holds

$$\begin{aligned} 2 T_{\alpha \beta } w^{\alpha } w^{\beta } + T \ge 0, \end{aligned}$$
(22)

where T is the trace of the energy momentum tensor. The energy conditions were investigated in details for imperfect fluids by Kolassis, Santos and Tsoubelis [67]. Following their work, we investigate the validity of the strong energy condition (\((\rho +3p) > 0\)) for our model. The strong energy condition can be violated only if the total energy density is negative or if there exists a large negative principal pressure of \(T^{\alpha \beta }\).

Since we have studied the solution for the scale factor from the scalar field evolution (7) equation straightaway, the field Eqs. (5) and (6) can be used to study the evolution of the constituent fluid density and pressure. We numerically study the strong energy condition and present it’s nature by the following plots.

Fig. 11
figure 11

Evolution of \((\rho +3p)\) with time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(t_{0} = 20\), \(a_{0} = 10\), \(\xi _{0} = 0.4\) and for different choices of \(C_{1}\) (\(C_{1} = 0.1 (Blue)\), \(C_{1} = 10 (Yellow)\), \(C_{1} = 100 (Green)\) and \(C_{1} = 10000 (Red)\))

Full size image

In Fig. 11, we plot the evolution of \((\rho +3p)\) as a function of time for a particular choice of potential \(V(\phi ) = \frac{\phi ^4}{4}\), positive \(a_{0}\) (ensuring the collapsing nature of the solution), \(\xi _{0} = 0.4\) and for different choices of the initial conditon \(C_{1}\) (\(C_{1} = 0.1 (Blue), 10 (Yellow), 100 (Green)\) and 10000(Red)). The evolution suggests that that \((\rho +3p) \ge 0\) is satisfied throughout the evolution of the collapsing body, before reaching the curvature singularity where \((\rho +3p)\) increases sharply, as both pressure and density diverges eventually.

Fig. 12
figure 12

Evolution of \((\rho +3p)\) with time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(t_{0} = 20\), \(a_{0} = 10\), \(\xi _{0} = 0.4\) and for very small values of \(C_{1}\) (\(C_{1} = 0.001 (Blue)\) and \(C_{1} = 0.0001 (Yellow)\)

Full size image

However, depending on the value of \(C_{1}\), we also find some cases where the strong energy condition may be violated during the collapsing evolution, before eventually diverging at the singular epoch. In Fig. 12, we plot the evolution of \((\rho +3p)\) as a function time for \(V(\phi ) = \frac{\phi ^4}{4}\), positive \(a_{0}\), \(\xi _{0} = 0.4\) and for very small positive values of \(C_{1}\) (\(C_{1} = 0.001 (Blue)\) and \(C_{1} = 0.0001 (Yellow)\).

However, the evolution of the strong-energy condition in time leading to the conclusion that this is often violated, is not really an unexpected outcome in theories of gravity that contain strong-curvature terms. In the present case one can argue that the Gauss-Bonnett term in priciple can create an effective energy-momentum tensor whose contribution to the total energy-momentum tensor can lead to the violation of the strong-energy condition. Therefore, it may not be the nature of the perfect fluid afterall that violates the strong-energy condition. We also note here that the energy conditions of general relativity are a mathematical way of making the notion of locally positive energy density by stating that various linear combinations of the components of the energy momentum tensor must stay non-negative and it is sometimes argued that subtle quantum effects can violate all of the energy conditions. Moreover, there are examples of classical systems that violate all the energy conditions as well (For instance, Lorentzian-signature traversable wormholes [70]). The simplest possible source of classical energy condition violations is from the contribution of scalar fields, in particular, non-minimally coupled scalar field contributions, as worked out in details by Visser and Barcelo [68, 70], Flanagan and Wald [69].

Fig. 13
figure 13

Evolution of \((\rho +3p)\) with time for \(V(\phi ) = \frac{\phi ^4}{4}\), \(t_{0} = 20\), \(a_{0} = 10\) and for different value of \(\xi _{0} = 0.04 (Blue)\) and \(\xi _{0} = 0.4 (Yellow)\)

Full size image

In Fig. 13, the evolution of \((\rho +3p)\) with time is plotted for \(V(\phi ) = \frac{\phi ^4}{4}\), positive \(a_{0}\), \(C_{1} = 10000\) and for different value of \(\xi _{0} = 0.04 (Blue)\) and \(\xi _{0} = 0.4 (Yellow)\). As can be seen from the graph, \((\rho +3p)\) maintains a positive signature throughout the collapse, however, for different value of \(\xi _{0}\) the energy condition goes to positive infinity at different time. This is quite consistent, since we have already discussed that the strength of coupling with the Gauss-Bonnet term \(\xi _{0}\) plays a crucial role in determining the time of formation of singularity.

7 Nature and visibility of the singularity

In order to investigate whether the singularity is a curvature singularity or just an artifact of coordinate choice, one must look into the behavior of the Kretschmann curvature (K) scalar at \(t\rightarrow t_{s} = t_{0} - \frac{2(3+n)\sqrt{\xi _{0}}}{\sqrt{3(1+n)}} ln [2\sqrt{3 a_{0} \xi _{0} (1+n)} (3+n)]\), i.e., when the scale factor a(t) goes to zero. For the metric presented in Eq. (4), K has the expression,

$$\begin{aligned} K = 6\bigg [\frac{\ddot{a}(t)^2}{a(t)^2} + \frac{\dot{a}(t)^4}{a(t)^4}\bigg ] \end{aligned}$$
(23)

Using the solution of a(t) Eq. (17), it is straightforward to deduce that the Kretschmann scalar diverges at the zero proper volume and thus the collapsing body discussed here ends up in a curvature singularity.

7.1 Formation of apparent horizon

Whether the curvature singularity is visible to an exterior observer or not, depends on the formation of an apparent horizon. The condition for such a surface is given by

$$\begin{aligned} g^{\mu \nu }R,_{\mu }R,_{\nu }=0, \end{aligned}$$
(24)

where R is the proper radius of the two-sphere, given by ra(t) in the present case.

Using the exact solution for collapse given by Eqs. (17) and (24), we deduce the condition for formation of apparent horizon as

$$\begin{aligned} \frac{\sqrt{3}x}{12(3+n)\sqrt{(1+n)\xi _{0}}} + \frac{\sqrt{3 \xi _{0} (1+n)} a_{0} (3+n)}{x} - \lambda = 0, \end{aligned}$$
(25)

where \(\lambda \) is a constant of separation and \(x = e^{\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)}\). One can solve the above equation to exactly deduce the time of formtion of an apparent horizon as

$$\begin{aligned} t_{ap}= & {} t_{0} \nonumber \\&- \frac{2(3+n) \sqrt{\xi _{0}}}{\sqrt{3(1+n)}} ln \Big [2(3+n) \sqrt{3(1+n)\xi _{0}} (\lambda \pm \sqrt{{\lambda }^2 - a_{0}})\Big ].\nonumber \\ \end{aligned}$$
(26)

Comparing Eq. (26) with the time of formation of the singularity (18) we arrive at the very important expression

$$\begin{aligned} (t_{s} - t_{ap}) = \frac{2 (3+n) \sqrt{\xi _{0}}}{\sqrt{3(1+n)}} ln \Big [\frac{(\lambda \pm \sqrt{{\lambda }^2 - a_{0}})}{\sqrt{a_{0}}} \Big ]. \end{aligned}$$
(27)

To comment on the ultimate visibility of the collapse outcome, it is necessary to see if nonspacelike trajectories emanate from the singularity and reach a faraway observer. The singularity is at least locally naked if there are future directed nonspacelike curves that reach faraway observers. This is possible if the formation of apparent horizon is delayed or if there is no formation of horizon at all. In the present case the time of formation of singularity is independent of r [given by Eq. (18)] and therefore it is only natural that the entire collapsing system (scalar field and perfect fluid) would reach the singularity simulteneously at \(t = t_{s}\). This kind of singularity is always expected to be covered by the formation of an apparent horizon as discussed by Joshi, Goswami and Dadhich [71]. Here, from Eq. (27), it can be said that the formation of an ultimate covered singularity is dependent over initial conditions such as \(\xi _{0}\), \(\lambda \) and \(a_{0}\) so that \((t_{s} - t_{ap}) > 0\). However, in principle, one can also ask about the possible end state if the initial conditions conspire to make the situation to end up otherwise, i.e., \((t_{s} - t_{ap}) < 0\). In such a case, there is no formation of apparent horizon at all, since all the collapsing shells labelled by different values of r, shrinks to zero proper volume and the physical quantities diverge at \(t = t_{s}\), which is reached before \(t_{ap}\). Therefore the singularity remains naked and the condition for such an end-state can be written from Eq. (27) as

$$\begin{aligned} a_{0} = \frac{4 \lambda ^2 \delta ^2}{(1 + \delta )^2}, \end{aligned}$$
(28)

where \(0< \delta < 1\).

8 Matching with an exterior Vaidya Spacetime

For the sake of completeness, proper junction conditions are to be examined carefully which allow a smooth matching of an exterior geometry with the collapsing interior. First of all it was extensively shown by Goncalves and Moss [42] that any sufficiently massive collapsing scalar field can be formally treated as collapsing inhomogeneous dust in general relativity. Moreover, astrophysical objects undergoing a gravitational collapse can be expected to be in an almost vacuum spacetime, and therefore the exterior spacetime around a spherically symmetric dying star is well described by the Schwarzschild geometry. From the continuity of the first and second differential forms, the matching of the sphere to a Schwarzschild spacetime on the boundary surface, \(\Sigma \), is extensively worked out in literature [72,73,74,75].

However, conceptually this leads to an inconsistency, since Schwarzschild has zero scalar field. Therefore, such a matching would lead to a discontinuity in the scalar field, and a delta function in the gradient of the scalar field. As a consequence, there will appear square of a delta function in the stress-energy, which is definitely an inconsistency. Since we have a scalar field distribution inside the collapsing sphere, it is more physical to match the interior with a Vaidya exterior solution (for a detailed analysis of vaidya matching in a scalar field collapse we refer to [76, 77]) across a boundary hypersurface defined by \(\Sigma \). The metric just inside \(\Sigma \) is,

$$\begin{aligned} d{s_-}^2=dt^2-a(t)^2dr^2-r^2 a(t)^2d{\Omega }^2, \end{aligned}$$
(29)

and the metric in the exterior of \(\Sigma \) is given by

$$\begin{aligned} d{s_+}^2=\left( 1-\frac{2M(r_v,v)}{r_v}\right) dv^2+2dvdr_v-{r_v}^2d{\Omega }^2. \end{aligned}$$
(30)

Matching the first fundamental form on the hypersurface we get

$$\begin{aligned} \Big (\frac{dv}{dt} \Big )_{\Sigma }=\frac{1}{\sqrt{1-\frac{2M(r_v,v)}{r_v}+\frac{2dr_v}{dv}}}, \end{aligned}$$
(31)

and

$$\begin{aligned}&(r_v)_{\Sigma } = r a(t) \\&\quad = r \Bigg [\frac{1}{6(1+n)} e^{\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)} \\&\qquad - 2a_{0}(3+n)^2 e^{-\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)}\Bigg ]. \end{aligned}$$

Matching the second fundamental form yields,

$$\begin{aligned} \Big (r a(t) \Big )_{\Sigma } = r_v\left( \frac{1-\frac{2M(r_v,v)}{r_v}+\frac{dr_v}{dv}}{\sqrt{1-\frac{2M(r_v,v)}{r_v}+\frac{2dr_v}{dv}}}\right) . \end{aligned}$$
(32)

Using Eqs. (31) and (32) one can write

$$\begin{aligned} \left( \frac{dv}{dt} \right) _{\Sigma } = \frac{3r a(t)^2 - r^2}{3(ra(t)^2 - 2Ma(t))}, \end{aligned}$$
(33)

where \(a(t) = \frac{1}{6(1+n)} e^{\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)} - 2a_{0}(3+n)^2 e^{-\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)}\).

From Eq. (32) one obtains

$$\begin{aligned} M_{\Sigma } = \frac{1}{4} \Bigg [ra(t) + \frac{r^3}{9 a(t)^3} + \sqrt{\frac{1}{r a(t)} + \frac{r^3}{81 a(t)^9} - \frac{2r}{9 a(t)^5}} \Bigg ]. \end{aligned}$$
(34)

Matching the second fundamental form we can also write the derivative of \(M(v,r_v)\) as

$$\begin{aligned} M{(r_v,v)}_{,r_v}=\frac{M}{r a(t)} - \frac{2r^2}{9 a(t)^{4}}. \end{aligned}$$
(35)

Equations (32), (33), (34) and (35) completely specify the matching conditions at the boundary of the collapsing scalar field.

9 Conclusion

The present work has been entirely dedicated towards a deeper understanding of a self-interacting scalar field collapse in Scalar-Einstein-Gauss-Bonnet gravity. The scalar field couples non-minimally with the Gauss-Bonnet term by a term quadratic in the scalar field (\(\xi _{0} \frac{\phi ^2}{2}\)). The possibility of the collapse reaching a zero proper volume singularity is seen to be dependent on the coupling parameter \(\xi _{0}\), choice of self-interaction potential and most importantly the initial condition \(a_{0}\) which can be related to the initial radius/volume of the collapsing star. We also comment on the allowed domain of the choice of the coupling parameter \(\xi _{0}\) to have a real evolution of the collapse. It is observed that the collapse ends up in a curvature singularity, where the Kretschmann curvature scalar blows up, alongwith density and pressure of the constituent fluid and the scalar field. The strong energy condition is showed to be valid throughout the collapse. However, depending on certain initial conditions defining the initial distribution of the scalar field, the energy conditions may be violated, which can be attributed to the introduction of a non-minimal coupling of a scalar field in the lagrangian, which is known to be a possible classical system that can violate almost all of the energy conditions [68, 70].

For the sake of completeness we match the interior collapsing solution with the an exterior Vaidya geometry on a boundary hypersurface, since the presence of a Gauss-Bonnet term non-minimally coupled to a scalar field generates a non-zero effective energy momentum tensor arising from spacetime curvature. Therefore matching with an exterior vacuum solution in the presence of Gauss-Bonnet term may lead to inconsistency. Very recently this has been investigated by Banerjee and Paul [54].

The time of formation of the singularity is independent of r, which suggests that all the collapsing shells labelled by different values of r collapses simulteneously when the zero proper volume is reached. Such a singularity is always expected to be covered by a horizon as far as similar studies under the domain of GR are concerned. We note here that in the present case of Scalar-Einstein-Gauss-Bonnet gravity with a polynomial coupling, the ultimate end-state of a covered singularity is conditionally consistent with the corresponding results in GR. For certain initial conditions [defined by Eqs. (27) and (28)] there is a probability of an end state where there is no possibility of formation of horizon at all.

We also observe some interesting results, for instance, depending on the signature of the aforementioned parameter \(a_{0}\), the evolution of scale factor suggests that all of the possible collapsing scenario may not lead to a curvature singularity. Rather, for \(a_{0} < 0\), the fluid undergoes contraction only unto a minimum non-zero cut-off radius, after which it goes into a rapidly expanding phase. It is also seen that the scalar field itself in such cases, decreases monotonically with time, exhibits certain periodic behavior, before becoming negligibly small. This can be somewhat compared with the phenomena of collapse and dispersal for a scalar field in general relativity which have drawn considerable interest in recent years, mainly from a numerical perspective, subtly pointing towards the existence of a critical phenomena is the whole collapsing picture (For details, we refer to the monograph by Gundlach [45, 46]). Recently Bhattacharya, Goswami and Joshi worked out a sufficient condition for the dispersal to take place for a collapsing scalar field in GR that initially begins with a contraction and showed that the transition of the collapsing body into expanding nature is crucially connected wiith the change of the gradient of the scalar field [78]. In the present case we can note that the signature of the parameter \(a_{0} > 0\) or \(a_{0} < 0\) is the defining factor of the transition; i.e., whether or not the collapse will lead into a zero proper volume or a dispersal of scalar field shall take place after the scale factor reaches a minimum cut-off volume. Since \(a_{0}\) comes from the expression defining the scale factor (\(a(t) = \frac{1}{6(1+n)} e^{\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)} - 2a_{0}(3+n)^2 e^{-\frac{\sqrt{3(1+n)}}{2(3+n)\sqrt{\xi _{0}}}(t_{0}-t)}\)), the initial volume of the collapsing system can be predicted as a defining factor connected to the value of \(a_{0}\).

We conclude with the note that this is indeed a simple case, in the sense that we have considered spatial homogeneity in the metric components as well as the scalar field. However, the solution found is simple enough to encourage further allied investigations in this direction, such as, the possibility of a collapse even when the energy conditions are violated may subtly direct one’s attention towards a possible clustering of dark energy distribution. Apart from that, this work also helps further expansion of the usefulness and scope of the particular method of integrability of anharmonic oscillator used here. The theorem which inspires this method is self sufficient as was discussed by Euler [65]. The same has been proved in the context of a massive scalar field minimally coupled to gravity by Chakrabari and Banerjee [61], where the solutions found by virtue of this theorem indeed solve the Klein Gordon type evolution equation once they are put back in the equation. In the current context, the solutions are far more complicated to say the least, as is evident from the expression of the scalar field, as worked out in details in section V. Since the criterion for integrability was investigated under the expectation that the proper volume would be very small and sharply decreasing in nature so that one can neglect the \(\frac{1}{H^2}\) term, to justify such an arguement one can put the solutions (17) and (21) together into the equation (14) and study for different initial conditions defining the scalar field. It can be confirmed that for relevant cases, the solutions put into the Klein-Gordon type Eq. (14) yields a number very close to zero (\(\sim 10^{-8}\) or lower). This approach of point transforming the klein-gordon type equation for the scalar field, to extract the solution out of a seemingly impossible non-linear system has worked really well in the past while investigating different setups of massive scalar field collapse and also inspired a handy reconstruction technique which helps one to assess the particular form of the lagrangian of a modified theory of gravity [64]. We note here that, although the assumption of integrability of the scalar field evolution equation is inspired only from a mathematical perspective, solutions found by means of this assumption are by no means unphysical. Used properly, this method can potentially be useful for allied investigations as well, for instance, a detailed definition of the possible bounds over the choice of coupling function \(\xi \) (somewhat similar to a possible Higgs-Kreschmann invariant coupling in white dwarfs and neutron stars, as shown by Wegner and Onofrio [79, 80]). Further investigation under the setup of a Scalar-Einstein-Gauss-Bonnet gravity in a more generalized scenario (inclusion of factors spatial inhomogeneity, pressure anisotropy, heat flux etc) can be done using this method properly and more rigorously and will be reported elsewhere in future.