Quark Stars in Massive Brans–Dicke Gravity with Tolman–Kuchowicz Spacetime

Abstract

In this paper, we construct anisotropic model representing salient features of strange stars in the framework of massive Brans–Dicke gravity. We formulate the field equations for Tolman–Kuchowicz ansatz by incorporating the MIT bag model. Junction conditions are applied on the boundary of the stellar model to evaluate the unknown constants in terms of mass and radius of the star. The radius of the strange star candidate PSR J1614-2230 is predicted by assuming maximum anisotropy at the surface of the star for different values of the coupling parameter, mass of the scalar field and bag constant. We examine various properties as well as the viability and stability of the anisotropic sphere. We conclude that the astrophysical model agrees with the essential criteria of a physically realistic model for higher values of the coupling parameter as well as mass of the scalar field.

1. Introduction

Self-gravitating systems such as stellar and galactic structures constitute a major part of the visible universe. The study of the evolution of these systems plays a crucial role in revealing various hidden aspects (composition, evolution, age) of the cosmos. One of the important evolutionary phases of a star is the process of gravitational collapse near its death which leads to the formation of dense compact objects. The neutron star is one of the intriguing outcomes of collapse with a dense core of about $1M_{⨀}$ to $3M_{⨀}$ ($M_{⨀}$ is the solar mass). The degeneracy pressure of neutrons resists the strong inward gravitational pull and halts the process of collapse. Since radiation emitted by neutron stars is less as compared to other cosmic objects, therefore, they are difficult to detect. Baade and Zwicky [1] predicted the existence of neutron stars in 1934. However, the first rapidly rotating pulsar was discovered in 1967 which pulsated at a regular interval of 1.37 s [2].

According to the Tolman–Oppenheimer–Volkoff limit, a stable neutron star must have a mass of at most 3$M_{⨀}$. If the mass exceeds this limit, it is hypothesized that neutrons lose their individuality under extreme pressure and breakdown into quarks. Quark matter can exist in the universe in two scenarios: the transition phase of quark-hadron in the initial era of the cosmos or the evolution of neutron stars into strange quark stars [3]. A quark star is smaller in size but ultra-dense as compared to the neutron star. However, increased pressure in its core stops quark stars from collapsing into black holes. It has been proposed that the amount of energy emitted in extremely luminous supernovae (observed once in every thousand supernovae) can be explained with the help of stable strange quark stars [4,5]. Moreover, estimates of radii of some stellar objects (LMC X-4, 4U 1820-30, Her X-1, etc.) suggest that their structure and characteristics may be similar to that of strange quark stars.

Interactions among tightly packed nuclear matter generate anisotropy in a system with high density [6]. Bowers and Liang proposed that extreme density at the core of celestial structures (like compact stars) makes them highly anisotropic in nature [7]. Anisotropy may be the result of different physical processes such as the presence of superfluid, phase transition or magnetic field [8,9,10]. Researchers have investigated salient properties of cosmic objects by introducing anisotropic pressure in matter source. Herrera and Santos [11] studied different aspects of anisotropic self-gravitating objects and identified the causes of local anisotropy. Harko and Mak [12] considered the anisotropy factor to obtain static interior solutions of relativistic objects. Hossein et al. [13] included cosmological constant to explore the stability of anisotropic systems. Viable solutions were derived for spherical systems in equilibrium by Paul and Deb [14]. Recently, Maurya et al. [15] presented a class of physically realistic solutions representing static anisotropic stars by incorporating Buchdahl’s ansatz.

Equation of state (EoS) relating different physical properties of the matter distribution is often introduced in the modeling of an astrophysical configuration to describe the realistic scenario within the stellar structure. An agreement on the most suitable EoS representing quark stars has not been achieved. Witten [3] conjectured that the quark star is composed of strange quark matter (SQM) containing strange, up and down quarks. The confined hadrons can exist in their true ground state in a hypothetical quark star [16,17]. As the EoS for neutron star does not accurately describe the compact objects like SAX J 1808.4-3658, Her X-1, PSR J1614-2230, etc. therefore, the MIT bag model (SQM EoS) [18] is chosen as the most suitable representative of quark stars. Moreover, recent statistics obtained from the collision of binary neutron stars (GW170817 [19] and GW190425 [20]) have provided an estimation of masses of neutron stars as well as strange stars which support the assumption of MIT bag model as an EoS for quark stars.

The difference between the energy density of true and false vacuums is incorporated in the MIT bag EoS through a bag constant ($\mathcal{B}$). The bag constant lies within the range 58.9–91.5 MeV/fm${}^3$ [21] for massless strange quarks whereas the range of $\mathcal{B}$ for massive quarks is of 56–78 MeV/fm${}^3$ [22]. However, researchers have also considered larger values of $\mathcal{B}$. Xu et al. [23] found that the cosmic body LMXB EXO 0748-676 can serve as a strange star candidate for two values of $\mathcal{B}$. Further, the experimental data of RHIC and CERN-SPS also allow higher values of the bag constant [24,25]. The MIT bag model has been employed in literature to discuss the structure and characteristics of quark stars [26,27,28]. Rahaman et al. [29] calculated the mass of a star of radius $9.9 \mathrm{km}$ using MIT bag model. Bhar [30] employed Krori and Barura ansatz to study a hybrid star made of quark matter. Anisotropic quark star models have been constructed in the presence of charge as well [31,32,33]. Deb et al. [34] examined different physical properties of strange stars for $\mathcal{B}=83$ MeV/fm${}^3$, 100 MeV/fm${}^3$, 120 MeV/fm${}^3$. Bhar [35] investigated the viability as well as the stability of quark models through the condition of embedding class-one.

