Collapse of an axion scalar field

The manuscript deals with an interacting scalar field that mimics the evolution of the so-called axion scalar dark matter or axion like particles with ultra-light masses. It is discussed that such a scalar along with an ordinary fluid description can collapse under strong gravity. The end state of the collapse depends on how the axion interacts with geometry and ordinary matter. For a self-interacting axion and an axion interacting with geometry the collapse may lead to a zero proper volume singularity or a bounce and total dispersal of the axion. However, for an axion interacting with the ordinary fluid description, there is no formation of singularity and the axion field exhibits periodic behavior before radiating away to zero value. Usually this collapse and dispersal is accompanied by a violation of the null energy condition for the ordinary fluid description.


Introduction
General theory of relativity (GR) describes gravity as an artifact of geometry and is full of counter-intuitive ideas and solutions. In particular, we focus on an idea that a massive stellar distribution written under the scope of GR can in principle collapse indefinitely to form a geodesic incompleteness, where curvature scalars diverge. This region is famously known as a Spacetime Singularity. An amendment is sometimes proposed in the theory in the form of a 'Cosmic Censorship Conjecture (CCC)' [1,2]. The conjecture predicts the singularity to remain covered through a development of trapped surface, in process producing a Black Hole end-state [3][4][5]. However, without any concrete proof of the conjecture, nothing stands on paper to favor the Black Holes over it's sibling outcome, the Naked Singularities, who can share information with a faraway observer [6][7][8][9]. Indeed, GR allows existence of such solutions where a collapsing star fails to form any trapped surface/horizon and keeps the ultimate sina e-mail: soumya.chakrabarti@saha.ac.in (corresponding author) gularity visible [10][11][12]. This is no doubt problematic for the predictability of events in the spacetime. The time evolution of a collapsing stellar body is the subject that inspires these inconclusive questions and the topic receives significant attention even today, in GR and it's viable modifications as well [13][14][15].
Modified gravity is a prospective avenue to play with the idea of CCC as here one can try to produce non-trivial geometric corrections, affecting the nature and understanding of the conjecture. Most of the modified gravity models are also motivated from cosmological requirements. They can describe the observed cosmic acceleration of our universe. A particularly popular candidate in this regard is a time evolving interacting scalar field that can produce an effective repulsive effect and fit in nicely as the so-called dark energy (DE) component of the universe. Any such model must also explain the large scale structures and the proper profile of relic radiation from big bang for the universe. This requires an additional feature which can drive the cosmic deceleration and describe the origin of the so-called dark matter (DM) component. The nature of this DM is one of the most studied cosmological question of modern era [16]. A pressureless 'cold dark matter' or the CDM can explain the large scale structure formation of the universe [17], however, numerical and observational evidences from galaxies and their clusters have imposed strong constraints on this description. Many alternative models work well to fit in as DM, for example the Weakly interacting massive particle models [18,19], the warm DM [20,21], self interacting dark matter models [22,23]). One particular approach has received recurring attention in the past decade, a Scalar Field Dark Matter (SFDM). This involves a time dependent scalar field endowed with an interaction potential that fits in the role of DM. It was first put forward Ji and Sin [24] and Sin [25], and received considerable treatment thereafter by Sahni and Wang [26], Hu et al. [27], Matos et al. [28], Matos and Urena-Lopez [29], Arbey et al. [30][31][32] and Matos and Arturo Ureña-López [33].
We are interested in a particular member of the SFDM family, the QCD inspired 'Axion', described by a field [34] or a system of axion-like particles (ALP) with ultra-light masses [35]. An axion can be described by a scalar field with a very low mass and a self-interaction potential. Astrophysical and cosmological observations and compact binary system observations have put some interesting constraints on the axionlike particles and their masses, which have been well studied in literature [36][37][38][39]. We discuss that these fields or ALP-s can collapse due to strong gravity and form infinitely dense points in spacetime or singularities, gravitationally bound objects or dynamical equilibrium states [40][41][42][43] through a scalar field collapse. On it's own, collapse of scalar fields is an interesting subject which tries to answer if the CCC is satisfied by such fundamental matter fields [44][45][46][47][48]. Moreover, it can be proved that an interacting scalar field can reproduce the evolution of different type of matter distributions [49]. Instability of such collapsing massive scalar fields and the probability of formation of a primordial black hole has also been an interesting topic of discussion over the years [50]. The subject gained additional popularity after the discovery of the 'Critical Phenomena' in the collapse of zero mass scalar fields [51][52][53][54]. The phenomenon indicates that a scalar field can collapse and form either a singular end-state or disperse away to zero depending on a finite number of critical parameters.
We study the collapse of the dark matter candidate axion or ALP-s. One can also consider these objects as collapsing axion stars. A dynamical equilibrium of such a collapsing system can originate from the interplay of kinetic pressure, self-gravity or the axion interactions [55]. We write the axion description as an interacting scalar field. Such an axion distribution has received some limited interest only recently where possible end states of the collapse is identified [56] numerically. We work with a spatially homogeneous Oppenheimer-Snyder like geometry and discuss three examples with step by step generalization of the scalar interaction. (a) A minimally coupled axion field with self-interaction potential. (b) An axion field coupling non-minimally with geometry alongwith a self-interaction potential. (c) An axion field interacting nonminimally with both geometry and ordinary matter. These considerations are motivated from a thought that during the final stages of a gravitational collapse, where the strength of gravity and the spacetime curvature is expected to be quite high, non-trivial interactions between scalar field and geometry and even normal matter can exist. A similar scalar configuration is popular from a cosmological purview, known as Chameleon Fields, who can provide a very smooth transition from a decelerated to an accelerated phase of expansion of the universe [57][58][59][60][61]. Apart from the axion field and it's interactions, we also consider a perfect fluid distribution without any restriction over the equation of state parameter at the outset. This is motivated from a requirement to describe the distribution of DE by the aforementioned fluid distribution. No concrete information on the distribution of the DE component is known apart from the speculation that it does not cluster below Hubble scale. The main focus of the present work is the time evolution of collapsing axion Dark Matter and formation of gravitationally bound end-states. However, the collapsing perfect fluid alongwith the scalar field motivates one to think about the possibility of clustering of DE.

