Quark stars in massive gravity might be candidates for the mass gap objects

We have investigated the structural properties of strange quark stars (SQSs) in a modified theory of gravity known as massive gravity. In order to obtain the equation of state (EOS) of strange quark matter, we have employed a modified version of the Nambu-Jona-Lasinio model (MNJL) which includes a combination of NJL Lagrangian and its Fierz transformation by using weighting factors ($1-\alpha $) and $\alpha$. Additionally, we have also calculated dimensionless tidal deformability ($\Lambda$) in massive gravity. To constrain the allowed values of the parameters appearing in massive gravity, we have imposed the condition $\Lambda_{1.4 {M}_{\odot }}\lesssim580 $. Notably, in the MNJL model, the value of $\alpha$ varies between zero and one. As $\alpha$ increases, the EOS becomes stiffer, and the value of $\Lambda$ increases accordingy. We have demonstrated that by softening the EOS with increasing the bag constant, one can obtain objects in massive gravity that not only satisfy the constraint $\Lambda_{1.4 {M}% _{\odot }}\lesssim580$, but they also fall within the unknown mass gap region ($2.5{M}_{\odot}-5{M}_{\odot }$). To establish that the obtained objects in this region are not black holes, we have calculated Schwarzschild radius, compactness, and $\Lambda_{{M_{TOV}}}$ in massive gravity.


I. INTRODUCTION
The gravitational wave (GW) events offer new insights into compact stars.The binary merger GW170817 [1] and its electromagnetic counterpart [2] have led to new constraints on the maximum mass of neutron stars (NSs).Based on GW170817, the upper bound of M T OV for NSs is predicted as ∼ (2.3 − 2.4)M ⊙ [3,4].However, there are other observations of pulsars and binary mergers with masses greater than these values such as PSR J0952-0607 [5], the secondary component of GW190814 [6,7] and the remnants of GW170817 [8] and GW190425 [9] which fall within the unknown mass gap region (2.5M ⊙ ≲ M ≲ 5M ⊙ ).This region was called a mass gap because analysis of observations of X-ray binaries revealed a small number of compact objects in this region [10][11][12].Recent observed systems within this interval, show a relative gap instead of an absolute gap [13].There are arguments that the presence or absence of mass gap objects can constrain the properties of the core-collapse supernova engine and the compact objects in the mass gap region can be created by merging lighter compact objects [14,15].There are two crucial reasons why objects in the mass gap region are unlikely to be neutron stars.The first reason is that for non-rotating NSs, the equation of state (EOS) must be very rigid to get such massive objects (the speed of sound is greater than 0.6 of the speed of light) [16,17].Such EOSs contrast dimensionless tidal deformability (Λ) constraints obtained from GW178017, which require softer EOSs.The second reason is that for a rotating NS, the rotation rate should be close to the keplerian limit [16,[18][19][20].Now, the question arises whether these objects are the smallest black holes or other forms of compact stars, like hybrid stars and self-bounded strange quark stars (SQSs).
Theoretically, strange quark matter (SQM) is formed in two classes of compact stars.i) Hybrid stars with a quark core and hadronic shells [21][22][23][24][25][26] and ii) SQSs composed of pure quark matter from the core up to the upper layers of the star [27][28][29][30][31]. Ever since Trazawa [32], Witten [33], and Bodmer [34] proposed SQM as the ground state of QCD, until now when the GWs have found so many binaries, the study of SQM in compact stars has always been of interest.Many studies on hybrid and quark stars have focused on recent discoveries of GWs.Refs.[35][36][37] study hybrid stars with constraints obtained from GW190814, GW170817, and PSR J0030+0451 [38,39].In Ref. [40], the properties of SQM in the bag model and non-newton gravity are investigated by the constraints obtained from GW170817 and PSR J0740+6620.In Ref. [41], the EOS of quark matter is constrained by using the bayesian statistical model and the mass and radius measurements of PSR J0030+0451 [42].Refs.[7,16,43] investigate whether the secondary component of GW190814 could be a SQS.In Refs.[44,45] the structural features of the SQSs have also been investigated using the data obtained from GWs.
GR is a successful theory of gravity.However, at large scales, it cannot explain why our universe is undergoing an accelerated cosmic expansion.This is one of the reasons why GR needs to be modified.Among various mod-ified theories of gravity, massive gravity can describe the late-time acceleration without considering dark energy [46][47][48][49][50].This theory of gravity modifies gravitational effects by weakening them at a large scale compared to GR, which allows the universe to accelerate.However, its predictions at small scales seem to be the same as those of GR.Massive gravity will result in the graviton having a mass of m g , which in the absence of the mass of the graviton (i.e., m g → 0), this theory reduces to GR.In Ref. [51], it has been indicated that the mass of a graviton is very small in the usual weak gravity environments but becomes much larger in the strong gravity regimes, such as black holes and other compact objects.In addition, recent observations by the advanced LIGO/Virgo collaborations have put a tight bound on the graviton mass [52,53].Also, other theoretical and empirical limits on the graviton's mass exist [54][55][56][57].In the astrophysics context and by considering massive gravity, the properties of compact objects such as relativistic stars [58,59], neutron stars [60,61], and white dwarfs [62], have been studied.In 2011, de Rham, Gabadadze and Tolley (dRGT) introduced a ghost-free theory of massive gravity [63,64] which is known dRGT massive gravity.This theory uses a reference metric to construct massive terms (see Refs. [63][64][65], for more details).These massive terms are inserted in the action to provide massive gravitons.Then, Vegh extended the dRGT massive gravity theory using holographic principles and a singular reference metric in 2013 [66], which is known as dRGT-like massive gravity.Recently, this extended theory of dRGT massive gravity attracted many people in various viewpoints of physics.
In this paper, we have used a modified version of the Nambu-Jona-Lasinio (MNJL) model to obtain the structural properties of SQS in massive gravity.In GR, the soft EOSs give small masses for quark stars, and the stiff EOSs do not satisfy the Λ constraint.We show that in massive gravity, it is possible to have quark stars that not only have masses higher than those in GR, but also satisfy the GW observational constraints well.We have used Λ constraint (Λ 1.4M⊙ ≲ 580) [67] to restrict the masses of SQSs in massive gravity.