Dirac’s hypothesis in 1937 suggested that the gravitational constant (G) varies with respect to cosmic time [36,37]. In 1961, Brans and Dicke [38] incorporated the dynamical gravitational constant in terms of a massless scalar field ($\varphi$) and modified general relativity (GR). Brans–Dicke (BD) theory, a scalar-tensor theory, is based on the Machian principle. The matter distribution is coupled to the dynamical scalar field through a tunable coupling parameter (${\omega }_{BD}$) in Jordan frame. The weak-field tests put a lower limit on the coupling parameter as ${\omega }_{BD}>40,000$ [39]. However, as higher values of scalar field support rapid expansion of the cosmos, the inflationary era is explained for small values of the coupling parameter [40]. The issue is resolved through a potential function that allows a massive scalar field ($\Psi$). Thus, the BD gravity is modified to Massive BD (MBD) theory in which the mass of scalar field ($m_{\Psi }$) provides a finite range of the scalar field in terms of its compton wavelength (${\lambda }_{\Psi }$). The constraints on the coupling parameter due to solar system observations vanish for $m_{\Psi }\gtrsim 2\times {10}^{-25}$ GeV and ${\omega }_{BD}>-\frac{3}{2}$ [41].

Various astrophysical phenomena have been discussed in scalar-tensor gravity. Sotani [42] discussed deviations from GR while studying neutron stars in the presence of a massless scalar field. Silva et al. [43] analyzed the inertia of anisotropic neutron stars. The effects of a massive scalar field on rapidly rotating neutron stars were studied by Doneva and Yazadjiev [44] and it was verified that deviations from GR increase for larger inertia. This work was extended by employing a self-interacting potential to examine the structure of slowly rotating neutron stars [45]. Researchers have constructed different models for quark stars in modified theories of gravity as well [46,47,48,49]. Recently, we have used the condition of embedding class-one to generate anisotropic model for quark stars in the background of MBD gravity [50].

In this paper, we formulate anisotropic model obeying MIT bag model with Tolman–Kuchowicz (TK) ansatz in the context of MBD theory. We further check the viability and stability of the resulting stellar model. The matter variables obeying the MIT bag model are obtained from MBD field equations in Section 2. The necessary matching conditions on the boundary of the model are applied in Section 3. In Section 4, we examine various physical features, viability and stability of the anisotropic solution. The results are summarized in the last section.

2. Massive Brans–Dicke Theory and Matter Variables

The action of MBD gravity in Jordan frame with $G_0=1$ is defined as

$$S=\int \sqrt{-g}(\mathcal{R}\Psi -\frac{{\omega }_{BD}}{\Psi }{\nabla }^{\gamma }{\nabla }_{\gamma }\Psi -V(\Psi )+L_m)d^{4}x,$$

(1)

where $g=|g_{\gamma \delta }|$, $\mathcal{R}$ is the Ricci scalar and $L_m$ represents the matter lagrangian. The MBD gravity can also be discussed in Einstein frame of reference. The conformal transformations ${\widehat{g}}_{\gamma \delta }={\mathcal{A}}^{-2}(\Psi )g_{\gamma \delta }$ and $\Psi ={\mathcal{A}}^{-2}(\widehat{\Psi })$ yield metric (${\widehat{g}}_{\gamma \delta }$) and scalar field ($\widehat{\Psi }$) for Einstein frame. The field and evolution equations are respectively, given as

$$\begin{array}{ccc}\hfill G_{\gamma \delta }& =& \frac{1}{\Psi }[T_{\gamma \delta }^{(m)}+T_{\gamma \delta }^{\Psi }]=\frac{1}{\Psi }[T_{\gamma \delta }^{(m)}+{\Psi }_{,\gamma ;\delta }-g_{\gamma \delta }\square \Psi +\frac{{\omega }_{BD}}{\Psi }({\Psi }_{,\gamma }{\Psi }_{,\delta }\hfill \\ \hfill & -& \frac{g_{\gamma \delta }{\Psi }_{,\mu }{\Psi }^{,\mu }}{2})-\frac{V(\Psi )g_{\mu \nu }}{2}],\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \square \Psi & =& \frac{T^{(m)}}{3+2{\omega }_{BD}}+\frac{1}{3+2{\omega }_{BD}}(\Psi \frac{dV(\Psi )}{d\Psi }-2V(\Psi )),\hfill \end{array}$$

(3)

where the matter source is described by the energy-momentum tensor $T_{\gamma \delta }^{(m)}$ and $T^{(m)}=T_{\gamma }^{\gamma (m)}$. Moreover, □ indicates the d’Alembertian operator.