Theoretical setup and methodology
We assume that the dark matter distribution is overall given by a scalar field with a self-interaction. We write the generic axion potential V (φ) as [56] m is the mass and f is the decay constant of the ALP-s. Usually the ALP dark matter model predicts the density parameter to be given by In the present work we take m and f as free parameters and assume f >> m. Under this assumption we write the axion potential in Eq. (1) as This form has an interesting resemblance with the Higgs potential, which has played significant role in the cosmological context as in inflation [62,63] and cosmological reconstruction of modified gravity [64][65][66]. In the rest of the manuscript, we study scalar field collapse where the selfinteraction of the scalar field is given by Eq. (3). Three different setups are considered, (a) a minimally coupled self-interacting axion, (b) a non-minimally coupled selfinteracting axion and (c) a third case where the axion interacts with the ordinary matter. For all three cases, we choose a spatially flat homogeneous metric of Oppenheimer-Snyder type, given by We restrict our study to collapsing models only. This means that radius of the two-sphere is a monotonically decreasing function of time, i.e.,ȧ < 0. We follow an unconventional methodology compared to the standard methods of studying exact solutions of scalar field cosmology or collapse. We deal with the Klein Gordon equation of the axion scalar field by incorporating a purely mathematical property of general second order differential equations with variable coefficients. The property involves point transforming the equations into an integrable form and is derived from the symmetry analysis of a general classical anharmonic oscillator equation system [67][68][69][70]. The general equation is written as f 1 , f 2 and f 3 are unknown functions of some variable. n is a constant. A transformation of this equation into an integrable form requires a pair of point transformations and the condition n / ∈ {−3, −1, 0, 1} to be satisfied. Moreover, the coefficients must satisfy the condition The point transformations are written as where C is a constant. Using this property, we solve the axion scalar evolution equation assuming it's integrability at the outset and use the other field equations to assess the constraints on the fluid energy density components. The main motivation of assuming this integrability comes out of a pure mathematical curiosity. However, by no means this produces unphysical solutions. The scope of this approach has been discussed at length quite recently, in the context of simple scalar field collapse [71,72], scalar-gauss-bonnet gravity [73] and cosmological reconstruction of modified theories of gravity [64][65][66].