II. MODIFIED NJL MODEL
We have used a modified self-consistent version of the (2 + 1)−flavor NJL model in mean field approximation [68] to derive the EOS.The modified Lagrangian is a combination of the original NJL Lagrangian and its Fierz transformation [69] by using weighting factors (1−α) and α as follows [68] where L N JL is the original Lagrangian in flavor space without Fierz transformation and L F is the Fierz transformation of L N JL in color space [68].ψ denotes a quark field with three flavor (u, d, s) and three color degrees of freedom.m is the corresponding mass matrix in flavor space.In addition, γ µ s (µ = 0 −→ 3) and λ i s (i = 1 −→ 8) are Dirac and Gell-Mann matrices, respectively.γ 5 = iγ 0 γ 1 γ 2 γ 3 and λ 0 = 2 3 I, where I is the identity matrix.Similar to λ 0 , λ C 0 is equal to but in color space [68].Also, G and K are the fourfermion and six-fermion interaction coupling constants, respectively.The parameter α varies from zero to one.By setting α = 0, the Lagrangian 1 converts to original NJL Lagrangian and for α = 1, it converts to Fierz transformation of L N JL .In Ref. [68], the authors have used the Lagrangian (1) and have shown that this Lagrangian gives the EOSs with M T OV ≃ 2M ⊙ in GR which satisfy structural constraints (mass, radius and Λ) obtained by GWs [68].They have used the first analysis of Λ with the upper bound of Λ 1.4M⊙ ≤ 800 [1] and have obtained a SQS with M T OV ≃ 2.057M ⊙ and Λ 1.4M⊙ = 634.In this paper, we show that by using the Lagrangian (1) and massive gravity, SQSs can be obtained with masses in gap region (2.5M ⊙ ≲ M ≲ 5M ⊙ ) and simultaneously, respect to the improved analysis of Λ with the upper bound of Λ 1.4M⊙ ≲ 580 [67].
In the following, we use the Lagrangian (1) and obtain the EOS of SQM in mean field approximation.As the first step, we calculate the quark number densities for different choices of α (to know the details of calculations, see the Ref. [68]).Fig. 1 shows quark number densities of up, down, and strange flavors as a function of chemical potential for different values of α.As one can see from the figure 1, the quark number densities of up and down quarks are zero for µ ≲ 200M eV , and the strange quark number density is zero for µ ≲ 320M eV .This figure shows that the quark number densities decrease with increasing α.Consequently, such behavior also occurs for baryon number density (for more discussions, see Ref. [68]).
The following conditions should be imposed to obtain the EOS of SQM, which are: The first two lines of Eq. ( 2) come from beta equilibrium and charge neutrality conditions, respectively.In the above conditions, µ i s (i = u, d, s) are quark chemical potentials with flavor i, and µ e is the electron chemical potential.n i s (i = u, d, s) are quark number densities with flavor i.In addition, n B is the baryon number density.Furthermore, n e = µ 3 e /(3π 2 ) is the electron number density.Then the pressure and the energy density can be obtained by the following relations where P 0 is the vacuum pressure, which is taken as the vacuum bag constant (−B).We consider two values for B including B 1/4 = 117M eV and B 1/4 = 130M eV .For these values of B, the energy per baryon of three-flavor quark matter is lower than that of two-flavor quark matter for α ≤ 0.94 [68].This condition is necessary for SQM to be the true ground state of the strong interaction [70].
In the following, we present our results for the structural properties of SQS in massive gravity for different values of B and α.To this end, we first explain the modified TOV equation in massive gravity and then derive the differential equations to calculate the Λ in this theory of gravity.