We consider a static sphere represented by the following line element

$$d{s}^2=e^{\nu (r)}d{t}^2-e^{\lambda (r)}d{r}^2-r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2).$$

(4)

The relativistic motion of particles in compact objects leads to a random distribution throughout the interior region. Consequently, compact stellar structures are characterized by variable energy density and anisotropic pressure. Therefore, we employ the energy-momentum tensor specifying the anisotropic interior in terms of density ($\rho$), radial ($p_r$) and tangential ($p_{\perp }$) pressures as follows

$$T_{\gamma \delta }^{(m)}=(\rho +p_{\perp })U_{\gamma }{U}_{\delta }-p_{\perp }{g}_{\gamma \delta }+(p_r-p_{\perp })S_{\gamma }{S}_{\delta },$$

(5)

where $U_{\gamma }=(e^{\frac{\nu }{2}},0,0,0)$ is the 4-velocity and $S_{\gamma }=(0,-e^{\frac{\lambda }{2}},0,0)$. Using Equations (2)–(5), the field equations are obtained as

$$\begin{array}{ccc}\hfill \frac{1}{r^2}-\left(\frac{1}{r^2}-\frac{{\lambda }^{\prime }}{r}\right)e^{-\lambda }& =& \frac{1}{\Psi }(\rho +e^{-\lambda }({\Psi }^{″}+(\frac{2}{r}-\frac{{\lambda }^{\prime }}{2}){\Psi }^{\prime }+\frac{{\omega }_{BD}}{2\Psi }\hfill \\ & \times & {\Psi }^{\prime 2}-e^{\lambda }\frac{V(\Psi )}{2})),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill -\frac{1}{r^2}+e^{-\lambda }\left(\frac{1}{r^2}+\frac{{\nu }^{\prime }}{r}\right)& =& \frac{1}{\Psi }(p_r-e^{-\lambda }(\left(\frac{2}{r}+\frac{{\nu }^{\prime }}{2}\right){\Psi }^{\prime }-\frac{{\omega }_{BD}}{2\Psi }{\Psi }^{\prime 2}\hfill \\ & -& e^{\lambda }\frac{V(\Psi )}{2})),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \frac{e^{-\lambda }}{4}\left(2{\nu }^{″}+{\nu }^{\prime 2}-{\lambda }^{\prime }{\nu }^{\prime }+2\frac{{\nu }^{\prime }-{\lambda }^{\prime }}{r}\right)& =& \frac{1}{\Psi }(p_{\perp }-e^{-\lambda }({\Psi }^{″}+\left(\frac{1}{r}-\frac{{\lambda }^{\prime }}{2}+\frac{{\nu }^{\prime }}{2}\right){\Psi }^{\prime }\hfill \\ & +& \frac{{\omega }_{BD}}{2\Psi }{\Psi }^{\prime 2}-e^{\lambda }\frac{V(\Psi )}{2})),\hfill \end{array}$$

(8)

where ${}^{\prime }$ denotes differentiation with respect to r. The evolution Equation (3) becomes

$$\begin{array}{ccc}\hfill \square \Psi & =& -e^{-\lambda }\left[\left(\frac{2}{r}-\frac{{\lambda }^{\prime }}{2}+\frac{{\nu }^{\prime }}{2}\right){\Psi }^{\prime }(r)+{\Psi }^{″}(r)\right],\hfill \\ \hfill & =& \frac{1}{3+2{\omega }_{BD}}\left[g_{\gamma \delta }{T}^{\gamma \delta }+\left(\Psi \frac{dV(\Psi )}{d\Psi }-2V(\Psi )\right)\right].\hfill \end{array}$$

(9)

Extreme temperature and pressure at the core of a neutron star may transform it into a quark star composed of up strange $(s)$, up $(u)$ and down $(d)$ flavors. We assume that the MIT bag model governs the state determinants (density and pressure) of these relativistic stars. It is further assumed that non-interacting and massless quarks occupy the interior of the configuration. The MIT bag model defines the individual quark pressure ($p^f$) as

$$p_r=\sum _f{p}^f-\mathcal{B},f=u,d,s.$$

(10)

The total external bag pressure (or bag constant) balances the individual pressure. The deconfined quarks within the bag model have the following energy density

$$\rho =\sum _f{\rho }^f+\mathcal{B},$$

(11)

where the density of a flavor (${\rho }^f$) is related to the corresponding pressure as $\frac{{\rho }^f}{3}=p^f$. Thus, Equations (10) and (11) are combined to formulate the following EoS of MIT bag model

$$p_r=\frac{1}{3}(\rho -4\mathcal{B}).$$

(12)

The metric potentials $\nu (r)=B{r}^2+2lnF$ and $\lambda (r)=ln(1+a{r}^2+b{r}^4)$ ($a,b,B$ and F are constants) define the Tolman–Kuchowicz (TK) spacetime [51,52]. Salient features of compact objects have been investigated by employing the TK metric in GR as well as modified theories of gravity [53,54,55]. The field Equations (6)–(8) are rewritten in terms of TK metric functions as