Minimally coupled axion
The action for a minimally coupled axion can be written as where L m is the Lagrangian density for ordinary fluid distribution. Energy-momentum contribution from the scalar field φ therefore is We assume that the scalar field is spatially homogeneous, i.e., φ = φ(t). Therefore we write the field equations as We also write the scalar evolution equation as Using the potential in Eq. (3) the scalar evolution equation is written as From this point on, we call the Klein-Gordon equation for the axion scalar field evolution as the axion evolution equation. The overhead dot denotes the derivative with respect to t. Now Eq. (14) is clearly a special case of Eq. (5) with n = 3 and f 1 = 3ȧ a , f 2 = m 2 and f 3 = − m 2 6 f 2 . Our progress involves the assumption that we can transform Eq. (13) into an integrable form, writing a transformation like Eq. (7). Therefore the criterion in Eq. (6) produces an equation governing the allowed dynamics of the radius of two-sphere a(t) and the Eq. (7) produces a solution for the scalar field. The equation of the radius of two-sphere is written beloẅ A simple first integral of this equation can be written aṡ where a 0 is a non-zero constant of integration. Forȧ < 0, the above equation is solved to write We plot the evolution of a(t) in Figs. 1 and 2 for different initial conditions. Figure 1 is the plot of Eq. (17) for a 0 > 0. It is clear that the stellar distribution experiences a uniform The graph on top shows the evolution for a fixed value of m, but different choices of a 0 . On the other hand, the graph on bottom shows the evolution for a fixed value of a 0 > 0, but different values of m. Figure 2 is the plot of Eq. (17) for a 0 < 0. In this case, the stellar distribution experiences an initial collapse. Eventually the collapse decelerates and goes through a signature flip ofȧ during it's evolution. This indicates a minimum cutoff for the radius after which the collapsing axion distribution bounces without reaching a zero proper volume.
We study the evolution of the scalar field as a function of time using Eq. (7) and transforming the evolution equation Eq. (13). A similar solution can be found by re-inserting the solution Eq. (17) in Eq. (13). The exact solution is complicated to write in a closed form. Therefore we resort to numerical study and plot the evolutions in Figs  is for initial condition da(t) dt t=t i > 0 and the graph on bottom are for initial condition da(t) Fig. 3 it is clear that for all a 0 > 0, the axion field diverges at a finite time, just about the time of formation of a zero proper volume. This indicates a formation of a spacetime singularity. On the other hand, if one chooses a negative value of the critical parameter a 0 , the scalar field initially experiences a sharp increasing time evolution. Interestingly, about the time the scale factor changes nature from collapsing to bouncing, the scalar field changes nature as well. It falls off sharply as a function of time, eventually dispersing away to zero value at an asymptotic future (see Fig. 4).
The curvature scalars can be written as R = −6 ä a +ȧ 2 a 2 and K = 6 ä 2 a 2 +ȧ 4 a 4 . Using Eq. (17) it is straightforward to check that for a 0 > 0 both Ricci and Kretschmann scalars diverge to infinity when a(t) → 0. Therefore, if the axion distribution somehow reaches zero proper volume it definitely ends up in a curvature singularity. Otherwise, it bounces indefinitely after reaching the cutoff radius and ends up radiating away to zero. The distribution of the additional fluid existing alongside the axion scalar can be studied, courtesy of the field Eqs. (11) and (12). We write and Therefore the null energy condition for the fluid description (NEC) is It is straightforward that the positivity of the NEC depends on the nature of the first term on the RHS, −2 d dt ȧ a , since −φ 2 is always negative. We assess the nature of − d dt ȧ a from The Raychaudhuri Equation which essentially is a geometric relation to write the dynamics of flows in terms of mean separation between a congruence of curves [74][75][76][77]. Quite recently the role of this purely kinematic equation of flows in gravitational collapse of fields or fluid distribution has drawn attention [78]. For a system of timelike geodesics whose tangent vector field is given by u μ , the Raychaudhuri equation is with τ an affine parameter and R αβ the Ricci tensor. The acceleration a α , expansion θ and shear σ αβ of the fluid are given by The rotation tensor ω αβ is defined as For a spatially homogeneous spacetime as chosen in the present manuscript, the equation simply leads to In general, to realize a collapsing nature dθ dτ must have a negative evolution, eventually reaching −∞ when the singularity forms. Such a criterion is realized only when It is interesting to note that the value of m makes no contribution in the predictability of the end state, as the signature of dθ dτ depends only on the value of a 0 . Therefore whether the collapsing axion dark matter keeps on collapsing until a singularity or bounces indefinitely dispersing away all the field distribution to zero value, depends on the parameter a 0 . This makes a critical parameter of the system, hinting at an underlying critical phenomena.
We note that before Eq. (20), we have no restriction over the equation of state parameter, as in what kind of matter makes ρ m and p m . However, carefully looking into the equation, if the matter distribution is another scalar field, then the left-hand-side (ρ m + p m ) is always positive. In such a case a positive right-hand-side is the standard understanding, leading to a collapsing geometry. However, a negative RHS leads to a non-trivial situation and can be realized if and only if the additional scalar field describing the matter distribution is a so-called ghost scalar field.