III. MODIFIED TOV EQUATION IN MASSIVE GRAVITY
Here, we consider dRGT-like massive gravity, which is introduced in Ref. [66].Notably, the dRGT-like massive gravity is similar to the dRGT massive theory of gravity, which is proposed by de Rham-Gabadadze-Tolley (dRGT) in 2011 [63].The dRGT-like massive gravity's action is given by [66] where R and m g are the Ricci scalar and the graviton mass, respectively.It is notable that κ 2 = 8π because we consider G = c = 1.Also, f and g are fixed symmetric and metric tensors, respectively.In the above action, I matter is related to the action of matter.Moreover, c i 's are free parameters of the massive theory, which are arbitrary constants.Notabley, their values can be determined according to theoretical or observational considerations [71][72][73].U i 's are symmetric polynomials of the eigenvalues of 4 × 4 matrix K µ ν = √ g µα f αν (where g µν is the dynamical (physical) metric, and f µν is the auxiliary reference metric, needed to define the mass term for gravitons) where they can be written in the following form where U i−y = 1, when i = y.Also, K stands for matrix square root, i.e., , and the rectangular bracket denotes the trace [K] = K µ µ .We consider the four-dimensional static spherical symmetric spacetime dr 2 e −2λ(r) −r 2 h ij dx i dx j , where i, j = 1, 2.
(6) where e 2Φ(r) and e −2λ(r) are unknown metric functions, and h ij dx i dx j = dθ 2 + sin 2 θdφ 2 .The equation of motion of dRGT-like massive gravity is given in the following form where G ν µ is the Einstein tensor, and T ν µ denotes the energy-momentum tensor which comes from the variation of I matter and χ µν is the massive term with the following explicit form Considering the quark stars as a perfect fluid, the nonzero components of the energy-momentum tensor are given by [61] T We follow the reference metric f µν for four-dimensional static spherical symmetric spacetime as it is considered for studying neutron stars and white dwarfs in the following form [61,62] in which J is known as a parameter of the reference metric, which is a positive constant.The explicit functional forms of U i 's are given by [61] Using the field equation and reference and dynamical metrics, the metric function e −2λ(r) , is obtained in the following form [61] in which M (r) = R 0 4πr 2 ϵ(r)dr or dM (r) dr = 4πr 2 ϵ(r) .
By definition 11), the modified TOV equation in massive gravity is given by [61] Here, we obtain Λ in massive gravity.To do this, we need to get a differential equation of metric function (H(r)) in massive gravity.In Appendix A, we have shown how we obtained this equation in massive gravity.Also, to check the correctness of the calculations, we have investigated the limit C 1 = C 2 = 0 and shown that in this limit, the differential equation obtained in massive gravity is converted to the differential equation in GR.The differential equation to get H(r) in massive gravity is where β = dH(r)/dr, and A i s are as follows where B 1 , B 2 , B 3 , and B 4 are also, f is defined as dϵ/dP .In addition, Λ is then obtained by where k 2 is dimensionless tidal Love number for l = 2, which is given as where σ = M/R and y is defined as y = Rβ(R)/H(R) − 4πR 3 ϵ 0 /M in which ϵ 0 is the energy density at the surface of the quark star [68].By simultaneously solving equations ( 12) and ( 13) along with dM/dr = 4πr 2 ϵ, we obtain the curves of star mass in terms of radius and Λ in terms of mass.The results are presented in the next section.