$$\begin{array}{ccc}\hfill \rho & =& \frac{1}{2r\Psi {\left(a{r}^2+b{r}^4+1\right)}^2}\left(2{\Psi }^2\left(a^2{r}^3+ar\left(2b{r}^4+3\right)+b{r}^3\left(b{r}^4+5\right)\right)\right.\hfill \\ & +& r\Psi V(\Psi ){\left(a{r}^2+b{r}^4+1\right)}^2-r{\omega }_{BD}\left(a{r}^2+\left.b{r}^4+1\right){\Psi }^{\prime 2}-2\Psi \left(r\left(a{r}^2\right.\right.\right.\hfill \\ & +& \left.\left.\left.b{r}^4+1\right){\Psi }^{″}(r)+\left(a{r}^2+2\right){\Psi }^{\prime }\right)\right),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill p_r& =& \frac{1}{4\Psi {\left(a{r}^2+b{r}^4+1\right)}^2}\left(\Psi \left(r{\Psi }^{\prime }\left(aB{r}^2+a+bB{r}^4+2b{r}^2+B\right)-\left(a{r}^2\right.\right.\right.\hfill \\ & +& \left.\left.b{r}^4+1\right){\Psi }^{″}\right)+2{\Psi }^2\left(aB{r}^2+a+bB{r}^4+2b{r}^2+B\right)-{\omega }_{BD}\left(a{r}^2\right.\hfill \\ & +& \left.\left.b{r}^4+1\right){\Psi }^{\prime 2}\right)-\mathcal{B},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill p_{\perp }& =& \frac{1}{4r\Psi {\left(a{r}^2+b{r}^4+1\right)}^2}\left(2r{\Psi }^2\left(2a^2{r}^2+a\left(4b{r}^4+2B^2{r}^4-B{r}^2+1\right)\right.\right.\hfill \\ & +& \left.2b^2{r}^6+2B^2\left(b{r}^6+r^2\right)-3bB{r}^4+B\right)+\Psi \left({\Psi }^{\prime }\left(a{r}^2\left(B{r}^2-7\right)\right.\right.\hfill \\ & +& \left.bB{r}^6-10b{r}^4+B{r}^2-4\right)-r\left(a{r}^2+b{r}^4+1\right)\left(4\mathcal{B}\left(a{r}^2+b{r}^4+1\right)\right.\hfill \\ & -& \left.\left.\left.3{\Psi }^{″}\right)\right)+3r{\omega }_{BD}\left(a{r}^2+b{r}^4+1\right){\Psi }^{\prime 2}\right).\hfill \end{array}$$

(15)

In this study, we choose $V(\Psi )=\frac{1}{2}{m}_{\Psi }^2{\Psi }^2$, where $m_{\Psi }$ denotes the mass of the scalar field. This form of potential function has been used to study neutron stars [44,56].

Junction Conditions

The constants ($a,b,B,F$) completely specify the solution in Equations (13)–(15) In order to evaluate these constants, we apply constraints at the hypersurface ($\mathsf{\Sigma }:r=R$) which ensure a smooth junction between interior and exterior geometries. The Schwarzschild spacetime describes the external vacuum as

$$d{s}^2=(1-\frac{2M}{r})d{t}^2-\frac{1}{(1-\frac{2M}{r})}d{r}^2-r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2),$$

(16)

where M is the mass. Using the procedure in [57], the scalar field associated with the external Schwarzschild metric is obtained as $\Psi =e^{(1-\frac{2M}{r})}$. The matching at the boundary surface ($h=r-R=0$) is continuous when the following conditions are satisfied

$$\begin{array}{cc}\hfill & {(d{s}_{-}^2)}_{\mathsf{\Sigma }}={(d{s}_{+}^2)}_{\mathsf{\Sigma }},{(K_{i{j}_{-}})}_{\mathsf{\Sigma }}={(K_{i{j}_{+}})}_{\mathsf{\Sigma }},\end{array}$$

(17)

$$\begin{array}{c}{(\Psi {(r)}_{-})}_{\mathsf{\Sigma }}={(\Psi {(r)}_{+})}_{\mathsf{\Sigma }},{({\Psi }^{\prime }{(r)}_{-})}_{\mathsf{\Sigma }}={({\Psi }^{\prime }{(r)}_{+})}_{\mathsf{\Sigma }}.\end{array}$$

(18)

Here $K_{ij}$ is curvature whereas interior and exterior spacetimes are represented by subscripts − and +, respectively. The metric defining the hypersurface is written as

$$d{s}^2=d{T}^2-R^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2),$$

(19)

where T represents the proper time on $\mathsf{\Sigma }$. The curvature is given by

$$K_{ij}^{\pm }=-\frac{{\partial }^2{z}_{\pm }^{\gamma }}{\partial {\eta }^i{\eta }^j}{n}_{\gamma }^{\pm }-{\Gamma }_{\delta \mu }^{\gamma }\frac{\partial z_{\pm }^{\delta }}{\partial {\eta }^i}\frac{\partial z_{\pm }^{\mu }}{\partial {\eta }^j}{n}_{\gamma }^{\pm },$$