Non-minimally coupled axion
While the previous section dealt with a simple minimally coupled scalar field, one can expect that further generalizations can play more interesting roles under strong gravity. This is precisely where non-minimally coupled theories come in the picture by introducing a direct interaction between scalar field and gravity [79,80] leading to a generalization of the Lagrangian. The setup enjoys motivations from various avenues. However, the most important motivation in the present context is that the non-minimal coupling can be thought of as a renormalization requirement for quantum corrected scalar field theories in a curved background [81][82][83][84]. A non-minimal coupling of scalar and geometry finds motivation in superstring theory [85] as well. Cosmologically non-minimal generalzations carry greater importance as minimally coupled interacting scalar field models are increasingly being disfavoured by observations [86,87]. Ideally, the non-minimal interaction profile must be reconstructed or estimated from the observations or from a more fundamental theory, however, the case of a quadratic interaction has stood out amongst many other considerations [88][89][90]. For a properly chosen interaction profile, a non-minimally coupled scalar field can drive an early inflation alongwith the late time acceleration [91][92][93][94][95][96][97][98][99][100][101][102][103]. The setup also plays significant role in particle physics in interacting Higgs scalar models [104][105][106].
Motivated by these, in this section we consider an axion field directly interacting with geometry, apart from having a self-interaction. In that sense it is a bit further generalization from a minimally coupled axion. We consider an action The field equations for a spatially homogeneous metric Eq. (4) are written as A variation with respect to the scalar field φ allows us to write the axion evolution equation as We choose a particular form of the non-minimal coupling function U (φ) at the outset as The axion potential is already given by Eq. (3). Together with these the axion evolution equation becomes We follow a similar path to extract a solution, as we did in the previous section. We assume that the axion evolution equation is integrable and employ the integrability analysis. For brevity we do not repeat the detailed procedure. The integrability criterion produces a second order differential equation of a(t) This produces a solution for a(t) as We note that this solution is identical with the minimally coupled axion case, as can be seen comparing with Eq. (17). The only change involves the mass parameter m which is replaced in the non-minimal axion solution by a rescaled mass parameter, including the non-minimal coupling parameter U 0 , as m 2 = p 2 (1 + 6U 0 ). This only scales the time evolution of the stellar evolution and the axion scalar field while the qualitative behavior remains exactly the same compared to the minimally coupled case. The evolutions are shown in Figs. 5 and 6. Evolution of the two-sphere in Fig. 5 shows the radius of the two sphere as a function of time for p = 1 √ 2 and U 0 = 2. Evolution of the axion field is shown in Fig. 6 for a similar choice of initial conditions. In both of the graphs, the curve on top is for a 0 > 0 and the curve on the bottom is for a 0 < 0.
Using field Eq. (26) for density and Eq. (27) for pressure we can study the nature of the Null Energy Conditions writing Although not as simple as the expression for Null Energy Condition for the minimally coupled case, we can assess the and U 0 = 2. For the graph on top, a 0 > 0 (the blue curve is for a 0 = 1, the yellow curve is for a 0 = 2 and the green curve is for a 0 = 3). For the graph on bottom, a 0 < 0 (the blue curve is for a 0 = −1, the yellow curve is for a 0 = −2 and the green curve is for a 0 = −3) nature of the collapsing distribution overall from this, considering some approximations. For a collapsing case, i.e., for all a 0 > 0, as the collapse evolves towards final phase, the radius of two sphere a(t) → 0 and φ,φ → ∞ at a finite time. The rate of collapse increases and this makesȧ 2 a sharply increasing and dominating function of time compared toä. Therefore, as one approaches the zero proper volume singularity, For collapseȧ < 0 and this makes the second and third term on the RHS negative. However, the dominating nature of φ 2ȧ 2 a 2 is expected to make the RHS positive, satisfying the energy condition.
However, for a collapse and dispersal scenario,ȧ > 0 and φ,φ → 0 during the dispersal. This allows us to write the energy condition in the form where we have chosen t 0 = 1. For a collapse and dispersal scenario, a 0 < 0, and all the other terms on the RHS of the above equation are strictly positive. Therefore, the dispersal of the axion field and perfect fluid system is usually accompanied by a violation of the NEC. We comment that, since the axion field distribution almost becomes negligible in this case, the violation of energy condition can be simply an artifact of just the perfect fluid distribution. Usually a dark energy fluid is thought to be violating energy conditions, and from that purview this opens up some intriguing allied questions regarding a clustering dark energy that can coexist with axion dark matter. Besides, the energy density remains non-negative which means the system is up for any relevant quantization scheme. Moreover we recall that there are examples of non-minimally coupled scalar field systems violating some or more of the energy conditions [107].