A. Modified Schwarzschild Radius
It is clear that by considering dRGT-like massive gravity, the Schwarzschild radius is modified.Using Eq. ( 11), we can get the Schwarzschild radius (R Sch ) for this theory of gravity.To find the modified Schwarzschild radius, we solve S(r) = 0. So, the modified Schwarzschild radius for dRGT-like massive gravity is given [61] where M = M (r = R) is the star's total mass (or the mass of the star on its surface).

B. Compactness
The compactness of a spherical object may be defined by the ratio of mass to radius of that object which may be indicated as the strength of gravity.In the following, we calculate the thermodynamic properties of SQM (EOS, speed of sound, and adiabatic index) for different values of α and B. We investigate the structural properties of SQS (mass, radius, dimensionless tidal deformability, compactness, and Schwarzschild radius) for a stiff EOS (B 1/4 = 117M eV ), and a soft EOS (B 1/4 = 130M eV ) in massive gravity.The value of Λ depends not only on the EOS but also on the gravity used.
In GR, the stiff EOSs cause the Λ constraint to be lost and the soft EOSs can not cover the objects fall within in the mass gap region.We demonstrate that in massive gravity, the soft EOSs can not only lead to masses in the mass gap region but also satisfy the Λ constraint well.

IV. EQUATION OF STATES OF SQM FOR
In this section, we present our results for the thermodynamic properties of SQM in a modified NJL model for B 1/4 = 117M eV and different values of α.As it was mentioned in section II, by choosing this value of bag constant, the parameter α should be equal or less than 0.94 to establish the stability condition of SQM.We derive the EOSs of SQM for different values of α including α = 0.5, α = 0.8, and α = 0.94.The results are presented in Fig. 2.This figure shows that the EOSs become stiffer by increasing the α parameter.We will show in section IV A that this behavior causes the mass of the star to increase.But it is possible that the structural features of these stars do not satisfy the constraints obtained from GW observations.Therefore, at the same time as the mass of the star increases, such constraints should also be checked.To check the causality, the speed of sound in terms of energy density is plotted in Fig. 3.The plots indicate that causality is established well for all α values.Also, with an increase in the value of α, the speed of sound increases, which means that the EOS becomes stiffer and subsequently, the mass of the star becomes greater.As a measure of dynamical stability, the adiabatic index Γ = dP dϵ P +ϵ P has been presented as a function of energy density in Fig. 4. It is known that Γ should be higher than 4/3 for dynamical stability [61,[74][75][76]78].As one can see from this figure, this condition is satisfied for all values of α.