where ${\eta }^i$ are the coordinates defined on $\mathsf{\Sigma }$. The normal ($n_{\gamma }^{\pm }$) to the boundary is defined as

$$n_{\gamma }^{\pm }=\pm \frac{dh}{d{z}^{\gamma }}{\left|g^{\delta \mu }\frac{dh}{d{z}^{\delta }}\frac{dh}{d{z}^{\mu }}\right|}^{\frac{-1}{2}},$$

with $n_{\gamma }{n}^{\gamma }=1$. Here $z_{\pm }^{\gamma }$ are the coordinates of internal/external regions. The unit normal vectors turn out to be

$$n_{\gamma }^{-}=(0,e^{\frac{\lambda }{2}},0,0),n_{\gamma }^{+}=(0,{(\frac{r-2M}{r})}^{\frac{-1}{2}},0,0).$$

(20)

Comparing the spacetimes defined by metrics (4) and (16) with (19), it follows that

$${(\frac{dt}{d\tau })}_{\mathsf{\Sigma }}={(e^{\frac{-\nu }{2}})}_{\mathsf{\Sigma }}={({(1-\frac{2M}{r})}^{\frac{-1}{2}})}_{\mathsf{\Sigma }},{(r)}_{\mathsf{\Sigma }}=R.$$

Moreover, the non-zero components of curvature are evaluated as

$$\begin{array}{ccc}\hfill K_{00}^{-}& =& {[\frac{-{\nu }^{\prime }{e}^{-\frac{\lambda }{2}}}{2}]}_{\mathsf{\Sigma }},K_{22}^{-}={sin}^{-2}(\theta )K_{33}^{-}={[e^{-\frac{\lambda }{2}}r]}_{\mathsf{\Sigma }},\hfill \\ \hfill K_{00}^{+}& =& {[-{(1-\frac{2M}{r})}^{\frac{-1}{2}}\frac{M}{r^2}]}_{\mathsf{\Sigma }},K_{22}^{+}={sin}^{-2}(\theta )K_{33}^{+}={[{(1-\frac{2M}{r})}^{\frac{1}{2}}r]}_{\mathsf{\Sigma }}.\hfill \end{array}$$

Employing the matching conditions on the hypersurface yields

$$e^{\nu (R)}=\frac{-2M+R}{R},e^{-\lambda (R)}=\frac{-2M+R}{R},{\nu }^{\prime }(R)=\frac{2M}{(-2M+R)R}.$$

(21)

The O’Brien and Synge [58] junction conditions (equivalent to the continuity of the second fundamental form) are given as ${(G_{\gamma \delta }{r}^{\delta })}_{\mathsf{\Sigma }}=0,$ where $r_{\gamma }$ is a unit radial vector. Using this condition in the field equations implies zero radial pressure at the hypersurface. The constants appearing in the TK metric are determined in terms of mass and radius through Equation (21) and $p_r(R)=0$ as

$$\begin{array}{ccc}\hfill a& =& \frac{24M^3-M^{2}R({\omega }_{BD}+34)-4R^5\mathcal{B}\sqrt{1-\frac{2M}{R}}+12M{R}^2}{R^2{(R-2M)}^2(2R-3M)},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill b& =& \frac{-12M^3+M^{2}R({\omega }_{BD}+20)+4R^5\mathcal{B}\sqrt{1-\frac{2M}{R}}-8M{R}^2}{R^4{(R-2M)}^2(2R-3M)},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill B& =& \frac{M}{R^2(R-2M)},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill F& =& e^{\frac{(2M-R)ln\left(\frac{R-2M}{R}\right)+M}{4M-2R}}.\hfill \end{array}$$

(25)

3. Physical Features of Compact Stars

In this section, we analyze the influence of the massive scalar field on the physical structure of strange stars through viability and stability constraints. In the present work, the numerical solutions have been obtained by considering $\mathcal{B}$ = 60 MeV/fm${}^3$, 83 MeV/fm${}^3$. These values are within the allowed limits [29]. Moreover, as per the results of Gravity Probe B experiment, all values of $m_{\Psi }>{10}^{-4}$ (in dimensionless units) are admissible [41,44]. Thus, we solve the wave equation numerically for $m_{\Psi }=0.001,0.3$ and ${\omega }_{BD}=5,8,10$ by imposing the initial conditions $\Psi (0)={\Psi }_c=$ constant and ${\Psi }^{\prime }(0)=0$. The constant ${\Psi }_c$ corresponding to the chosen values of the parameters $\mathcal{B},{\omega }_{BD}$ and $m_{\Psi }$ are presented in

Table 1. Physical properties of PSR J1614-2230 ($M=1.97\pm 0.04M_{⨀}$) with $m_{\Psi }=0.001$.

$\mathcal{B}$ = 60 MeV/fm${}^3$
${\omega }_{BD}$	${\Psi }_c$	Predicted	${\rho }_c$ (gm/cm${}^3$)	${\rho }_s$ (gm/cm${}^3$)	$p_c$ (dyne/cm${}^2$)
		Radius (km)