Axion interacting with normal matter distribution
In this section, we work on an idea that under strong gravity, the self-interacting axion field can interact with ordinary mat-ter, which finds it's motivation from the string inspired dilaton gravity. A similar configuration of scalar field is already popular in a cosmological context, known as Chameleon Fields who can describe a smooth deceleration to acceleration transition of the universe and satisfy the so-called fifth force constraints [108,109]. These models also inspire ideas of interaction between dark energy and dark matter [110,111]. A time evolving solution in such a mathematical setup is our primary focus that can reveal the nature of the scalar field during gravitational collapse. We write an action where the axion scalar field has a non-minimal coupling with Ricci scalar as well as with the ordinary fluid. The action is similar to the Brans-Dicke theories [112], a prototype of scalar-tensor theories. The generalized action in Jordan frame is whereR is the Ricci scalar, φ is a scalar field, ω bd is the Brans-Dicke parameter, V (φ) is the self-interaction potential and f (φ) is an analytic function of the scalar field. We study the system in Einstein frame and employ a conformal transformation to transform the action. We define the transformation as with = √ Gφ. We redefine the scalar field as φ 0 ∼ G −1 , φ > 0 and ω bd > − 3 2 to have real solutions. Therefore in Einstein Frame the action in Eq. (36) becomes where σ = 8 π 2ω bd +3 . The self-interaction potential in Einstein frame is therefore written as We now focus entirely on the system in Einstein frame as in Eq. (39) from this point onwards. We do not think further about it's origin in an extended BD theory sense. We treat ψ as the axion scalar field. The self-interaction U (ψ) of ψ is the axion potential, which we choose to be as in Eq. (3). We also define h(ψ) = e −σ ψ/M p f (ψ). Variation of the action (39) with respect to the metric and the axion field gives the following equations and The energy momentum contributions of the axion field and the ordinary fluid are (43) and We write differentiation with respect to ψ as a prime. A calculation of ∇ μ G μν using Eq. (41) helps us write where T m = g μν T m μν . This means, covariant derivative of T m μν does not vanish, indicating an exchange of energy between matter and the axion scalar field, depending on the choice of matter Lagrangian density L m . We take a perfect fluid energy-momentum tensor. We also assume that it has an equation of state parameter w m .
u μ is a comoving four velocity. The perfect fluid Lagrangian density in the context of GR can be chosen as L m = p m or L m = −ρ m . Whether or not these two choices are equivalent is an intriguing question and has received significant attention in literature over time [113][114][115][116][117][118][119]. We choose L m = −ρ m and write the field equations and the axion evolution equation using Eqs. (41) and (42) for metric Eq. (4) and To realize the underlying physics we also write the matter conservation equation aṡ We integrate Eq. (50) once and write, The function is defined as ρ 0 is an integration constant. Interestingly, the Eqs. (51) and (52) indicate that the energy density evolves in a non-trivial manner as compared to GR, due to scalar interaction. For < 0, we can say that there is matter creation and energy is fed into matter. Similarly, for > 0 it can be thought that matter is annihilated. There is an outwards energy transfer from the matter.
The axion evolution equation Eq. (49) has a contribution from the non-trivial scalar-matter interaction, which is evident from the terms ρ m and h (ψ) on the LHS. For simplicity we assume the ordinary matter distribution inside the collapsing distribution to be dust, i.e., w m = 0. A careful look into Eq. (52) reveals that, during the final stages of a collapse, a → 0, implies βdψ ln a → 0. This implies that ρ m ∝ 1 a 3 during the final phases of the collapse. Using this simplification restricts the allowed time evolution, nevertheless, it allows us to solve the set of equations to extract a solution of significant importance. We assume the axion interaction term to be h(ψ) = ψ 0 ψ 2 2 and assume that the axion evolution equation is integrable. Moreover, we write the axion self-interaction potential as in Eq. (3). With all these, the axion evolution equation finally becomes The integrability criterion produces a second order differential equation of a(t) as follows The first integral from the equation is written aṡ An exact solution of this equation can not be written in closed form, however, we note that during the final stages of the collapse, a(t) → 0 and thus 2ρ 0 ψ 0 a + a 0 a 2 clearly dominates over m 2 a 2 2 . With this in mind, we can ignore the first term on the RHS and solve the equation to write The evolution of a(t) is shown in Figs. 7 and 8 as a function of time. In Fig. 7, the evolution is shown for different choices of ρ 0 (top graph) and ψ 0 (bottom graph) while all other parameters ar kept fixed. On the other hand, in Fig.  8, the plot is for different choices of the parameter a 0 . The graph on top shows the evolution for all a 0 > 0 and the graph on bottom shows the evolution for a 0 < 0 while all the other parameters ar kept fixed. It is interesting to see that the radius of two-sphere a(t) never reaches a zero proper volume for any combination of choices of initial conditions. The collapse always halts after a finite time, at a non-zero minimum cutoff value, after which it bounces.
We study the evolution of the axion scalar as a function of time using Eq. (7) and the evolution equation Eq. (54). The same solution is found by re-inserting the solution in Eq. (57) in Eq. (54) for a consistency check. The solution is complicated to write in a closed form. Therefore we resort to numerical study and plot the evolutions in Figs. 9 and 10. The initial conditions chosen for the numerical solutions are the values of a(t) and da(t) dt for an initial time t i . The graphs on top are for initial condition da(t) dt t=t i > 0 and the graphs on bottom are for initial condition da(t) dt t=t i < 0. It is easy to note that for different choice of these initial condition, the evolution of the scalar field only flips in signature, keeping the qualitative behavior same. From Figs. 9 and 10 we note that the axion field initially goes through a sharp increasing/decreasing time evolution, depending on the choice of da(t) dt t=t i . Just about the time the scale factor changes nature from collapsing to bouncing, the magnitude of the scalar field falls off sharply as a function of time and the field starts showing an oscillatory profile. The amplitude of the oscillations eventually dies down and the axion field disperses away to zero value at an asymptotic future. In Figs. 9 and 10, the factor determining the timescale of oscillation is an interesting arena of further research. Although we do not go into exclusive details of those aspects, we naively comment that it may be determined byȧ a .
Using the field Eqs. (47) and (48) we write the NEC as The axion field ψ rapidly decays throughout the collapse and therefore the role ofψ 2 eventually becomes less dominant compared to the first term of the numerator. We note that the choice of ψ 0 , indicates the choice of the axion-ordinary matter interaction, and plays a crucial part in determining the positive or negative signature of the LHS, and therefore, the validity of the NEC.