A. Structural Properties of SQS in Massive
Gravity for B 1/4 = 117M eV In this section, we obtain the structural properties of SQS, including the M − R and Λ − M relations in dRGT- like massive gravity.We then compare the results with the constraints derived from astronomical observations.In addition, we obtain modified Schwarzschild radius and compactness of quark stars in this theory of gravity.Practically, the modified TOV equation ( 12) is solved simultaneously with the equation, M (r) = R 0 4πr 2 ϵ(r)dr.To solve these equations, we have to consider the modified NJL EoS with the obtained pressure at the center of the star (i.e.P (r = 0) = P c , which P c is the central pressure of the quark star and given by the modified NJL EoS).Notably, the mass of the star at the center of the star is zero, i.e.M (r = 0) = 0.Then, we continue to solve the equations up to the star's surface, where the pressure is zero (P (r = R) = 0).The result is the mass and corresponding radius of the star (i.e.M = M (R)).We do this for different values of central pressures allowed by the mentioned EOS.
1. α = 0.5 Now, we investigate the structural properties of SQS in massive gravity.We first use the obtained EOSs in section IV and the modified TOV equation (12) to derive the mass-radius (M − R) relation and Λ − M relation in massive gravity.As we saw in section IV, we used three different values of α including 0.5, 0.8, and 0.94 for the EOSs in the MNJL model.Now, we first examine the case α = 0.5.The equation (12) shows that there are two parameters, including C 1 and C 2 , in massive gravity (these parameters are related to the dRGT-like massive gravity's parameters).In the beginning, we set a fixed value for C 1 and investigate the structural features of SQS for different values of C 2 .Figures 5 and 6 show the M − R and Λ − R diagrams, respectively, for α = 0.5, C 1 = 10 −5 , and different values of C 2 .As one can see from Fig. 5, there are different colored regions in Fig. 5.The gray and orange areas are the M − R constraints from GW170817.The black boundary corresponds to PSR J0740+6620 [79].The pink area represents the secondary component of GW190814 [7], and the green area shows the remnant mass of GW190425 [9].As one can see from this figure, we have not drawn graphs for the values |C 2 | > 0.14, because Fig. 6 shows for these values of C 2 , the Λ condition (Λ 1.4M⊙ ≲ 580) is not satisfied.In Fig. 6, the horizontal black line shows the Λ = 580, and the vertical red line represents M = 1.4M ⊙ .As one can see from this figure, the constraint Λ 1.4M⊙ ≲ 580, is satisfied for |C 2 |≲ 0.14 which corresponds to the M T OV ≲ 2.13M ⊙ .Indeed, the maximum mass in this case cannot be more than 2.13M ⊙ .More details for structural properties of SQS for B 1/4 = 117M eV , α = 0.5, C 1 = 10 −5 , and different values of C 2 are presented in Table .I. Previously, in the Ref. [68] the value of M T OV for B 1/4 = 117M eV and α = 0.5 has been calculated in GR.The obvious difference between their result (M T OV = 1.70M ⊙ ) with the result obtained in this paper shows the important effect of the massive gravity in the mass of the SQS.Table.I demonstrates that by increasing the absolute value of C 2 , the values of M T OV , Λ 1.4M⊙ , and the compactness of the star are increased too.Additionally, The results calculated for R Sch show that the obtained masses can not be black holes.Now, we turn to the next case by choosing It is shown in Refs.[61,76] that changing the C 1 parameter has a negligible effect on the structural properties of the star.Therefore, in this section, we show this feature and for other values of α and the bag constant, we only change the value of C 2 .Figure 7  The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135 [77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451 [38,39].
changing the C 1 parameter has an insignificant effect on the maximum mass of the star.We have not plotted the results for values C 1 ≳ 1.3 × 10 −5 , because as the figure 8 shows, the condition Λ   Here, we increase the value of α and investigate its effect on the structure properties of the SQS in dRGTlike massive gravity.The figures.9 and 10 show the M − R and Λ − M diagrams for α = 0.8, C 1 = 10 −5 , and different values of C 2 , respectively.The green and black hatched regions represent pulsars J2215+5135 [77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.
The brown and the purple regions show two different reports of the pulsar J0030+0451 [38,39].
Comparison of figures 9 and 5 shows that the maximum mass of the star increases with the increase of α value.This result is obvious because we showed in Fig. 2 that the increase of α makes the EOS stiffer.However, Fig. 10 shows that this feature increases Λ and the constraint Λ 1.4M⊙ ≲ 580 is lost for |C 2 | ≳ 10 −2 which corresponds to the M T OV ≲ 2.33M ⊙ .Indeed, the maximum mass in this case cannot be more than 2.33M ⊙ .The structural properties of SQS for B 1/4 = 117M eV , α = 0.8, C 1 = 10 −5 , and different values of C 2 are presented in Table .III.
There are the same behaviors for the compactness and the modified Schwarzschild radius.They increase by increasing the value of |C 2 |.Also, the radii of these compact objects are more than the modified Schwarzschild radii (i.e., R > R Sch ).For further investigation, we check the effect of massive gravity on the structural properties of SQS for α = 0.94.We do not consider values of α greater than 0.94.Because, as mentioned in section II, the stability condition of SQM will be lost for α ≳ 0.94 [68].BY setting B 1/4 = 117M eV , α = 0.94, C 1 = 10 −5 , and different values of C 2 , we obtain the Fig. 11, and 12 as the M − R and Λ − M diagrams, respectively.Here, we also see that with the increase of α, the maximum gravitational mass increases, but as the EOS becomes stiffer, the Λ constraint is no longer satisfied.Table .IV shows the structural properties of SQS for B 1/4 = 117M eV , α = 0.94, C 1 = 10 −5 , and different values of C 2 .Although the values obtained for the M T OV for α = 0.94 are higher than those with α = 0.5 and α = 0.8, Table.IV demonstrates that none of the results satisfy the constraint Λ 1.44M⊙ ≲ 580.As an important result of our calculation, we could not describe the mass gap region and observational objects such as PSRJ0740+6620, the secondary component of GW190814, and the remnant mass of GW190425 by SQS when we considered the mentioned EOSs of SQM with B 1/4 = 117M eV and different values of α.In other words, the maximum mass of SQS with the mentioned EOSs in dRGT-like massive gravity is in the range M T OV < 2.5M ⊙ .
To obtain massive SQSs that satisfy the Λ constraint, we make the EOS soft.Because the value of Λ depends on both the EOS and the gravity used, simultaneously.Therefore, in the next section, we change the value of the B constant (as we expected) and investigate the results for thermodynamic properties of SQM and structural properties of SQS in dRGT-like massive gravity.

V. EQUATION OF STATES OF SQM FOR
In the previous section, we saw that by increasing the value of α, the EOS becomes stiffer, and subsequently the mass of the star increases, but the constraint Λ 1.44M⊙ ≲ 580 is not satisfied for the SQSs fall within the mass gap region and consequently, a small range of the C 2 parameter is reachable.In this section, we increase the value of the bag constant, B, to make the EOS softer, and the condition of Λ to be satisfied in a larger range of C 2 parameter.Figure 13 shows a comparison between two EOSs obtained from two different bag constants (B 1/4 = 117M eV and B 1/4 = 130M eV ).From this figure, it can be seen that the EOSs with the same value of α have the same slope and consequently have the same speed of sound.But the diagrams with B 1/4 = 130M eV lie lower than those with B 1/4 = 117M eV .It shows that the EOSs with B 1/4 = 130M eV , are softer than those with B 1/4 = 130M eV .Although this result is obvious, we will indicate that it has a significant effect on the structural properties of SQS in massive gravity.This makes it possible to obtain quark stars that fall within the mass gap region and respect to the constraint Λ 1.4M⊙ ≲ 580, simultaneously.Now, we set B 1/4 = 130M eV and investigate its effect on the structural properties of SQS in massive gravity.In this situation, we consider another constraint too.We assume that the value of Λ for M T OV should not be zero.Because otherwise the obtained mass is a black hole [80].

A. Structural Properties of SQS in Massive
Gravity for B 1/4 = 130M eV In the following, we study the structural properties of SQS in the presence of dRGT-like massive gravity for Here, we derive the structural properties of SQS in massive gravity for α = 0.5 and B 1/4 = 130M eV .Similar to the previous section, we first set a fixed value for C 1 and investigate M − R and Λ − M diagrams for different values of C 2 .Figures 15 and 16 show M − R and Λ − M diagrams, respectively.We can see from the fig.16 that all the diagrams satisfy the constraint Λ 1.4M⊙ ≲ 580.The results show that by softening the EOS, the masses obtained in massive gravity not only satisfy the Λ 1.4M⊙ ≲ 580 constraint, but also some of them are located in the mass gap region.As we can see from these figures, we did not use the values |C 2 | ≳ 0.84, because they led to zero value for Λ M T OV .The results for structural properties are given in Table .V.
As an important result, our findings indicate that the maximum mass of SQS can be in the range M T OV ≲ 4M ⊙ .Indeed, by using the MNJL model with B 1/4 = 130M eV , and α = 0.5, SQSs fall within the mass gap region.
2. α = 0.8 As a further investigation, we change the value of α.Figures.17    The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135 [77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451 [38,39].The light blue band denotes the remnant mass of GW190425.the increase of α, which leads to the stiffer EOS, the constraint Λ 1.4M⊙ ≲ 580 is still maintained.Here, the obtained masses lie in a larger area of the mass gap region compared to the previous case (B  The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135 [77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451 [38,39].The light blue band denotes the remnant mass of GW190425. It is notable that by considering the MNJL model with B 1/4 = 130M eV , and α = 0.8, the maximum mass of SQSs can be in the range 1.68M ⊙ ≲ M T OV ≲ 4.2M ⊙ , which fall within the mass gap region.for C 1 = 10 −5 and different values of C 2 , respectively.Figure 20 shows that for values |C 2 | ≳ 0.54, the constraint Λ 1.4M⊙ ≲ 580 is no longer satisfied.For this reason, compared to the case α = 0.8, a smaller area of the C 2 parameter is accessible, which results from the sharpness of the EOS.The structural properties of SQS for B 1/4 = 130M eV , α = 0.94, C 1 = 10 −5 and different values of C 2 are presented in Table .VII.Briefly, the maximum mass of SQSs with B 1/4 = 130M eV , and α = 0.94, cannot be more than 3.45M ⊙ .However, the mentioned values of the parameters of our system can satisfy the mass gap region (see Table .VII, for more details).Besides, there are the same behaviors of the compactness and the Schwarzschild radius when the value of |C 2 | increases.In other words, the strength of gravity augments by increasing the value of |C 2 |.
As a final point, by considering the MNJL model with B 1/4 = 130M eV and different values of α (0.5 < α < 0.94), we indicated that the mass of SQSs could be in the mass gap region.So, it seems reasonable to consider SQSs in massive gravity as candidates for the mass gap region.

VI. CONCLUSION
We investigated the structural properties of quark star in massive gravity.In this study, we considered three-  flavor quark matter (up, down and strange flavors) in MNJL model [68] which is a combination of NJL Lagrangian and its Fierz transformation by using weighting factors (1 − α) and α.By imposing charge neutrality and chemical equilibrium, we obtained the EOSs of SQM for three different values of α, (α = 0.5, α = 0.8, and α = 0.94) and two different values of Bag constant (B 1/4 = 117M eV and B 1/4 = 130M eV ).We ensured that the EOSs satisfied causality and dynamical stability conditions.We then presented a brief introduction on the modified TOV equation in massive gravity.We explained the parameters appeared in TOV equation in massive gravity, including C 1 and C 2 .We also derived the differential equation for the metric function H(r) in massive gravity to calculate the tidal deformability.
To assess the structural properties of SQS, we initially considered B 1/4 = 117M eV in the EOS and then by setting a fixed value for C 1 (the mass-radius relation is independent of C 1 [61]) and selecting different values of C 2 , we calculated the structural properties of SQS.The results analyzed to ensure compliance with the constraint Λ 1.4M⊙ ≲ 580.We observed that increasing the α parameter led to an increase in the maximum gravitational mass but also resulted the aforementioned constraint to be satisfied in a smaller range of massive gravity parameters.Consequently, the allowed mass range of the star decreased.For the case of B 1/4 = 117M eV , we obtained the maximum allowed M T OV approximately 2.13M ⊙ .
Subsequently, we altered, the value of B and considered B 1/4 = 130M eV .This change resulted in a softer EOS and a broader range of massive parameters satis-fying the astronomical constraints.We obtained masses that fall within the mass gap region (M ∼ (2.5 − 5)M ⊙ ) that simultaneously satisfied all the observational constraints.Moreover, we calculated Schwarzschild radius to show that compact objects studied in this paper were not black holes.
In general relativity, soft EOSs yield small masses for quark stars, while stiff EOSs fail to satisfy the the Λ constraint.However, our study showed that massive gravity allows for the existence of quark stars with higher masses than those predicted by general relativity, while still satisfying the constraints imposed by gravitational wave observations.Finally in Appendix B, by considering a constant value of α (α = 0.8), we increased the value of the bag constant to show that our results in massive gravity still cover the mass gap region.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135 [77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451 [38,39].The light blue band denotes the remnant mass of GW190425.For further investigation, we increase the value of B 1/4 to 148M eV .First of all, we have to check the stability condition of the SQM. Figure .24 shows that this condition is well established for all baryon number densities.In the next step, we examine the Λ behavior.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135 [77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451 [38,39].The light blue band denotes the remnant mass of GW190425.

3 )FIG. 1 :
FIG.1: Quark number density (n) versus chemical potential (µ).Notably, n u,d devotes to up and down quark number densities, and ns is related to the strange quark number density.

- 5 FIG. 7 :
FIG.7: Mass-Radius diagram for B 1/4 = 117M eV, α = 0.5, C2 = −0.1, and different values of C1.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].

15 FIG. 9 :
FIG. 9: Mass-Radius diagram for B 1/4 = 117M eV, α = 0.8, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].

15 FIG. 11 :
FIG.11: Mass-Radius diagram for B 1/4 = 117M eV , α = 0.94, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].

Figure. 14
represents the adiabatic index for different values of α and B. As this figure shows, for any choice of α, dynamical stability is increased by increasing B.

14 FIG. 15 :
FIG.15: Mass-Radius diagram for B 1/4 = 130M eV , α = 0.5, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].The light blue band denotes the remnant mass of GW190425.

82 FIG. 17 :
FIG. 17: Mass-Radius diagram for B 1/4 = 130M eV , α = 0.8, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].The light blue band denotes the remnant mass of GW190425.

54 FIG. 19 :
FIG.19: Mass-Radius diagram for B 1/4 = 130M eV , α = 0.94, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].The light blue band denotes the remnant mass of GW190425.

