Lie Symmetries, Painlevé analysis and global dynamics for the temporal equation of radiating stars

Ordinary differential equations play an important role in the study of physical systems. There are various approaches for the study of the properties of physical systems described by differential equations as well as to construct and determine exact and analytic solutions. Relativistic astrophysics is an active area employing various solution-generating techniques to solve highly nonlinear systems of differential equations. From finding exact solutions of the Einstein field equations or their modifications (Einstein-Gauss Bonnet gravity, $f(R)$ gravity or Brans-Dicke theory, etc.) describing stellar objects through to the study of causal thermodynamics, where solutions of governing differential equations have played a key role.

The end-states of gravitational collapse of bounded configurations have held the attention of astrophysicists since the pioneering work of Oppenheimer and Snyder¹ in which they studied idealised collapse of a dust sphere. The Weak Cosmic Censorship Conjecture, first articulated by Penrose, forbids the existence of naked singularities arising from continued gravitational collapse. However, there have been many counter-examples put forward within the framework of Einstein’s classical general relativity ^{2, 3, 4, 5}. The confirmation of classical general relativity as a cornerstone of gravitational theory was borne out in 2019 when the photograph of the shadow of a black hole ⁶ was obtained, heralding a new frontier of theoretical predictions and observations. The discovery of the Vaidya solution ⁷ paved the way for researchers to study dissipative collapse of stars. The boundary of the collapsing object divides spacetime into two distinct regions, ${\cal M^{-}}$, the interior region and ${\cal M^{+}}$, the exterior region described by the Vaidya solution. Early work on dissipative gravitational collapse can be attributed to Herrera and co-workers (see ^{8, 9} and references therein). In their investigations they studied spherically symmetric, shear-free stellar objects undergoing dissipative gravitational collapse in the form of a radial heat flux. The junction conditions for the smooth matching of the interior spacetime to the exterior Vaidya solution were derived by Santos ¹⁰. The Santos junction conditions demonstrated that the pressure at the boundary of a collapsing, radiating star is nonzero and is proportional to the magnitude of the outgoing radial heat flux. This junction condition represents the conservation of momentum across the boundary of the collapsing star. Recently, the Santos junction conditions have been extended to include a dynamically unstable core with a general energy-momentum tensor describing an imperfect fluid with heat flux and null radiation with the exterior being described by the generalised Vaidya solution ¹¹. Over the next two decades, the study of radiating stars has provided us with a rich insight into the end-states of continued gravitational collapse, particularly with regards to time of formation of the horizon¹², temperature profiles and relaxational effects related to causal heat flux¹³. The shear-free models were subsequently extended to include the effects of shear viscosity. It has been demonstrated that the Chandrasekhar stability criterion for isotropic fluid spheres is modified in the presence of shear viscosity. Furthermore, in both the Newtonian and relativistic regimes, the shear viscosity decreases the instability of the stellar fluid¹⁴.

The inclusion of shear, anisotropy, electromagnetic field, and rotation in the slow approximation have been fruitful in studying the thermodynamics of such systems. The impact of shear on the kinematics and dynamics of the collapse process has been addressed by several authors. The instability of the shear-free condition has been demonstrated by Herrera et al. ¹⁵ in which they showed that the shear-free condition may hold for a limited epoch of the collapse process. The presence of pressure anisotropies, density inhomogeneities, and dissipative fluxes can mimic shear-like effects ¹⁶. The inclusion of shear during dissipative collapse has led to interesting results when contrasted to the shear-free case. Shearing effects lead to higher core temperatures it has been shown that horizon formation is delayed when shear is present. In an attempt to model shearing, radiating stars, Ivanov introduced the so-called horizon function which simplifies the boundary condition representing the temporal behavior of the model^{17, 18}. The horizon function has a physical attribute in the sense that it is directly related to the surface redshift of the collapsing sphere. Once this function is determined from the boundary condition, the end-state or possible outcome of continued gravitational collapse can be studied. The avoidance of the singularity was demonstrated by using a simple model of shear-free collapse in which the rate of collapse is balanced by the rate of energy emission to the exterior spacetime ¹⁹. This so-called horizon-free model, or Banerjee, Chatterjee & Dadhich (BCD) ¹⁹ model, forms the basis of the work contained in this paper. The horizon-free collapse in the presence of shear was studied in various models in which the gravitational potentials were highly simplified. The so-called Euclidean stars, in which the areal radius is equal to the proper radius were shown to undergo horizon-free collapse ²⁰. The BCD model has two inherent assumptions: the interior spacetime is shear-free and the gravitational potentials are separable in space and time. Works by Chan and co-workers attempted to generalise the shear-free model to include shear by demanding that the metric functions be separable. These models have several drawbacks such as the proper radius being independent of time or the resulting boundary condition rendered solvable via numerical methods^{21, 22}.

In this piece of work, we investigate the dynamics for the master equation of the temporal equation of radiating stars. The equation that we are interested describes the Einstein field equations for a spherically symmetric line element

for a pressure isotropy fluid with an exterior solution the Vaidya metric ⁷ such that to describe a radiative solution. Banerjee et al. ²³ suggested the metric ansatz $A\left(t,r\right)=1+\zeta_{0}r^{2},~{}B\left(t,r\right)=R(t)~{}$where $\zeta_{0}$ is a positive constant. In this case, the field equations reduce to the second-order differential equation

where $\alpha$ and $\beta$ are constants and the dot means derivative with respect to the time $t$. A special solution of the latter equation is the solution $R=-Ct$ where $C>0$ is a constant, which describes a collapsing star ²³. Recently, in ²⁴ new families of exact solutions for the master equation (2) were found using Lie symmetries. This body of work forms part of several investigations about dissipative collapse via Lie symmetries. Radiating stars with shear, in which the particle trajectories within the stellar fluid were geodesics, were studied using Lie symmetries²⁵. Several new solutions were obtained while other solutions were reduced to well-known cases studied earlier in the literature. The geodesic case with shear was extended to the most general shearing matter distribution radiating energy to the exterior spacetime. The method of Lie symmetries proved useful in obtaining four new classes of solutions ²⁶, two of which included the horizon function and the Euclidean condition.

In the following, we study the dynamics and the stability properties for the master equation (2) and also, we show in detail how the method of Lie point symmetries and the singularity analysis can be used for the derivation of new analytic and exact solution for equation (2).

The theory of Lie symmetries ²⁷ is a powerful approach for the investigation of invariant functions for differential equations. The novelty of Lie theory is that the infinitesimal representations of the finite transformations of continuous groups are considered, by moving from the group to a local algebraic representation, and to studying the invariance properties under them. The resulting invariant functions can be used for the reduction of order for a given ordinary differential equation and consequently for the derivation of solutions. These exact solutions which follow from the application of Lie symmetries are known as similarity solutions. The theory of Lie symmetries cover a range of applications in physical science and gravitational theory, see for instance ^{28, 29, 30, 31, 32, 33} and references therein.

The singularity analysis, which is mainly associated with the school of Painlevé, that is why also it is called Painlevé analysis ³⁴, is another powerful mathematical approach for the derivation of solutions of differential equations. When a differential equation passes the singularity the analysis we say that the differential equation possesses the Painlevé property. Nowadays the application of the singularity analysis is summarised in ARS algorithm, from the initials of Ablowitz, Ramani, and Segur ^{35, 36, 37} who established a systematic method for investigation of analytic solutions, inspired by the approach applied by Kowalevskaya ³⁸ for the determination of the third integrable case of Euler’s equations for a spinning top. The basic characteristic for a differential equation to possesses the Painlevé property is the existence of at least a movable singularity. In the singularity analysis the analytic solution for a given differential equation is expressed by Painlevé Series, and specifically in our consideration we shall write the analytic solution in the Puiseux series.

On the other hand, stability analysis provides us with an important tool to investigate the evolution of the given dynamical system. Indeed, from the stability analysis, we can investigate if a given solution is stable or not, while we can determine families of initial conditions such that specific behavior for the dynamical system to be stable. In addition, we can define constraints for the free parameters of the differential equations according to the stability of exact solutions. Hence, we can extract important information about the nature of the free parameters. In our consideration, for equation (2)  we can understand how the free parameters $\alpha$ and $\beta$ affects the dynamics. The plan of the paper is as follows.

In Section 2 we investigate the stability properties for the master equation (2) for the exact solution found in ²³, while we perform a detailed study of the global dynamics for the master equation. From the dynamical analysis, a new family of exact solutions is constructed. In Section 3, we investigate equation (2) by applying Lie’s theory, as we study if equation (2) possesses the Painlevé property. We find that the differential equation admits two Lie point symmetries which form the $A_{2,2}$ Lie algebra. Consequently, the application of the Lie symmetries provides two similarity transformations. Moreover, the application of the ARS algorithm indicates that the master equation (2) satisfies the Painlevé test, and the analytic solution is expressed by a right Painlevé series. However, for a specific relation between the two free variables $\alpha$, $\beta$; a new exact solution is determined. Finally, in Section 4 we discuss our results and draw our conclusions.

For the master equation (2),  by convenience we assume $\alpha>0$. A special and simple solution to this equation is $R=-Ct$ defined for $-\infty<t\leq 0$, where $C>0$ is a constant given by the positive roots of

That is, for $\beta>0$, we have only one positive root $C=\frac{1}{2}\left(\sqrt{\alpha^{2}+4\beta}+\alpha\right)$. For $\beta<0$ and $\alpha^{2}+4\beta>0$, we have two different positive roots: $C_{-}=\frac{1}{2}\left(\alpha-\sqrt{\alpha^{2}+4\beta}\right)$ and $\quad C_{+}=\frac{1}{2}\left(\sqrt{\alpha^{2}+4\beta}+\alpha\right)$.

To avoid ambiguities we prefer to use the relation $\beta=C^{2}-C\alpha$, and in the analysis separate the cases $\beta<0$ and $\beta\geq 0$.

For the analysis of stability of the scaling solution $R(t)=-Ct$, with $C=\frac{1}{2}\left(\sqrt{\alpha^{2}+4\beta}+\alpha\right)>0$ in the interval $-\infty<t<0$ we use similar methods as in Liddle & Scherrer ³⁹ and Uzan ⁴⁰. For this purpose, the new logarithmic time variable $\tau$ related to $t$ through

such that $t\rightarrow-\infty$ as $\tau\rightarrow-\infty$ and $t\rightarrow 0$ as $\tau\rightarrow+\infty$ is defined.

Then, replacing the relation between $t$ and $\tau$ given by (4) in the expression $R(t)$, we obtain a function of $\tau$ from $R(t)$ denoted and defined by

Using again the transform (4), the scaling solution $R(t)=-Ct$ can be re-expressed as a function of $\tau$ as $R_{s}(\tau)=Ce^{-\tau}$.

To compare a generic solution $\bar{R}(\tau)$ with the scaling solution $R_{s}(\tau)=Ce^{-\tau}$, we define a dimensionless function of $\tau$ as the ratio

Thus, the solution $R_{s}(\tau)$ will corresponds to the equilibrium point $u=1$ in the new formulation of the master equation (2). The physical region corresponds to $u\geq 0$. Recall, the physical solution $R=-Ct$ is defined for $-\infty<t\leq 0$, where $C>0$ is a constant.

To reformulate the master equation (2) in terms of $u(\tau)$, and its derivatives with respect to $\tau$, we use the chain rule to obtain several differentiation rules. Let be denoted the time derivative with respect to $\tau$ be denoted by a prime and the time derivative with respect to $t$ by a dot. That is, for $f(\tau)$, let

and for $g(t)$, let

Then,

and

Therefore, by using the rules (9), and substituting expression (4), the equation (2) becomes

The next step is to rewrite equation (11) in terms of $u(\tau)$ and its derivatives.

Solving equation (6) for $\bar{R}$, and then using the derivatives rules for products of functions and the chain rule, we obtain, subsequently,

Substituting (12), (13) and (14) in equation (11), and dividing by the overall factor $e^{-\tau}$, we obtain

where we have used the relation $\beta=C^{2}-C\alpha$.

Defining the new function

we obtain the dynamical system

The scaling solution $R_{s}(\tau)=Ce^{-\tau}$ corresponds to the equilibrium point $u=1,v=0$.

Defining $\varepsilon$ through $u=1+\varepsilon$, we obtain the final dynamical system

where the equilibrium point is translated to the origin.

The linearisation matrix is defined by

Then, linearising around the equilibrium point, $\varepsilon=0,v=0$, we obtain

The linearisation matrix

has eigenvalues $\left\{1,1-\frac{\alpha}{2C}\right\}$.

Assume first $\beta\geq 0$. In this case, the origin is always unstable as $\tau\rightarrow+\infty$ due to $2C=\alpha+\sqrt{\alpha^{2}+4\beta}>\alpha$. That is, the origin is stable as $\tau\rightarrow-\infty$.

An additional equilibrium point is

Evaluating the linearisation matrix $J\left(\frac{\alpha}{C}-2,0\right)$, we obtain the eigenvalues $\left\{1,\frac{2C-\alpha}{2(C-\alpha)}\right\}$. Due to $2C-\alpha>0$ it is unstable as $\tau\rightarrow\infty$. Indeed for $0<\frac{\alpha}{2}<C<\alpha$ it is a saddle, whereas for $C>\alpha>0$ is an unstable node. The last conditions is forbidden due to the physical condition $u\geq 0$ evaluated at $(u,v)=\left(\frac{\alpha}{C}-1,0\right)$ implies $\alpha\geq C$.

Now, let us study the case $\beta<0,\alpha<-2\sqrt{-\beta}$ or $\beta<0,\alpha>2\sqrt{-\beta}$. Henceforth, we have two solutions

where

Observe that $2C_{+}>\alpha$, implies that $R_{s+}(\tau)$ is an unstable solution (unstable node) as $\tau\rightarrow\infty$. Due to $\alpha>2C_{-}>0$, $R_{s-}(\tau)$ is an unstable (saddle) solution.

Figure 1 shows the phase plot of system (17), (18) for some $\alpha>0,\beta>0$ and the choice $2C:=\alpha+\sqrt{\alpha^{2}+4\beta}$. The origin represents the solution $R_{s}(\tau)=Ce^{-\tau}$, which is unstable (node).

Figure 2 shows the phase plot of system (17), (18) for some $\alpha>0,\beta<0$ and the choice $2C_{-}:=\alpha-\sqrt{\alpha^{2}-4|\beta|}$. The origin represents the solution $R_{s-}(\tau)=C_{-}e^{-\tau}$, which is unstable (saddle).

Figure 3 shows the phase plot of system (17), (18) for some $\alpha>0,\beta<0,\alpha^{2}+4\beta\geq 0$ and the choice $2C_{+}:=\alpha+\sqrt{\alpha^{2}-4|\beta|}$. The origin represents the solution $R_{s+}(\tau)=C_{+}e^{-\tau}$, which is unstable (node).

As can be seen in figures 1, 2 and 3, there are non-trivial dynamics as $(u,v)$ are unbounded.

In the following, we discuss the application of Lie’s theory and the singularity analysis for the derivation of exact and analytic solutions for the master equation (2).