Matching with an exterior Vaidya spacetime
A collapsing stellar distribution must be surrounded by a suitably chosen exterior spacetime geometry. It is not unphysical to think that the exterior is almost vacuum. For a spherically symmetric case, a Schwarzschild solution is the most popular choice to fit in as the exterior and the metric and the extrinsic curvature is matched across the boundary hypersurface following the standard Israel-Durmois methodology [120][121][122]. With a scalar field however, this is less straightforward, as a Schwarzschild solution supports no scalar field [73]. We find it more reasonable to match the interior solutions in the present case with an exterior Vaidya metric across a boundary hypersurface [123,124]. We write the interior metric as The exterior Vaidya metric is In all the three examples discussed in the present manuscript, the basic metric structure of the interior follows the structure of Eq. (59), only the time dependence of a(t) changes. Therefore we deal with the matching condition for a general a(t). Continuity of metric or the first fundamental form gives and Eq. (62) is the first matching condition. Continuity of the extrinsic curvature on the boundary hypersurface gives We write, combining Eqs. (61), (62) and (63) dv dt = 3ra(t) 2 − r 2 3(ra(t) 2 − 2Ma(t)) . Equation (64) is the second matching condition. We write from Eq. (63) Continuity of extrinsic curvature leads to the derivative of Equations (65) and (66) are the third and fourth matching conditions. The three different cases presented in the manuscript are three specific cases of this general conditions.