5	0.1	${8.8024}_{-0.1813}^{+0.1815}$	$2.3051\times {10}^{16}$	$7.2604\times {10}^{14}$	$3.9692\times {10}^{36}$
8	0.125	${11.5653}_{-0.2426}^{+0.243}$	$9.0518\times {10}^{15}$	$4.0282\times {10}^{14}$	$1.4958\times {10}^{36}$
10	0.15	${13.4483}_{-0.2869}^{+0.2878}$	$6.8805\times {10}^{15}$	$3.6403\times {10}^{14}$	$1.1523\times {10}^{36}$
$\mathcal{B}$= 83 MeV/fm${}^{\mathbf{3}}$
${\omega }_{BD}$	${\Psi }_c$	Predicted	${\rho }_c$ (gm/cm${}^3$)	${\rho }_s$ (gm/cm${}^3$)	$p_c$ (dyne/cm${}^2$)
		Radius (km)
5	0.1	${8.828}_{-0.1829}^{+0.1832}$	$2.5432\times {10}^{16}$	$6.0511\times {10}^{14}$	$4.3768\times {10}^{36}$
8	0.18	${11.6418}_{-0.2478}^{+0.2482}$	$1.2336\times {10}^{16}$	$5.2002\times {10}^{14}$	$2.1006\times {10}^{36}$
10	0.25	${13.5868}_{-0.296}^{+0.2975}$	$9.0211\times {10}^{15}$	$4.3507\times {10}^{14}$	$1.5319\times {10}^{36}$

and

Table 2. Physical properties of PSR J1614-2230 ($M=1.97\pm 0.04M_{⨀}$) with $m_{\Psi }=0.3$.

$\mathcal{B}$ = 60 MeV/fm${}^3$
${\omega }_{BD}$	${\Psi }_c$	Predicted	${\rho }_c$ (gm/cm${}^3$)	${\rho }_s$ (gm/cm${}^3$)	$p_c$ (dyne/cm${}^2$)
		Radius (km)
5	0.1	${8.8024}_{-0.1813}^{+0.1815}$	$1.5412\times {10}^{16}$	$1.2272\times {10}^{15}$	$2.4962\times {10}^{36}$
8	0.125	${11.5653}_{-0.2426}^{+0.243}$	$9.3622\times {10}^{15}$	$8.3976\times {10}^{14}$	$1.509\times {10}^{36}$
10	0.15	${13.4483}_{-0.2869}^{+0.2878}$	$6.4885\times {10}^{15}$	$9.9402\times {10}^{14}$	$9.4198\times {10}^{35}$
$\mathcal{B}$= 83 MeV/fm${}^{\mathbf{3}}$
${\omega }_{BD}$	${\Psi }_c$	Predicted	${\rho }_c$ (gm/cm${}^3$)	${\rho }_s$ (gm/cm${}^3$)	$p_c$ (dyne/cm${}^2$)
		Radius (km)
5	0.1	${8.828}_{-0.1829}^{+0.1832}$	$2.5687\times {10}^{16}$	$9.3087\times {10}^{14}$	$4.4261\times {10}^{36}$
8	0.18	${11.6418}_{-0.2478}^{+0.2482}$	$1.3305\times {10}^{16}$	$1.3258\times {10}^{15}$	$2.1018\times {10}^{36}$
10	0.25	${13.5868}_{-0.296}^{+0.2975}$	$9.7515\times {10}^{15}$	$1.0625\times {10}^{15}$	$1.5094\times {10}^{36}$

. The quantities with subscripts s and c are derived at the surface and center of the star, respectively. The results are displayed graphically for PSR J1614-2230 ($M=1.97\pm 0.04M_{⨀}$) [59].

Literature survey reveals that anisotropy vanishes at $r=0$ and then increases to attain a maximum value at the surface of the compact star [12,13,60]. Deb et al. [34] evaluated the extremum of anisotropy at the boundary and showed that for a viable and stable spherical structure anisotropy must be maximum at $r=R$. Adopting the procedure in [34], we predict radius of the strange star candidate PSR J1614-2230 by maximizing anisotropy at the boundary. The metric potentials must be well-behaved and positive throughout the interior region to ensure a singularity-free geometry [61]. As shown in Figure 1, the TK metric functions are monotonically increasing as well as regular functions of the radial coordinate. Energy density and pressure play a dominant role in determining the behavior of extremely dense strange stars. The matter variables must be positive and decrease monotonically towards the surface of the star. Figure 2 and Figure 3 show that the state determinants are maximum at the center and decrease away from it for the considered values of $m_{\Psi }$, coupling parameter and bag constant.

The presence of radial/transverse pressure leads to anisotropy ($\Delta =p_{\perp }-p_r$) within the system. The anisotropy is negative when $p_r>p_{\perp }$ and positive otherwise. Moreover, particles are tightly packed together in dense celestial objects which limits the particles’ movement in the radial direction. Consequently, the radial force or pressure is less than the tangential force leading to positive anisotropy. Thus, the positive anisotropy generates an outward repulsive force increasing the stability and compactness of the star. The anisotropy of the current setup, displayed in Figure 4, is calculated through Equations (14), (15) and (22)–(25). The anisotropy is positive throughout the stellar region for the considered values of $m_{\Psi },{\omega }_{BD}$ and $\mathcal{B}$.

