Charged stellar model with three layers

Formulation of stellar models in a core-envelope setting has become an area of interest to many researchers in astrophysics and other related subjects. The physical phenomena in such models are still complicated in terms of the solutions generated. The reason is that the interior structure of the stellar body has several layers with distinct physical properties. It is a challenging question for researchers to identify the nature of the layers possessed by superdense stellar objects. The compact bodies with such properties have strong magnetic and gravitational fields. These include neutron stars, quark stars, gravastars, pulsars, white dwarfs and quasars. It has been noted that superdense stellar objects possess high pressures and extreme density leading to a variety of physical properties. It has been demonstrated that the overall masses and densities of superdense objects are approximately 1.0 – 2.0 M_⊙ and 10¹⁵ g cm⁻³ respectively (Pant et al. 2020; Bisht et al. 2021). According to Jasim et al. (2018), extremely superdense stellar bodies may alter the spacetime geometry while the radii and masses remain between 11.0 – 15.0 km and 1.4 – 2.0 M_⊙ respectively. A variety of mathematical approaches has been used to solve the system of field equations to investigate the physical features of superdense stellar objects (Tiwari 2010).

The concept of a stellar object with a core and an envelope is embedded in the general theory of relativity discovered by Albert Einstein between 1907 and 1916. The theory has been applied to describe the structure of stellar bodies and understanding the Universe at large (Misner et al. 1973). In addition to Einstein's theory, various theories of gravitation like vector-tensor theory and scalar-tensor theory were developed with some modifications to Einstein's theory. These theories may be utilized to formulate relativistic stellar models governed by a set of differential equations (Collins 2003). Stellar models may also be built with corrections from the Ricci and Riemann tensors in modified gravity.

The Einstein field equations are the most basic tools in the general theory of relativity. These equations depict the relationship between metric potentials, pressures and energy density in the form of differential equations (Stephani et al. 2003). The field equations are purely geometric. It is possible to solve the field equations if some reasonable physical assumptions are imposed. Importantly, the field equations are utilized to investigate different structures, behaviors and properties of stellar spheres relevant to astrophysical studies. They have drawn the attention and interest of researchers to develop relativistic stellar models in static or non-static spherically symmetric spacetime (Fulara & Sah 2018). A variety of mathematical approaches has been used to solve the system of field equations to investigate the physical features of superdense stellar objects (Tiwari 2010).

The electric field introduced in the Einstein-Maxwell field equations has to capture the reality of physical properties that charged stellar objects in the Universe should possess. Stellar spheres should loose electric charge in time to become uncharged. A realistic physical model for charged spheres should therefore have a vanishing electric field in some limit. The presence of an electric field in a stellar object creates repulsion with the pressure gradient to prevent gravitational collapse. This phenomenon makes the object stable and works against collapse (Jasim et al. 2018). The effect of electric field on the critical mass, stability and redshift in single-layered models has been described in Dev & Gleiser (2002); Ivanov (2002); Chaisi & Maharaj (2005); Chaisi & Maharaj (2006); Sunzu et al. (2014). The effect of the electric field on physical properties in two layered models has been discussed in the work developed by Mafa Takisa & Maharaj (2016). It is important to generate a three layer model with an electric field present for a better understanding of the model.

The importance of pressure anisotropy has been demonstrated in different stellar models. According to Malaver (2018), the presence of pressure anisotropy in the interior of the stellar objectmay be due to the existence of a solid core, a particular phase transition and electric charge. The detailed analysis of pressure anisotropy has also been highlighted in the studies performed by Bower & Liang (1974); Bijalwan (2011); Cosenza et al. (1981); Sokolov (1980); Usov (2004); Komathiraj & Maharaj (2007); Ivanov (2010); Maharaj et al. (2014); Maurya et al. (2018) and Sunzu et al. (2019).

Models with two layers describing a core layer and the corresponding envelope layer for superdense stellar objects have been formulated. This has been illustrated by Sharma & Mukherjee (2002) who formulated a coreenvelope model with an inner layer containing quark matter and the envelope layer which is less compact. The model which describes a parabolic density profile, developed in Negi et al. (1990), indicates the continuity of all variables at the interface between the core and the envelope. The approach developed by Thomas et al. (2005) describes a superdense matter configuration which is characterized by an isotropic fluid for the core and an anisotropic fluid at the envelope layer. Other models with two layers include the works of Mafa Takisa & Maharaj (2016); Mafa Takisa et al. (2019); Pant et al. (2019); Gedela et al. (2019); Hansraj et al. (2016); Paul & Tikekar (2005) and Tikekar & Jotania (2009). It is also important to consider the physics of a three-layered model with both electric field and anisotropy present.

Models with three interior layers are rarely found in general relativity. The recent models studied by Bisht et al. (2021) and Pant et al. (2020) indicate that a sublayer exists between the core-envelope boundary. Themodels analyzed the physical properties and features of a neutron star as a dense stellar object. However, the models generated were uncharged with the core obeying a linear equation of state (EoS), and intermediate and envelope layers satisfying a quadratic EoS. In the presentmodel, we include the electric field which is missing in Bisht et al. (2021) and Pant et al. (2020).

The motivation of this paper is to generate a new stellar model comprising three regions. We apply a separate EoS in each region depending on the nature of the material it contains. It is reasonable to have a stellar model that incorporates three interior layers with distinct EoSs. The core layer is described by a linear EoS sufficient to describe quark matter, and the intermediate and envelope layers both satisfy different quadratic EoSs. In our model, we choose one of the gravitational potentials and the electric fields. We generate a new anisotropic charged model containing three regions which generalizes the neutral models developed by Bisht et al. (2021) and Pant et al. (2020).

We describe the interior with a static spherically symmetric spacetime with metric given in Schwarzschild coordinates (x^j ) = (t, r, θ, ϕ) as

$$\begin{eqnarray}\begin{array}{ll}d{s}^{2}= & -{e}^{2\nu (r)}d{t}^{2}+{e}^{2\lambda (r)}d{r}^{2}\\ & +{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}),\end{array}\end{eqnarray} \tag{ 1 }$$

where ν(r) and λ(r) define the metric potentials. The exterior line element is given by the Reissner-Nordstrom spacetime as

$$\begin{eqnarray}\begin{array}{ll}d{s}^{2}= & -\left(\begin{array}{c}1-\displaystyle \frac{2M}{r}+\displaystyle \frac{{Q}^{2}}{{r}^{2}}\end{array}\right)d{t}^{2}\\ & +{\left(\begin{array}{c}1-\displaystyle \frac{2M}{r}+\displaystyle \frac{{Q}^{2}}{{r}^{2}}\end{array}\right)}^{-1}d{r}^{2}\\ & +{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}),\end{array}\end{eqnarray} \tag{ 2 }$$

where Q stands for the total charge and M represents the total mass of the sphere. For charged stellar objects, the energy momentum tensor is written as

$$\begin{eqnarray}\begin{array}{ll}{T}_{ij}= & {\rm{diag}}\bigg(-\rho -\displaystyle \frac{1}{2}{E}^{2},\,\,{p}_{r}-\displaystyle \frac{1}{2}{E}^{2},\,\,{p}_{t}+\displaystyle \frac{1}{2}{E}^{2},\\ & \,{p}_{t}+\displaystyle \frac{1}{2}{E}^{2}\bigg),\end{array}\end{eqnarray} \tag{ 3 }$$

where ρ defines the energy density, E is the electric field, and p_r and p_t are radial and tangential pressures respectively. Setting G = c = 1, the system of field equations becomes

$$\begin{eqnarray}\begin{array}{lll} & & 8\pi \rho +\displaystyle \frac{1}{2}{E}^{2}=\displaystyle \frac{1}{{r}^{2}}(1-{e}^{-2\lambda })+\displaystyle \frac{2{\lambda }^{\prime}}{r}{e}^{-2\lambda },\end{array}\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}\begin{array}{lll} & & 8\pi {p}_{r}-\displaystyle \frac{1}{2}{E}^{2}=-\displaystyle \frac{1}{{r}^{2}}(1-{e}^{-2\lambda })+\displaystyle \frac{2{\nu }^{\prime}}{r}{e}^{-2\lambda },\end{array}\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}\begin{array}{lll} & & 8\pi {p}_{t}+\displaystyle \frac{1}{2}{E}^{2}={e}^{-2\lambda }\left(\begin{array}{c}{\nu }{^{\prime}{\prime} }+{\nu }^{{^{\prime} }^{2}}-{\nu }^{\prime}{\lambda }^{\prime}+\displaystyle \frac{{\nu }^{\prime}}{r}-\displaystyle \frac{{\lambda }^{\prime}}{r}\end{array}\right)\\ & & \,\,\,\,+\displaystyle \frac{2{\lambda }^{\prime}}{r}{e}^{-2\lambda },\end{array}\end{eqnarray} \tag{ 4c }$$

$$\begin{eqnarray}\begin{array}{lll} & & \sigma =\displaystyle \frac{1}{{r}^{2}}{e}^{-\lambda }{({r}^{2}E)}^{\prime},\end{array}\end{eqnarray} \tag{ 4d }$$

where primes (') represent derivatives with respect to radial coordinate (r).

We introduce the transformation variables from Durgapal & Bannerji (1983) given by

$$\begin{eqnarray}&&x={r}^{2},\,Z(x)={e}^{-2\lambda },\,{e}^{2\nu }={A}^{2}{y}^{2}(x).\end{eqnarray} \tag{ 5 }$$

Using Equation (5), the system (4) becomes

$$\begin{eqnarray}\begin{array}{lll} & & 8\pi \rho =\displaystyle \frac{1-Z}{x}-2\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{2}{E}^{2},\end{array}\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}\begin{array}{lll} & & 8\pi {p}_{r}=-\displaystyle \frac{1}{x}(1-Z)+4Z\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}+\displaystyle \frac{1}{2}{E}^{2},\end{array}\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}\begin{array}{lll} & & 8\pi {p}_{t}=4xZ\displaystyle \frac{1}{y}\displaystyle \frac{{d}^{2}y}{d{x}^{2}}+\left(\begin{array}{c}4Z+2x\displaystyle \frac{dZ}{dx}\end{array}\right)\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}\\ & & \,\,+\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{2}{E}^{2},\end{array}\end{eqnarray} \tag{ 6c }$$

$$\begin{eqnarray}\begin{array}{lll} & & \Delta =8\pi {p}_{t}-8\pi {p}_{r}\\ & & \,\,\,=4xZ\displaystyle \frac{1}{y}\displaystyle \frac{{d}^{2}y}{d{x}^{2}}+\left(\begin{array}{c}4xZ+2x\displaystyle \frac{dZ}{dx}-4Z\end{array}\right)\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}\\ & & \,\,+\displaystyle \frac{dZ}{dx}+\displaystyle \frac{1}{x}(1-Z)-{E}^{2},\end{array}\end{eqnarray} \tag{ 6d }$$

$$\begin{eqnarray}\begin{array}{lll} & & \sigma =2\left(\begin{array}{c}x\displaystyle \frac{dE}{dx}+E\end{array}\right)\sqrt{\displaystyle \frac{Z}{x}}.\end{array}\end{eqnarray} \tag{ 6e }$$

The system (6) has eight unknown variables. To find a solution of the system (6), we need to specify any two parameters so as to find expressions for the others. In our model we choose to specify metric potential Z and electric field intensity E². We apply the metric potential Z used in the study of Pant et al. (2020). This is expressed in the form

$$\begin{eqnarray}&&Z(x)=\,{e}^{-2\lambda }=1-hx+g{x}^{2},\end{eqnarray} \tag{ 7 }$$

where h and g are arbitrary real constants. The metric function Z in Equation (7) is real, continuous and regular at the center which allow an extensive range of values of h and g. It is in acceptable condition to avoid the singularity at the center and within the stellar sphere. It is possible to solve the system (6) once the electric field is specified. The particular choice of electric field is specified in the form

$$\begin{eqnarray}&&{E}^{2}=\,mxZ=mx(1-hx+g{x}^{2}),\,m\ge 0,\end{eqnarray} \tag{ 8 }$$

where m stands for an arbitrary real constant. The choice of E² given in Equation (8) vanishes at the center (x = 0) of the stellar object, but remains real and continuous throughout the interior points of the stellar body. It is physically reasonable and, thus, consistent for the study of the charged stellar models. Substituting m = 0 in Equation (8), we regain an uncharged model in Pant et al. (2020).

We have three interior regions, namely: the core (η), the intermediate (ι) and the envelope (ε) as follows:

• Core layer (Region 1); 0 ≤ r ≤ R_η ,
• Intermediate layer (Region 2); R_η ≤ r ≤ R_ι and Envelope layer (Region 3); R_ι ≤ r ≤ R_ε .

For the classified boundaries, line element (1) becomes

$$\begin{eqnarray}&&d{s}^{2}{|}_{1}=-{e}^{2{\nu }_{\eta }}d{t}^{2}+{e}^{2{\lambda }_{\eta }}d{r}^{2}+{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}),\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&d{s}^{2}{|}_{2}=-{e}^{2{\nu }_{\iota }}d{t}^{2}+{e}^{2{\lambda }_{\iota }}d{r}^{2}+{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}),\end{eqnarray} \tag{ 9b }$$

$$\begin{eqnarray}&&d{s}^{2}{|}_{3}=-{e}^{2{\nu }_{\epsilon }}d{t}^{2}+{e}^{2{\lambda }_{\epsilon }}d{r}^{2}+{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}),\end{eqnarray} \tag{ 9c }$$

in the three regions.

3.1. The Core

We represent the core layer to be quark matter satisfying the linear EoS in the form

$$\begin{eqnarray}&&{p}_{{r}_{\eta }}=\alpha {\rho }_{\eta }-\beta,\end{eqnarray} \tag{ 10 }$$

where α and β are arbitrary real constants. Combining Equations (6a) and (10) yields

$$\begin{eqnarray}&&{p}_{{r}_{\eta }}=\alpha \left(\begin{array}{c}\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right)-\beta.\end{eqnarray} \tag{ 11 }$$

Equating Equation (6b) and Equation (11) gives

$$\begin{eqnarray}\begin{array}{ll}\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}= & (\alpha +1)\left(\begin{array}{c}\displaystyle \frac{1-Z}{4xZ}\end{array}\right)-\displaystyle \frac{\alpha }{2Z}\displaystyle \frac{dZ}{dx}\\ & -\displaystyle \frac{{E}^{2}(\alpha -1)}{8Z}-\displaystyle \frac{2\pi \beta }{Z}.\end{array}\end{eqnarray} \tag{ 12 }$$

Substituting Equation (7) and Equation (8) into Equation (12) yields

$$\begin{eqnarray}\begin{array}{ll}\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}= & [h(2-m{x}^{2}(\alpha -1)-6\alpha)mx(\alpha -1)\\ & \times gx(-2+m{x}^{2}(\alpha -1)+10\alpha)+16\pi \beta]\\ & \times {[8x(h-gx)-8]}^{-1}.\end{array}\end{eqnarray} \tag{ 13 }$$

Integrating Equation (13) gives

$$\begin{eqnarray}\begin{array}{ll}y & =\exp \bigg[\displaystyle \frac{1}{16}\bigg[-m{x}^{2}(\alpha -1)+4(h(\alpha -1))\\ & -16\pi \beta \arctan (2gx-h)(\sqrt{4g-h})\bigg]\times {\bigg[\sqrt{4g-{h}^{2}}\bigg]}^{-1}\\ & +\displaystyle \frac{1}{16}(2(1-5\alpha)\mathrm{log}(1-hx+g{x}^{2}))\\ & \times {(\sqrt{4g-{h}^{2}})}^{-1}+{k}_{0}\bigg].\end{array}\end{eqnarray} \tag{ 14 }$$

Putting Equations (7), (8), (13) and (14) into the system (6), we obtain

$$\begin{eqnarray}\begin{array}{lll} & & {e}^{2{\lambda }_{\eta }}={[1-hx+g{x}^{2}]}^{-1},\end{array}\end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray}\begin{array}{lll} & & {e}^{2{\nu }_{\eta }}={A}^{2}{[\exp ({Q}_{0}(x))]}^{2},\end{array}\end{eqnarray} \tag{ 15b }$$

$$\begin{eqnarray}\begin{array}{lll} & & rh{o}_{\eta }=[h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]\\ & & \,\,\times {[16\pi]}^{-1},\end{array}\end{eqnarray} \tag{ 15c }$$

$$\begin{eqnarray}\begin{array}{lll} & & {p}_{{r}_{\eta }}=\alpha [h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]\\ & & \,\,\times {[16\pi]}^{-1}-\beta,\end{array}\end{eqnarray} \tag{ 15d }$$

$$\begin{eqnarray}\begin{array}{lll} & & {p}_{{t}_{\eta }}={Q}_{1}(x)+{Q}_{2}(x)+{Q}_{3}(x)+{Q}_{4}(x)\\ & & \,\,+{Q}_{5}(x)+{Q}_{6}(x),\end{array}\end{eqnarray} \tag{ 15e }$$

$$\begin{eqnarray}\begin{array}{lll} & & {\Delta }_{\eta }={p}_{{t}_{\eta }}-{p}_{{r}_{\eta }},\end{array}\end{eqnarray} \tag{ 15f }$$

$$\begin{eqnarray}\begin{array}{lll} & & {\sigma }_{\eta }=2(mx(1-hx+g{x}^{2})+\displaystyle \frac{1}{2}(3+[hx-2]\\ & & \,\,\times {[1-hx+g{x}^{2}]}^{-1}))\\ & & \,\,\times \sqrt{[1-hx+g{x}^{2}]{x}^{-1}},\end{array}\end{eqnarray} \tag{ 15g }$$

$$\begin{eqnarray}\begin{array}{lll} & & {E}_{\eta }^{2}=mx(1-hx+g{x}^{2}),\end{array}\end{eqnarray} \tag{ 15h }$$

where for simplicity, we have set

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{0}(x)=\displaystyle \frac{1}{16}[-m{x}^{2}(\alpha -1)+4(h(\alpha -1))\\ & & \,-16\pi \beta \arctan (2gx-h)(\sqrt{4g-h})]\\ & & \,\times {[\sqrt{4g-{h}^{2}}]}^{-1}+\displaystyle \frac{1}{16}(2(1-5\alpha)\\ & & \,\times \mathrm{log}(1-hx+g{x}^{2})){(\sqrt{4g-{h}^{2}})}^{-1}+{k}_{0},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{1}(x)=[{h}^{2}x(-{m}^{2}{x}^{4}{(\alpha -1)}^{2}\\ & & \,-4m{x}^{2}(4+3(3\alpha -4)))-64(gx-2\pi \beta)]\\ & & \,\times {[128\pi (x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{2}(x)=[2h(16-24\alpha +x(mx(12+m{x}^{2}{(\alpha -1)}^{2}\\ & & \,-26\alpha +6{\alpha }^{2})))]{[128\pi (x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{3}(x)=[2hg{x}^{2}(40+{m}^{2}{x}^{4}{(\alpha -1)}^{2}-92\alpha \\ & & \,+60{\alpha }^{2}+2m{x}^{2}(9+\alpha (8\alpha -21)))]\\ & & \,\times {[128\pi (x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{4}(x)=[32hx\pi (-4+m{x}^{2}(\alpha -1)+6\alpha)\beta \\ & & \,+64g\pi {x}^{2}(1-5\alpha)\beta -256{\pi }^{2}{\beta }^{2}x]\\ & & \,\times {[128\pi (x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{5}(x)=[x(-{m}^{2}{x}^{2}{(1+g{x}^{2})}^{2}{(\alpha -1)}^{2}\\ & & \,+4g(40\alpha +g{x}^{2}(-13+5(6-5\alpha)\alpha)))]\\ & & \,\times {[128\pi (x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {Q}_{6}(x)=4mx(1+g{x}^{2})(2-4\alpha \\ & & \,+x(gx(5+\alpha (5\alpha -12))+8\pi (\alpha -1)\beta)),\end{array}\end{eqnarray}$$

and k₀ stands for a constant of integration. The total mass of the stellar object for the core is expressed in the form

$$\begin{eqnarray}\begin{array}{lll}{M}_{\eta }(r) & = & 4\pi \displaystyle \int {r}^{2}{\rho }_{\eta }dr,\,{r}^{2}=x\\ & = & [{x}^{3}(h(6+m{x}^{2})-x(m+g(10+m{x}^{2})))]\\ & & \,\,\times {[96\pi]}^{-1}.\end{array}\end{eqnarray} \tag{ 16 }$$

3.2. The Intermediate Region

We consider the middle layer to be less dense than the core. A modified Bose-Einstein condensate EoS is seen to satisfy the intermediate layer characterized by Bose-Einstein condensate matter to describe the physical features of the dense stellar object. Thus we describe this region by the equation

$$\begin{eqnarray}&&{p}_{{r}_{\iota }}=\mu {\rho }_{\iota }^{2}-a,\end{eqnarray} \tag{ 17 }$$

where μ and a are arbitrary real constants. Substituting Equation (6a) into Equation (17) we have

$$\begin{eqnarray}&&{p}_{{r}_{\iota }}=\mu {\left(\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\right)}^{2}-a.\end{eqnarray} \tag{ 18 }$$

From Equation (6b) and Equation (18) we get

$$\begin{eqnarray}\begin{array}{ll}\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}= & \displaystyle \frac{2\pi \mu }{Z}{\left(\begin{array}{c}\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right)}^{2}\\ & -\displaystyle \frac{2\pi }{Z}\left(\begin{array}{c}\displaystyle \frac{z-1}{8x\pi }+\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}+a\right).\end{array}\end{eqnarray} \tag{ 19 }$$

Putting Equation (7) and Equation (8) into Equation (19) we obtain

$$\begin{eqnarray}\begin{array}{ll}\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}= & [-16\pi (16a\pi +h(2-m{x}^{2})+x(m-2g+gm{x}^{2}))\\ & +{(h(6+m{x}^{2})-x(m+g(10+m{x}^{2})))}^{2}\mu]\\ & \times {[128\pi (1-hx+g{x}^{2})]}^{-1}.\end{array}\end{eqnarray} \tag{ 20 }$$

Integrating Equation (20) gives

$$\begin{eqnarray}\begin{array}{ll}y= & \exp \bigg[\bigg[100gx\mu +\displaystyle \frac{1}{3}m(20g+m){x}^{3}\mu -\displaystyle \frac{1}{4}h{m}^{2}{x}^{4}\mu \\ & +\displaystyle \frac{1}{5}g{m}^{2}{x}^{5}\mu -2m{x}^{2}(4\pi +3h\mu)\\ & +2(8\pi -5h\mu)\mathrm{log}(1-hx+g{x}^{2})\bigg]{\bigg[128\pi \bigg]}^{-1}\\ & +\bigg[4(-8h\pi -128a{\pi }^{2}-50g\mu +13{h}^{2}\mu)\\ & \times \arctan (\bigg[2gx-h\bigg]{\bigg[\sqrt{4g-{h}^{2}}\bigg]}^{-1})\bigg]\\ & \times {\bigg[128\pi \sqrt{4g-{h}^{2}}\bigg]}^{-1}+{k}_{2}\bigg].\end{array}\end{eqnarray} \tag{ 21 }$$

Substituting Equations (7), (8), (20) and (21) into the system (6) provides the matter variables

$$\begin{eqnarray}\begin{array}{lll} & & {e}^{2{\lambda }_{\iota }}={[1-hx+g{x}^{2}]}^{-1},\end{array}\end{eqnarray} \tag{ 22a }$$

$$\begin{eqnarray}\begin{array}{lll} & & {e}^{2{\nu }_{\iota }}={A}^{2}{y}^{2}(x)={A}^{2}{[\exp ({\tau }_{0}(x))]}^{2},\end{array}\end{eqnarray} \tag{ 22b }$$

$$\begin{eqnarray}\begin{array}{lll} & & {\rho }_{\iota }=[h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]\\ & & \,\,\times {[16\pi]}^{-1},\end{array}\end{eqnarray} \tag{ 22c }$$

$$\begin{eqnarray}\begin{array}{lll} & & {p}_{{r}_{\iota }}=\mu ([h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]\\ & & {\times {[16\pi]}^{-1})}^{2}-a,\end{array}\end{eqnarray} \tag{ 22d }$$

$$\begin{eqnarray}\begin{array}{lll} & & {p}_{{t}_{\iota }}={\tau }_{1}(x)+{\tau }_{2}(x)+{\tau }_{3}(x)+{\tau }_{4}(x)+{\tau }_{5}(x),\end{array}\end{eqnarray} \tag{ 22e }$$

$$\begin{eqnarray}\begin{array}{lll} & & {\Delta }_{\iota }={p}_{{t}_{\iota }}-{p}_{{r}_{\iota }},\end{array}\end{eqnarray} \tag{ 22f }$$

$$\begin{eqnarray}\begin{array}{lll} & & {\sigma }_{\iota }=2(mx(1-hx+g{x}^{2})+\displaystyle \frac{1}{2}(3+[hx-2]\\ & & \,\times {[1-hx+g{x}^{2}]}^{-1}))\sqrt{[1-hx+g{x}^{2}]{x}^{-1}},\end{array}\end{eqnarray} \tag{ 22g }$$

$$\begin{eqnarray}\begin{array}{lll} & & {E}_{\iota }^{2}=mx(1-hx+g{x}^{2}),\end{array}\end{eqnarray} \tag{ 22h }$$

where for simplicity we have set

$$\begin{eqnarray}&&{\tau }_{0}(x)=[100gx\mu +\displaystyle \frac{1}{3}m(20g+m){x}^{3}\mu -\displaystyle \frac{1}{4}h{m}^{2}{x}^{4}\mu \\ &&\,+\displaystyle \frac{1}{5}g{m}^{2}{x}^{5}\mu -2m{x}^{2}(4\pi +3h\mu)+2(8\pi -5h\mu)\\ &&\,\times \mathrm{log}(1-hx+g{x}^{2})]{[}^{128}\\ &&\,+[4(-8h\pi -128a{\pi }^{2}-50g\mu +13{h}^{2}\mu)\\ &&\,\times \arctan ([2gx-h]{[\sqrt{4g-{h}^{2}}]}^{-1})]\\ &&\,\times {[128\pi \sqrt{4g-{h}^{2}}]}^{-1}+{k}_{2},\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {\tau }_{1}(x)=[-4096(h-2gx)-2048mx(1-hx+g{x}^{2})]\\ & & \,\,\,\times {[32768\pi]}^{-1}+[(h(6+m{x}^{2})\\ & & \,\,\,-x(m+g(10+m{x}^{2})){)}^{2}\mu]\\ & & \,\,\,\times {[32768{\pi }^{3}(1-hx+g{x}^{2})]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {\tau }_{2}(x)=[32\pi (4-6hx+8g{x}^{2})(-16\pi (16a\pi \\ & & \,\,\,+h(2-m{x}^{2})+x(m-2g+gm{x}^{2})))]\\ & & \,\,\,\times {[32768{\pi }^{3}(1-hx+g{x}^{2})]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {\tau }_{3}(x)=[x(-16\pi (16a\pi +h(2-m{x}^{2})\\ & & \,\,\,+x(gm{x}^{2}+m-2g))+(h(6+m{x}^{2})\\ & & \,\,\,-x(m+g(10+m{x}^{2})){)}^{2}\mu {)}^{2}]\\ & & \,\,\,\times {[32768{\pi }^{3}(1-hx+g{x}^{2})]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {\tau }_{4}(x)=[x(h-2gx)(-16\pi (16a\pi +h(2-m{x}^{2})\\ & & \,+x(gm{x}^{2}+m-2g)))+(h(6+m{x}^{2}))]\\ & & \,\times {[32768{\pi }^{3}(1-hx+g{x}^{2})]}^{-1}\\ & & \,-[x{(m+g(10+m{x}^{2}))}^{2}]{[256{\pi }^{2}(1-hx+g{x}^{2})]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}&&{\tau }_{5}(x)=[2x(m-2hmx+g(10+3m{x}^{2}))\\ &&\,\,\times (-h(6+m{x}^{2})+x(m+g(10+m{x}^{2})))\mu]{[}^{32768}\\ &&\,\,+[x(-16\pi (m-2hmx+g(3m{x}^{2}-2)))]{[}^{32768}.\end{eqnarray}$$

The total mass of the stellar object in the intermediate layer is expressed as

$$\begin{eqnarray}\begin{array}{lll}{M}_{\iota }(r) & = & 4\pi \displaystyle \int {r}^{2}{\rho }_{\iota }dr,\,{r}^{2}=x\\ {M}_{\iota } & = & [{x}^{3}(h(6+m{x}^{2})-x(m+g(10+m{x}^{2})))]\\ & & \,\times {[96\pi]}^{-1}.\end{array}\end{eqnarray} \tag{ 23 }$$

3.3. The Envelope

With consideration of the density profile, the extended envelope region may have less dense material than the core and the intermediate layers. We consider a general quadratic EoS to be convenient to describe the outermost layer influenced by neutron fluid and Coulomb liquids. The quadratic EoS is taken to describe the envelope region so as to ensure the radially addressed pressure in the outer layer is less compact than the pressure in the inner layer and the intermediate layer [Mafa Takisa et al. (2019)]. This is written in the form

$$\begin{eqnarray}&&{p}_{{r}_{\epsilon }}=b{\rho }_{\epsilon }^{2}+d{\rho }_{\epsilon }-f,\end{eqnarray} \tag{ 24 }$$

where b, d and f are arbitrary real constants. From Equation (6a) and Equation (24) we have

$$\begin{eqnarray}\begin{array}{lll}{p}_{{r}_{\epsilon }} & = & b{\left(\begin{array}{c}\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right)}^{2}\\ & & +d\left(\begin{array}{c}\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right)-f.\end{array}\end{eqnarray} \tag{ 25 }$$

Equating Equations (6b) and (25) yields

$$\begin{eqnarray}\begin{array}{lll}\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx} & = & \displaystyle \frac{2\pi b}{Z}{\left(\begin{array}{c}\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right)}^{2}\\ & & +\displaystyle \frac{2\pi d}{Z}\left(\begin{array}{c}\displaystyle \frac{1-Z}{8\pi x}-\displaystyle \frac{1}{4\pi }\displaystyle \frac{dZ}{dx}-\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right)\\ & & -\displaystyle \frac{2\pi }{Z}\left(\begin{array}{c}f+\displaystyle \frac{Z-1}{8\pi x}+\displaystyle \frac{1}{16\pi }{E}^{2}\end{array}\right).\end{array}\end{eqnarray} \tag{ 26 }$$

Putting Equation (7) and Equation (8) into Equation (26), the differential equation becomes

$$\begin{eqnarray}&&\displaystyle \frac{1}{y}\displaystyle \frac{dy}{dx}=[16f\pi +(2g+10dg+m+dm)x+(1+d)gm{x}^{3}\\ &&\,-h(2+6d+m{x}^{2}+dm{x}^{2})]\times {[}^{8}\\ &&\,-[b{(}^{h}]\\ &&\,\times {[}^{128}.\end{eqnarray} \tag{ 27 }$$

Integrating Equation (27) yields

$$\begin{eqnarray}&&y=\exp [[100bgx-2m(3bh+4\pi +4\pi d){x}^{2}\\ &&\,\,+\displaystyle \frac{1}{3}bm(20g+m){x}^{3}-\displaystyle \frac{1}{4}bh{m}^{2}{x}^{4}+\displaystyle \frac{1}{5}bg{m}^{2}{x}^{5}\\ &&\,\,+2(5bh+8\pi +40\pi d)\mathrm{log}(1-hx+g{x}^{2})]{[}^{128}\\ &&\,\,+[4(13b{h}^{2}-50bg+8h\pi +8dh\pi -128f{\pi }^{2})\\ &&\,\,\times \arctan ([2gx-h]{[\sqrt{4g-{h}^{2}}]}^{-1})]\\ &&\,\,\times {[128\pi \sqrt{4g-{h}^{2}}]}^{-1}\times {[}^{128}+{k}_{3}].\end{eqnarray} \tag{ 28 }$$

Then the matter variables for the envelope layer become

$$\begin{eqnarray}&&{e}^{2{\lambda }_{\epsilon }}={[1-hx+g{x}^{2}]}^{-1},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{e}^{2{\nu }_{\epsilon }}={A}^{2}{y}^{2}(x)={A}^{2}{[\exp ({G}_{0}(x))]}^{2},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\rho }_{\epsilon }=[h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]{[16\pi]}^{-1},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}\begin{array}{ll}{p}_{{r}_{\epsilon }}= & b{([h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]{[16\pi]}^{-1})}^{2}\\ & +d([h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]\\ & \times {[16\pi]}^{-1})-f,\end{array}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}\begin{array}{ll}{p}_{{t}_{\epsilon }}= & {G}_{1}(x)+{G}_{2}(x)+{G}_{3}(x)+{G}_{4}(x)-{G}_{5}(x)\\ & +{G}_{6}(x)+{G}_{7}(x)+{G}_{8}(x)-{G}_{9}(x),\end{array}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\Delta }_{\epsilon }={p}_{{t}_{\epsilon }}-{p}_{{r}_{\epsilon }},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}\begin{array}{ll}{\sigma }_{\epsilon }= & 2(mx(1-hx+g{x}^{2})+\displaystyle \frac{1}{2}(3+[hx-2]\\ & \times [1-hx+g{x}^{2}{]}^{-1}))\sqrt{[1-hx+g{x}^{2}]{x}^{-1}},\end{array}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{E}_{\epsilon }^{2}=mx(1-hx+g{x}^{2}),\end{eqnarray} \tag{ 36 }$$

where for simplicity we have set

$$\begin{eqnarray}&&{G}_{0}(x)=[100bgx-2m(3bh+4\pi +4\pi d){x}^{2}\\ &&\,\,\,+\displaystyle \frac{1}{3}bm(20g+m){x}^{3}-\displaystyle \frac{1}{4}bh{m}^{2}{x}^{4}+\displaystyle \frac{1}{5}bg{m}^{2}{x}^{5}\\ &&\,\,\,+2(5bh+8\pi +40\pi d)\mathrm{log}(1-hx+g{x}^{2})]{[}^{128}\\ &&\,\,\,+[4(13b{h}^{2}-50bg+8h\pi +8dh\pi -128f{\pi }^{2})\\ &&\,\,\,\times \arctan ([2gx-h]{[\sqrt{4g-{h}^{2}}]}^{-1})]\\ &&\,\,\,\times {[128\pi \sqrt{4g-{h}^{2}}]}^{-1}{[}^{128}+{k}_{3},\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{1}(x)=[-4096{\pi }^{2}(h-2gx)-2048m{\pi }^{2}\\ & & \,\,\,\times x(1-hx+g{x}^{2})]{[3276{\pi }^{3}]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{2}(x)=[32\pi (4-6hx+8g{x}^{2})(16\pi (16f\pi \\ & & \,+(2g+10dg+m+dm)x+(1+d)gm{x}^{3}))]\\ & & \,\times {[32768{\pi }^{3}(x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{3}(x)=[-h(2+6d+m{x}^{2}+dm{x}^{2})\\ & & \,\,-b{(h(6+m{x}^{2})-x(10g+m+gm{x}^{2}))}^{2}]\\ & & \,\,\times {[32768{\pi }^{3}(x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{4}(x)=[x(1-hx+g{x}^{2})(16\pi (2(h(h+3dh-8f\pi)\\ & & \,\,\,-(2g+10dg+m+dm)x){)}^{2})]\\ & & \,\,\,\times {[32768{\pi }^{3}{(x(h-gx)-1)}^{2}]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{5}(x)=[x(1-hx+g{x}^{2})(1+d)hm{x}^{2}(1+d)gm{x}^{2}\\ & & \,\,+b{(h(6+m{x}^{2})-x(10g+m+gm{x}^{2}))}^{2}]\\ & & \,\,\times {[32768{\pi }^{3}{(x(h-gx)-1)}^{2}]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{6}(x)=[128\pi x(1-hx+g{x}^{2})(2gx-h)\\ & & \,\,\times (16\pi (16f\pi +(2g+10dg+m+dm)x))]\\ & & \,\,\times {[32768{\pi }^{3}{(x(h-gx)-1)}^{2}]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{7}(x)=[128\pi x(1-hx+g{x}^{2})(1+d)gm{x}^{3}\\ & & \,\,\,-h(2+6d+m{x}^{2}+dm{x}^{2})-b(h(6+m{x}^{2})\\ & & \,\,\,-x(10g+m+gm{x}^{2}){)}^{2}]\\ & & \,\,\,\times {[32768{\pi }^{3}{(x(h-gx)-1)}^{2}]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{8}(x)=[128\pi x(1-hx+g{x}^{2})(16\pi (2(5+5d)g\\ & & \,\,\,+3(1+d)gm{x}^{2}+(1+d)m(1-2hx)))]\\ & & \,\,\,\times {[32768{\pi }^{3}(x(h-gx)-1)]}^{-1},\end{array}\end{eqnarray}$$

$$\begin{eqnarray}\begin{array}{lll} & & {G}_{9}(x)=[128\pi x(1-hx+g{x}^{2})(2b(m-2hmx\\ & & \,\,+g(10+3m{x}^{2}))(h(6+m{x}^{2})\\ & & \,\,+x(m+g(10+m{x}^{2}))))]\\ & & \,\,\times {[32768{\pi }^{3}(x(h-gx)-1)]}^{-1}.\end{array}\end{eqnarray}$$

The total mass of the stellar sphere for the envelope layer is written in the form

$$\begin{eqnarray}\begin{array}{lll}{M}_{\epsilon }(r) & = & 4\pi \displaystyle \int {r}^{2}{\rho }_{\epsilon }dr,\,{r}^{2}=x,\\ {M}_{\epsilon } & = & [{x}^{3}(h(6+m{x}^{2})-x(m+g(10+m{x}^{2})))]\\ & & \times {[96\pi]}^{-1}.\end{array}\end{eqnarray} \tag{ 37 }$$

Here we match the radial pressure, gravitational potentials and electric field at the interfaces as follows:

Junction conditions at the core-intermediate boundary.

This is expressed by

$$\begin{eqnarray}\begin{array}{lll}{e}^{-2{\lambda }_{\eta }}({R}_{\eta }) & = & {e}^{-2{\lambda }_{\iota }}({R}_{\eta }),\end{array}\end{eqnarray} \tag{ 38a }$$

$$\begin{eqnarray}\begin{array}{lll}{e}^{2{\nu }_{\eta }}({R}_{\eta }) & = & {e}^{2{\nu }_{\iota }}({R}_{\eta }),\end{array}\end{eqnarray} \tag{ 38b }$$

$$\begin{eqnarray}\begin{array}{lll}{p}_{{r}_{\eta }}({R}_{\eta }) & = & {p}_{{r}_{\iota }}({R}_{\eta }),\end{array}\end{eqnarray} \tag{ 38c }$$

$$\begin{eqnarray}\begin{array}{lll}{E}_{\eta }^{2}({R}_{\eta }) & = & {E}_{\iota }^{2}({R}_{\eta }).\end{array}\end{eqnarray} \tag{ 38d }$$

Junction conditions at the intermediate-envelope interface.

This is written as

$$\begin{eqnarray}\begin{array}{lll}{e}^{-2{\lambda }_{\iota }}({R}_{\iota }) & = & {e}^{-2{\lambda }_{\epsilon }}({R}_{\iota }),\end{array}\end{eqnarray} \tag{ 39a }$$

$$\begin{eqnarray}\begin{array}{lll}{e}^{2{\nu }_{\iota }}({R}_{\iota }) & = & {e}^{2{\nu }_{\epsilon }}({R}_{\iota }),\end{array}\end{eqnarray} \tag{ 39b }$$

$$\begin{eqnarray}\begin{array}{lll}{p}_{{r}_{\iota }}({R}_{\iota }) & = & {p}_{{r}_{\epsilon }}({R}_{\iota }),\end{array}\end{eqnarray} \tag{ 39c }$$

$$\begin{eqnarray}\begin{array}{lll}{E}_{\iota }^{2}({R}_{\iota }) & = & {E}_{\epsilon }^{2}({R}_{\iota }).\end{array}\end{eqnarray} \tag{ 39d }$$

Junction conditions at the envelope-surface boundary. The interior and exterior line elements (1) and (2) should match smoothly at the surface r = R_ε . Therefore

$$\begin{eqnarray}\begin{array}{lll}{e}^{2{\lambda }_{\epsilon }}({R}_{\epsilon }) & = & {\left(\begin{array}{c}1-\displaystyle \frac{2M}{{R}_{\epsilon }}+\displaystyle \frac{{Q}^{2}}{{R}_{\epsilon }^{2}}\end{array}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 40a }$$

$$\begin{eqnarray}\begin{array}{lll}{e}^{2{\nu }_{\epsilon }}({R}_{\epsilon }) & = & \left(\begin{array}{c}1-\displaystyle \frac{2M}{{R}_{\epsilon }}+\displaystyle \frac{{Q}^{2}}{{R}_{\epsilon }^{2}}\end{array}\right),\end{array}\end{eqnarray} \tag{ 40b }$$

$$\begin{eqnarray}\begin{array}{lll}{p}_{{r}_{\epsilon }}({R}_{\epsilon }) & = & 0,\end{array}\end{eqnarray} \tag{ 40c }$$

$$\begin{eqnarray}\begin{array}{lll}{E}_{\epsilon }^{2}({R}_{\epsilon }) & = & \displaystyle \frac{{Q}^{2}}{{R}_{\epsilon }^{4}}.\end{array}\end{eqnarray} \tag{ 40d }$$

They yield

$$\begin{eqnarray}\begin{array}{lll} & & (1-h{R}_{\epsilon }^{2}+g{R}_{\epsilon }^{2n})=\bigg(1+\displaystyle \frac{2{P}_{o}(x)}{{R}_{\epsilon }}\\ & & \,\,\,\,{+\displaystyle \frac{m(1-h{R}_{\epsilon }+g{R}_{\epsilon }^{2})}{{R}_{\epsilon }^{3}}\bigg)}^{-1},\end{array}\end{eqnarray} \tag{ 41a }$$

$$\begin{eqnarray}\begin{array}{lll} & & {A}^{2}{[\exp ({G}_{0}(x))]}^{2}=\bigg(1+\displaystyle \frac{2{P}_{o}(x)}{{R}_{\epsilon }}\\ & & \,\,\,\,\,+\displaystyle \frac{m(1-h{R}_{\epsilon }+g{R}_{\epsilon }^{2})}{{R}_{\epsilon }^{3}}\bigg),\end{array}\end{eqnarray} \tag{ 41b }$$

$$\begin{eqnarray}\begin{array}{lll} & & \,0=b{({P}_{1}(x))}^{2}+d([h(6+m{x}^{2}-x(m\\ & & \,\,\,+g(10+m{x}^{2})))]{[16\pi]}^{-1})-f,\end{array}\end{eqnarray} \tag{ 41c }$$

$$\begin{eqnarray}\begin{array}{lll} & & \,{E}^{2}=\displaystyle \frac{{Q}^{2}}{{R}_{\epsilon }^{4}}=m{R}_{\epsilon }(1-h{R}_{\epsilon }+g{R}_{\epsilon }^{2}),\end{array}\end{eqnarray} \tag{ 41d }$$

where for simplicity we have set

$$\begin{eqnarray}\begin{array}{ll}{P}_{o}(x)= & [{x}^{3}(h(6+m{x}^{2})-x(m+g(10+m{x}^{2})))]\\ & \times {[96\pi]}^{-1}+{k}_{1},\\ {P}_{1}(x)= & ([h(6+m{x}^{2}-x(m+g(10+m{x}^{2})))]\\ & \times {[16\pi]}^{-1}).\end{array}\end{eqnarray}$$

We observe from system (41) that there are a sufficient number of free parameters in four equations. This indicates that the matching conditions are easily satisfied. In system (41) above, it is complicated to solve the value of R_ε analytically. However the value of R_ε can be obtained geometrically from the plots, especially at the vanishing point of radial pressure which is clearly shown to be R_ε = 9.96 km. The same choice of electric field E² is applied for all three layers. From systems (31) and (32) we observe that at the core-intermediate junction at R_η and at the intermediate-envelope junction at R_ι , the electric field exhibits continuity throughout the interior of the stellar sphere.

For a well behaved physical model, the matter variables and other physical quantities have to satisfy the following conditions.

• (i)  
  The energy density, radial and tangential pressures should be regular, finite and continuous, and greater than or equal to zero.
• (ii)  
  The metric potentials e^2λ and e^2ν should be greater than or equal to zero.
• (iii)  
  The radial sound speed should be less than the speed of light so as to obey the causality condition. These conditions in each layer are written as

  $$\begin{eqnarray}\begin{array}{lll}{\nu }_{\eta } & = & A,\end{array}\end{eqnarray} \tag{ 42a }$$

  $$\begin{eqnarray}\begin{array}{lll}{\nu }_{\iota } & = & [(h(6+m{x}^{2})-x(m+g(10+m{x}^{2})))\mu]\\ & & \times {[8\pi]}^{-1},\end{array}\end{eqnarray} \tag{ 42b }$$

  $$\begin{eqnarray}\begin{array}{lll}{\nu }_{\epsilon } & = & [(8\pi d+bh(6+m{x}^{2})-bx(m\\ & & +g(10+m{x}^{2})))]{[8\pi]}^{-1}.\end{array}\end{eqnarray} \tag{ 42c }$$
• (iv)  
  The model should satisfy the energy conditions: strong energy condition (SEC), weak energy condition (WEC) or null energy condition (NEC). i.e. SEC: ρ – p_r – 2p_t ≥ 0, WEC: ρ – 3p_t ≥ 0, ρ – 3p_r ≥ 0, NEC: ρ – p_r ≥ 0, ρ – p_t ≥ 0.
• (v)  
  In general relativity, stability of the stellar sphere should satisfy the condition $\Gamma =\displaystyle \frac{\rho +{p}_{r}}{{p}_{r}}\displaystyle \frac{d{p}_{r}}{d\rho }\ge \displaystyle \frac{4}{3}$. This condition for each region becomes

  $$\begin{eqnarray}\begin{array}{ll}{\Gamma }_{\eta }= & [(\alpha (-h(6+m{x}^{2})(1+\alpha)+x(m+g(10+m{x}^{2}))\\ & \times (1+\alpha)+16\pi \beta))][-h((6+m{x}^{2})\alpha \\ & +x(m+g(10+m{x}^{2}))\alpha +16\pi \beta){]}^{-1},\end{array}\end{eqnarray} \tag{ 43 }$$

  $$\begin{eqnarray}\begin{array}{ll}{\Gamma }_{\iota }= & [{B}_{1}(x)(-a+({B}_{2}(x)(16\pi +h(6+m{x}^{2})\mu \\ & -x(m+g(10+m{x}^{2}))))+8\pi (-a+((h(6+m{x}^{2})\\ & -x(10g+m+gm{x}^{2}){)}^{2}\mu)))]{[256{\pi }^{2}]}^{-1},\end{array}\end{eqnarray} \tag{ 44 }$$

  $$\begin{eqnarray}&&{\Gamma }_{\epsilon }=[{B}_{3}(x)(b{(}^{h}\\ &&\,+16\pi (-{B}_{5}(x)))][(2048{\pi }^{2}(-\pi f+\displaystyle \frac{1}{16}d(h(6+m{x}^{2})\\ &&\,-x(m+g(10+m{x}^{2})))+{B}_{4}(x))){]}^{-1},\end{eqnarray} \tag{ 45 }$$

  where

  $$\begin{eqnarray}&&{B}_{1}(x)=((h(6+m{x}^{2})-x(m+g(10+m{x}^{2}))))\mu,\\ &&{B}_{2}(x)=(h(6+m{x}^{2})-x(10g+m+gm{x}^{2})),\\ &&{B}_{3}(x)=((8\pi d+bh(6+m{x}^{2})\\ &&\,\,-bx(m+g(10+m{x}^{2}))),\\ &&{B}_{4}(x)=[(b(h(6+m{x}^{2})-x(m\\ &&\,\,+g(10+m{x}^{2})){)}^{2})]{[}^{256},\\ &&{B}_{5}(x)=16\pi f+(1+d)h(6+m{x}^{2})\\ &&\,\,-(1+d)x(10g+m+gm{x}^{2})),\\ &&{B}_{6}(x)=(32\pi (4-6hx+8g{x}^{2})(-16\pi (16\pi a\\ &&\,+h(2-m{x}^{2})+x(-2g+m+gm{x}^{2}))\\ &&\,+{(}^{h}\mu)).\end{eqnarray}$$
• (vi)  
  The surface redshift for the relativistic stellar objects has been discussed by many authors in the past (Baraco & Hamity (2002)). Bohmer & Harko (2006) showed that for an imperfect fluid the value of the redshift should not exceed 5. The value of the surface redshift was later extended in the treatment by Ivanov (2002) who suggested the maximum value be z_s ≤ 5.211 for an anisotropic relativistic stellar object. The surface redshift is calculated using the form

  $$\begin{eqnarray}&&{z}_{s}=\displaystyle \frac{1}{\sqrt{1-\displaystyle \frac{2M(r)}{r}}}-1,\end{eqnarray} \tag{ 46 }$$

  where M stands for the total mass of the stellar object. In our study we obtain

  $$\begin{eqnarray}\begin{array}{lll}{z}_{s} & = & \displaystyle \frac{1}{\sqrt{1+\displaystyle \frac{{x}^{5/2}(-h(6+m{x}^{2})+x(m+g(10+m{x}^{2})))}{48\pi }}}\\ & & -1.\end{array}\end{eqnarray} \tag{ 47 }$$
• (vii)  
  A thorough analysis of the mass radius ratio for a relativistic stellar object is described in Jasim et al. (2021). It has been noted that the compactification factor for the relativistic object is categorized into different classes: (i) normal star: $\displaystyle \frac{M(r)}{r}\sim {10}^{-5}$, (ii) white dwarf: $\displaystyle \frac{M(r)}{r}\sim {10}^{-3}$, (iii) neutron star: ${10}^{-1}\lt \displaystyle \frac{M(r)}{r}\lt \displaystyle \frac{1}{4}$, (iv) ultra-dense compact star: $\displaystyle \frac{1}{4}\lt \displaystyle \frac{M(r)}{r}\lt \displaystyle \frac{1}{2}$ and (v) black hole: $\displaystyle \frac{M(r)}{r}=\displaystyle \frac{1}{2}$. The compactification factor is obtained in the form

  $$\begin{eqnarray}&&\mu (r)=\displaystyle \frac{2M(r)}{r}.\end{eqnarray} \tag{ 48 }$$

  In our treatment we get

  $$\begin{eqnarray}&&\mu (x)=\displaystyle \frac{{x}^{5/2}(-h(6+m{x}^{2})+x(m+g(10+m{x}^{2})))}{48\pi }.\end{eqnarray} \tag{ 49 }$$

We discuss the physical features of the matter variables and other quantities which were generated by using the Python programming language. The plots generated include energy density (Fig. 1), radial pressure (Fig. 2), tangential pressure (Fig. 3), gravitational potentials (Figs. 4 and 5), mass (Fig. 6), energy conditions (Figs. 7, 8, 9 and 10), the measure of anisotropy (Fig. 11), adiabatic index (Fig. 12), radial sound speed (Fig. 13), charge density (Fig. 14), electric field (Fig. 15), mass radius ratio (Fig. 16) and surface redshift (Fig. 17). We generate graphs utilizing the following specified values of the constants: a = ± 1.187 × 10⁻⁵, b = 319.2, d = 0.7283, f = ± 1.547 × 10⁻⁵, g = ± 1.5 × 10⁻⁵, h = ± 0.004587, k_o = 0.5, k₁ = 0.0342, k₂ = 0.2, k₃ = 0.2, m = ± 7.15 × 10⁻⁵, A = 0.1, α = 0.19897, β = 4.9166 × 10⁻⁵ and μ = ± 481.124. All graphs are plotted versus radial distance in the specified domain of radius as per Pant et al. (2020).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our model, Figures 1 and 2 affirm that the energy density and radial pressure are continuous functions, have their maximum at the center, and monotonically decrease towards the surface. We also observe that these trends are similar to the findings in Pant et al. (2020), Bisht et al. (2021), Maharaj et al. (2014), Pant et al. (2019), Gedela et al. (2019) and Sunzu & Danford (2017). It is clearly shown that the radial pressure vanishes at the surface (R_ε = 9.96 km). Figure 3 indicates that tangential pressure is an increasing function towards the surface. This physical behavior is also found in the models by Thirukkanesk & Ragel (2014), Maharaj & Mafa Takisa (2013) and Ngubelanga & Maharaj (2015). Figures 4 and 5 demonstrate that metric potentials are continuously increasing functions. Figure 6 indicates that the mass is a monotonically increasing function with the increase of radial distance. The energy conditions in Figures 7, 8, 9 and 10 are decreasing functions that are positive. These physical features are also found in the works by Pant et al. (2020), Maurya et al. (2018) and Maurya & Govender (2017). The measure of anisotropy is an increasing function with a discontinuous point between the intermediate layer and the envelope layer (Fig. 11). This demonstrates that tangential pressure is greater than radial pressure. A similar profile is also obtained by Pant et al. (2020), Thirukkanesk & Ragel (2014), Bhar et al. (2017), Ngubelanga et al. (2015) and Murad (2016). In Figure 12, the stability condition is found to satisfy the criteria within the stellar interior $(\Gamma \ge \displaystyle \frac{4}{3})$. This shows a gradual increase in Γ towards the surface. In our treatment, we obtained the minimum adiabatic index Γ = 1.9769, affirming that the object is stable from gravitational collapse. We can also see a similar physical feature in the treatments by Pant et al. (2020), Bisht et al. (2021), Jasim et al. (2018) and Bhar et al. (2017). Figure 13 shows that the speed of sound is less than the speed of light. It is found to be in the range (0.25523 ≤ ν ≤ 0.69613). Similar physical features are also found by Pant et al. (2020), Pant et al. (2019), Gedela et al. (2019), Maharaj & Mafa Takisa (2013), Bhar et al. (2017), Murad (2016) and Maharaj & Mafa Takisa (2012). We observe in Figures 14 and 15 that charge density and electric field are monotonically increasing functions. From Figure 16, we observe that the trend of mass radius ratio is an increasing function, and reaches its maximum value μ(x) = 0.3933 which satisfies the condition for an ultra-dense compact star. The same physical profile is also observed in the treatments of Pant et al. (2020), Bisht et al. (2021) and Jasim et al. (2018). In Figure 17 we observe that the surface redshift is continuous and monotonically increasing in nature with radial coordinate (r), attaining a maximum value at z_s = 1.05188 which is physical and in an acceptable range. The physical analysis demonstrates that the presence of charge continues to yield a physically reasonable model of a three-layer stellar object.

In this paper we formulated a new class of exact solutions for superdense stellar spheres comprising three interior layers. In our model we used the Einstein-Maxwell field equations incorporating different EoSs. The core region obeys a linear EoS, the intermediate layer has a quadratic EoS (with no linear term) and the envelope layer is outfitted with a general quadratic EoS. In our model we made a specific choice for the electric field and one for the gravitational potentials. Our model contains the electric field which is absent in Pant et al. (2020) and Bisht et al. (2021). Setting electric field to zero (E = 0), we regain the anisotropic uncharged model developed by Pant et al. (2020), and other earlier models in the core-envelope setting. The plots generated indicate that the gravitational potentials, matter variables and other physical quantities are well behaved and compatible with astrophysical studies. In this case, they are regular, continuous and free from singularity throughout the interior of the stellar objects. The physical analysis demonstrates that our model is well behaved. The results obtained in this paper are significant and allow us to describe the physical structures, features and properties of charged anisotropic superdense stellar objects with three layers. In future work, other results can be generated by considering different forms of gravitational potential, EoSs and electric field.

SDM thanks the University of KwaZulu-Natal for financial support during preparation of this work. JMS expresses his thanks to the University of Dodoma-Tanzania for creating a favorable environment for research works. ASL and EWM express their sincere thanks to the University of Dar es Salaam and the Ministry of Education, Science and Technology-Tanzania (MOEST) for continuous support.