Limitations of the models
In this section we mention two limitations of the models discussed in this manuscript. The limitations, although quite important, do not rule out the mathematical setup or the physics explored.
The first limitation comes in the form of violation of energy conditions. The energy conditions are an essential tool in General Relativity to ensure a concept of positive energy density. For a particularly chosen system, different energy conditions can be written depending on the energy momentum distribution under consideration. The null energy condition amongst these is particularly significant, as it was originally believed that it is not possible to construct a NEC violating system. This is more emphasized from the fact that NEC is a necessary requirement for the famous singularity theorem of GR to hold [4]. The theorem says that for a spherically symmetric homogeneous system, an NEC obeying fluid can not avoid a formation of singularity if it is initially collapsing. This is also realized using the Raychaudhuri Equation, as is discussed in this manuscript [74][75][76][77], atleast for most of the initial conditions. However, there comes a cosmological requirement as for a spatially flat cosmological metric (as in this manuscript, the only difference being that we are interested in a spatially homogeneous star, not the universe) the Hubble parameter can not grow if the NEC is to be satisfied. This is in absolute contradiction with observations and it is only natural to think that the so-called dark energy or the exotic fluid must violate the NEC. It is quite a familiar practice to consider this exotic fluid in the form of one or more interacting scalar fields. In order to accomodate for a proper theory that allows a violation of NEC, interacting scalar field theories demand particular attention. As discussed in literature [125], minimally coupled scalar theories generally produce homogeneous solutions obeying the NEC, with limited exceptions. One way to get around this is to consider vector fields [126,127] or higher-derivative Lagrangians, for instance, the Horndeski theories and it's recent generalizations [128][129][130][131][132][133][134][135][136]. It is also possible to describe an NEC violating system by playing with the idea of covariant energymomentum conservation, which has also been discussed in Sec. 5.
The second limitation comes accompanying the method of integrability of classical anharmonic oscillator system. In order to write the Klein-Gordon equation in the form of the integrability Eq. (5), it is assumed at the outset that the axion potential in Eq. (1) can be series expanded and written as in Eq. (3). It works as long as f >> φ, however, it can not be denied that we can not accomodate all terms by truncating the potential. As the scalar field collapses and density grows, after a certain point the sextic and higher order terms in the series expanded trigonometric potential becomes significant. Therefore, truncating the potential to quartic order no longer remains a good approximation. One way to amend this limitation is to try and look for numerical solution of the system of equations for the full cosine potential. On the other hand, writing a solution of the sytem purely from the boundary matching conditions and using the field equations as constraints can be another interesting prospect.