3.1. Energy Conditions

The anisotropic configuration is physically viable if it complies with four constraints, i.e., null (NEC), weak (WEC), strong (SEC) and dominant (DEC). In the background of MBD gravity, these conditions are stated in terms of following inequalities [62].

$$\begin{array}{c}\mathrm{NEC}:\rho \ge 0,\hfill \\ \mathrm{WEC}:\rho +p_r\ge 0,\rho +p_{\perp }\ge 0,\hfill \\ \mathrm{SEC}:\rho +p_r+2p_{\perp }\ge 0,\hfill \\ \mathrm{DEC}:\rho -p_r\ge 0,\rho -p_{\perp }\ge 0.\hfill \end{array}$$

The positive trend of state determinants presented in Figure 2 and Figure 3 readily satisfies the first three inequalities. The plots of DEC in Figure 5 show that DEC is violated near the boundary of the surface for $m_{\Psi }=0.001$. However, the setup corresponding to $m_{\Psi }=0.3$ satisfies the DEC except for ${\omega }_{BD}=5$ as shown in Figure 6. The critical values of $m_{\Psi }$ for which the stellar model becomes viable are mentioned in

Table 3. Critical values of $m_{\Psi }$ for different values of $\mathcal{B}$ and ${\omega }_{BD}$ corresponding to PSR J1614-2230 ($M=1.97\pm 0.04M_{⨀}$).

$\mathcal{B}$ = 60 MeV/fm${}^3$		$\mathcal{B}$ = 83 MeV/fm${}^3$
${\omega }_{BD}$	$m_{\Psi }$	${\omega }_{BD}$	$m_{\Psi }$
5	0.7	5	0.65
8	0.3	8	0.25
10	0.2	10	0.15

.

3.2. Effective Mass, Compactness and Redshift

Gravitational mass (m) is one of the key features that determine the structure and compactness of stellar objects. The mass of a spherical celestial object (measured using Kepler’s law when a satellite orbits the star) is evaluated as

$$m=\frac{1}{2}{\int }_0^r\rho r^{2}dr.$$

(26)

The mass of the anisotropic stellar model is obtained by numerically solving the above equation along with the wave equation under the condition $m(0)=0$. The size and mass of a cosmic body provide a measure of its compactness ($\mu (r)=\frac{m}{r}$). Buchdal [63] proposed that a compact object is stable if $\mu (r)<\frac{4}{9}$. The value of the compactness function at the surface of the star is shown in

Table 4. Values of ${\mu }_s$ and $Z_s$ for different values of $m_{\Psi },\mathcal{B}$ and ${\omega }_{BD}$ corresponding to PSR J1614-2230 ($M=1.97\pm 0.04M_{⨀}$).

	$m_{\Psi }=0.001$
	$\mathcal{B}$ = 60 MeV/fm${}^3$		$\mathcal{B}$ = 83 MeV/fm${}^3$
${\omega }_{BD}$	${\mu }_s$	$Z_s$	${\mu }_s$	$Z_s$
5	0.0291	0.0303	0.0326	0.0338
8	0.0311	0.0326	0.0441	0.0468
10	0.0386	0.0406	0.0522	0.0571
	${\mathbf{m}}_{\Psi }=\mathbf{0}.\mathbf{3}$
	$\mathcal{B}$ = 60 MeV/fm${}^3$		$\mathcal{B}$ = 83 MeV/fm${}^3$
${\omega }_{BD}$	${\mu }_s$	$Z_s$	${\mu }_s$	$Z_s$
5	0.0354	0.0370	0.0360	0.0381
8	0.0387	0.0412	0.0613	0.0667
10	0.0474	0.0514	0.0976	0.1146

for $m_{\Psi }=0.001$ and $m_{\Psi }=0.3$. In both cases, the compactness parameter increases corresponding to an increase in $m_{\Psi },{\omega }_{BD}$ or $\mathcal{B}$. Moreover, it attains values less than $\frac{4}{9}$. Under the influence of a star’s gravitational field, the electromagnetic radiation loses some of the energy resulting in an increase in its wavelength, i.e., the radiation is redshifted. The effect of the gravitational force is measured through a redshift parameter defined as

$$Z=\frac{1}{\sqrt{1-2u(r)}}-1.$$

For an anisotropic source, the redshift must not exceed the value $5.211$ [64]. The redshift parameter obeys the limit proposed in [64] (refer to Table 4).

3.3. Stability of Anisotropic Stellar Model

In this section, the anisotropic model is checked for stability via three approaches: causality condition [65], Herrera’s cracking approach [66] and adiabatic index. Causality condition enforces the relation between cause and effect within the stellar region. The condition holds when sound waves travel at a speed less than that of light, i.e., $0<v_r^2<1$ and $0<v_{\perp }^2<1$, where $v_r^2=\frac{d{p}_r}{d\rho }$ and $v_{\perp }^2=\frac{d{p}_{\perp }}{d\rho }$ are the radial and tangential components of sound speed, respectively. On the other hand, the cracking approach provides a criterion for potentially stable systems. The phenomenon of internal perturbations changing the sign of the radial forces is known as cracking. According to this approach, the condition for a region free from cracking is $0<|v_{\perp }^2-v_r^2|<1$. Figure 7 indicates that the system corresponding to $m_{\Psi }=0.001$ is unstable as tangential velocity becomes negative near the stellar surface. However, the system is compatible with Herrera’s cracking approach for higher vales of the coupling parameter. The plots in Figure 8 adhere to the criteria of causality and cracking approaches. Thus, a stable structure is obtained for $m_{\Psi }=0.3$.

If an increase in density results in an effective increase in pressure, the system obeys a stiff EoS. A structure associated with a stiff EoS is harder to compress and more stable as compared to a setup with a soft EoS. The stiffness of EoS is measured through adiabatic index given as

$$\Gamma =\frac{p_r+\rho }{p_r}\frac{d{p}_r}{d\rho }=\frac{p_r+\rho }{p_r}{v}_r^2.$$

The lower limit of adiabatic index for a stable anisotropic distribution is $\frac{4}{3}$ [67]. The adiabatic index for different values of $m_{\Psi },{\omega }_{BD}$ and $\mathcal{B}$ is displayed in Figure 9. The anisotropic configuration is stable according to this criterion when

• $m_{\Psi }=0.001,{\omega }_{BD}=5,8,\mathcal{B}$ = 60 MeV/fm${}^3$.
• $m_{\Psi }=0.001,{\omega }_{BD}=8,10,\mathcal{B}$ = 83 MeV/fm${}^3$.
• $m_{\Psi }=0.3,{\omega }_{BD}=8,\mathcal{B}$ = 60 MeV/fm${}^3$.
• $m_{\Psi }=0.3,{\omega }_{BD}=8,10,\mathcal{B}$ = 83 MeV/fm${}^3$.

4. Concluding Remarks

Astrophysicists have studied the dynamics of stellar structures and their remnants to gain insight into the mechanism of the universe. One of the hypothesized remnants are the quark stars that emerge from the collapse of neutron stars. The interaction among constituent quarks (up, down and strange) of these dense celestial structures is represented by the MIT bag model. In this paper, we have investigated the existence and stability of quark stars in the framework of MBD gravity. For this purpose, we have constructed field equations using well-behaved TK metric potentials and MIT bag model. The unknown constants appearing in the anisotropic model have been specified by matching the interior TK spacetime to exterior Schwarzschild spacetime. The solution has been determined by choosing the potential function as $V(\Psi )=\frac{1}{2}{m}_{\Psi }^2{\Psi }^2$. Moreover, the wave equation has been solved numerically for $m_{\Psi }=0.001,0.3,{\omega }_{BD}=5,8,10$ and $\mathcal{B}$ = 60 MeV/fm${}^3$, 83 MeV/fm${}^3$. Finally, we have explored the behavior of matter variables as well as viability and stability of the resulting model in the presence of a massive scalar field.

The observed mass of the strange star candidate PSR J1614-2230 has been used to maximize the anisotropy at $r=R$ and predict the radius of the stellar structure. The predicted radius presented in Table 1 and Table 2 increases with an increase in the values of the bag constant and coupling parameter. The behavior of matter variables shows that the cosmic object becomes denser as the value of the bag constant increases. Moreover, these variables attain a maximum value at the center and decrease towards the surface. The anisotropy increases with an increase in the values of the bag constant. However, an increase in $m_{\Psi }$ or ${\omega }_{BD}$ reduces the anisotropic pressure within the star. Three energy conditions (NEC, WEC, SEC) are readily satisfied for all values of the involved parameters. The model violates the DEC near the boundary for $m_{\Psi }=0.001$. On the other hand, the anisotropic structure becomes viable for $m_{\Psi }=0.3$ except for ${\omega }_{BD}=5$.

We have also analyzed the effective mass, compactness and redshift parameter of the constructed model. An increase in the mass of scalar field and value of the bag constant leads to massive and more compact structures with higher redshift parameter (refer to Table 4). It is worth mentioning here that compactness and redshift parameters obey their respective upper bounds [63,64]. The stability of the stellar model has been investigated through three criteria. The tangential velocity becomes negative at the surface of the star when $m_{\Psi }=0.001$ for the chosen values of the bag constant and coupling parameter. However, the model is stable according to causality and cracking criteria for $m_{\Psi }=0.3$. Finally, the plots in Figure 9 indicate that the adiabatic index is less than $\frac{4}{3}$ for some values of the parameters. It is concluded that viable and stable solutions have been obtained in MBD theory for $m_{\Psi }=0.3$ when

• ${\omega }_{BD}=8$ and $\mathcal{B}$ = 60 MeV/fm${}^3$.
• ${\omega }_{BD}=8,10$ and $\mathcal{B}$ = 83 MeV/fm${}^3$.

It is worthwhile to mention here that all our results reduce to GR for ${\omega }_{BD}\to \infty$.