82 FIG. 23 :
FIG.23: Mass-Radius diagram for B 1/4 = 140M eV , α = 0.8, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].The light blue band denotes the remnant mass of GW190425.

82 FIG. 26 :
FIG.26: Mass-Radius diagram for B 1/4 = 148M eV , α = 0.8, C1 = 10 −5 , and different values of C2.The gray and orange regions are the mass -radius constraints from the GW170817 event.The black region shows pulsar J0740+6620, The green and black hatched regions represent pulsars J2215+5135[77] and PSR J095-0607, respectively.The red region amounts to the secondary component of GW190814.The brown and the purple regions show two different reports of the pulsar J0030+0451[38,39].The light blue band denotes the remnant mass of GW190425.
1.4M⊙ ≲ 580 is not established for these values of C 1 .The structural properties of SQS for B 1/4 = 117M eV , α = 0.5, C 2 = −0.1, and different values of C 1 are in Table.II.On the other hand, our results in Table.II indicate that the obtained objects cannot be black holes because the radii of these compact objects are bigger than the Schwarzschild radii.Besides, the strength of gravity increases by increasing the value of |C 2 |.

TABLE V :
Structural properties of SQS for B 1/4 = 130M eV , α = 0.5, C1 = 10 −5 , and different values of C2.The results falling in mass gap have been separated in different column.
1/4= 130M eV and α = 0.5).Moreover, we checked the Λ of M T OV for different values of C 2 .We observed that for |C 2 | ≳ 0.84 the values of Λ M T OV became zero.The results for structural properties are given in Table.VI.On the other hand, our calculation for the compactness and the modified Schwarzschild radius show that these compact objects cannot be black holes because their radii are bigger than the Schwarzschild radii, i.e., R > R Sch (see Table.VI, for more details).Also, the compactness increases by increasing the value of |C 2 |.In other words, the strength of gravity augments when the value of |C 2 | increases (see Table.VI, for more details).

TABLE VI :
Structural properties of SQS for B 1/4 = 130M eV , α = 0.8, C1 = 10 −5 , and different values of C2.The results falling in mass gap have been separated in different column.

TABLE VII :
Structural properties of SQS for B 1/4 = 130M eV , α = 0.94, C1 = 10 −5 , and different values of C2.The results falling in mass gap region have been separated in a different column.