Conclusion
The manuscript deals with the subject of stellar collapse under the scope of classical gravity. The topic has received rigorous and recurring attempts over the years for possible explanations of some of the inconclusive questions of gravitational physics. For example, the correct time evolution of a collapsing stellar distribution or the predictability of a collapse end-state has baffled physicists. The main difficulty remains the high non-linearity of the field equations from which an extraction of time-dependent solutions is never a cakewalk. Therefore, we find it sensible to try and understand the dynamics in small and gradual steps, using relatively simpler setups that can describe the underlying physics. We consider the gravitational collapse of an interacting timedependent scalar field distribution alongwith a perfect fluid and gradually generalize and include non-trivial interactions of the scalar with geometry and even ordinary matter.
To read through the non-linearity, we point transform the Klein-Gordon type equation for the scalar field evolution into a totally integrable form using the Euler-Duarte-Moreira analysis for anharmonic oscillator equation systems [68][69][70]. Application of this method leads us to some new interesting solutions of the radius of the two-sphere for the metric tensor and the scalar field. We deal with a simple, spatially homogeneous Oppenheimer-Snyder type metric. For a minimally coupled self-interacting scalar field, we find an equally likely probability that the collapse can form a singular end-state or a bounce, depending on an initial collapsing parameter. The parameter comes out as a constant of integration and can be related to the initial value of the radius of two-sphere or the initial volume of the stellar distribution (the parameter a 0 ). We predict this parameter to be a critical parameter and catch a hint of an unexplored critical phenomena of the system. We predict a similar end result for a second example when the scalar field interacts with the geometry through a non-minimal coupling. In a third example, the scalar field interacts with geometry as well as the normal matter distribution through a non-minimal coupling function. For no functional form of the self-interaction potential or the nonminimal interaction functions or for no value of the parameters of the metric coefficients, any zero proper volume singularity is reached. This serves as an example of a rare case of non-singular gravitational collapse, where the collapse can only go unhindered until a minimum accessible volume of the spherical geometry. After this point, the evolution changes nature and bounces indefinitely. The scalar field has a periodic evolution as a function of time with changing amplitude and frequency. Eventually the scalar field disperses away all of it's field strength and asymptotically reaches zero.
We note that, if and when a singular end state develops, the visibility of the same is dependent on the existence of nonspacelike trajectories coming out of the singular epoch and reaching a faraway observer. In the present case, for all singularity reaching solutions, the time at which a zero proper volume is reached is independent of the radial coordinate. This indicates that all the collapsing shells consisting of the scalar field and perfect fluid fall onto the singularity simultaneously. This can only form a covered singularity or a Black Hole [137].
The investigation stems from an idea to study the evolution of the QCD inspired 'axion' field under strong gravity. The axion field can actually serve as a very good fit for the cosmological 'dark' matter distribution of the universe, the driver of cosmic deceleration. The scalar field and the self-interaction potential used throughout this manuscript describes the evolution of axion dark matter in spherical symmetry. Our discussion effectively leads one to realize the possibility of axion dark matter to collapse and form singularities, gravitationally bound objects or explode and disperse away into zero field strength, under strong gravity. This result is not very far from the outcomes predicted by a few existing novel works on axion collapse. For instance, it was shown quite recently that a collapsing axion distribution can not collapse to zero volume due to the repulsive higher order terms of the trigonometric potential. It has also been proved that energy emission through heat flux is actually quite normal during the end stage of such a collapse, leading rapid to emission of axions [56,138]. In the present manuscript, we find a dependence on initial conditions for the final outcome. (a) 'Axion Black Holes', where a zero proper volume singularity forms in the end. Such a singularity is always accompanied by the formation of an apparent horizon as it is supposed to be for a spatially homogeneous geometry. In this case, the initial conditions ensure that the Null Energy Condition is satisfied. (b) 'A Collapse and Dispersal of axion field', where the collapse ends at a non zero minimum radius of two-sphere and the stellar distribution begins to bounce indefinitely, dispersing away all the axion scalar field strength. This scenario is always accompanied by a violation of null energy condition. The second case can be roughly compared with the phenomena of a stellar distribution radiating energy away from itself through bursts of axions, mimicking the 'Bosenova phenomena' of cold-atom physics [139]. The accompanying perfect fluid distribution in the models can somewhat be connected with a Dark Energy fluid as we show that even for perfectly reasonable physical conditions, it is natural for this fluid to violate atleast the Null Energy Conditions. This serves as a hint towards a couple of scenarios. (a) A collapsing fluid violating energy conditions under any circumstances is suggestive of a possible clustering nature of dark energy. (b) A stellar distribution involving a Scalar dark matter and a fluid dark energy combination under strong gravity can lead one to understand the nature of gravitational interaction in a cosmic era where both of these dark components were dominant and interacting with one another.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .