Quark Stars in $ f(T,{\cal{T}})$ Gravity: Structure, Stability, and Observational Constraints

Compact stars are the densest objects that can be directly observed, and therefore serve two purposes: they probe the behaviour of matter at densities inaccessible to laboratories, and they test gravity in the strong-field regime. Over the past decade, the observational situation has improved dramatically. The gravitational-wave event GW170817 provided the first direct measurement of the tidal deformability of neutron-star matter [1, 2]. Radio timing of PSR J0740$ + $6620 established that any acceptable model must support at least two solar masses [3], and NICER X-ray observations provided a simultaneous mass–radius measurement for the same pulsar [4, 5]. At the extremes of the mass distribution are the “black widow” pulsar PSR J0952$ - $0607 with $ M\approx2.35\, M_{\odot} $ [6], the $ 2.6\, M_{\odot} $ secondary of GW190814, whose nature remains debated [7], and the unusually light central compact object in HESS J1731$ - $347 with $ M\approx0.77\, M_{\odot} $ [8]. Together, these observations define a set of targets that any theory of dense matter, or of gravity, must address.

One long-standing possibility for the interiors of such objects is quark matter. If strange quark matter is the true ground state of hadronic matter, as conjectured by Witten and by Farhi and Jaffe [9, 10], then self-bound quark stars can exist [11, 12]. The simplest quantitative description of such stars is the MIT bag model [13], in which confinement is represented by a constant vacuum pressure B. More refined treatments are based on perturbative QCD at high density, a programme that goes back to the pioneering calculations of Freedman and McLerran, Baluni, and Toimela [14−17] and has now reached partial N$ ^3 $LO accuracy [18−23]. Interpolations between the perturbative regime and stellar densities yield quark-matter equations of state that can be confronted directly with pulsar and gravitational-wave data [24−27]. At the densities relevant to stellar cores, quark matter is also expected to pair, most symmetrically in the colour-flavour-locked phase [28−32], and quark or hybrid configurations have been proposed as explanations of both GW190814 and HESS J1731$ - $347 [33−36]. Quark stars have likewise been used as test beds for modified gravity, for instance in 4D Einstein–Gauss–Bonnet gravity [37] and in Rastall–Rainbow gravity [38]. Related programmes have treated quark stars in the curvature-based counterpart $ f(R,T) $ gravity [39] and in its matter-Lagrangian extension $ f(R,L_m,T) $ gravity [40]; the present work extends this line of trace-coupled constructions to the torsional sector. The trace-coupled curvature theory has now been applied to a broad range of stellar matter, from strange and compact configurations [41, 42] to density-dependent bag descriptions [43] and realistic nuclear equations of state, including magnetized neutron-star matter [44] and relativistic mean-field models confronted with GW170817, pulsar timing, and NICER data [45].

On the gravity side, teleparallel theories provide an attractive alternative starting point. In the teleparallel equivalent of general relativity, the dynamical variable is the tetrad and gravity is encoded in torsion rather than curvature; promoting the torsion scalar T to an arbitrary function $ f(T) $ gives a modified theory whose field equations remain second order [46−49]. Relativistic stars are known to exist in $ f(T) $ gravity [50], although the choice of tetrad requires care: for nonlinear $ f(T) $, a diagonal tetrad in spherical coordinates generates spurious constraints, and one must instead use a “good” (rotated) tetrad or, equivalently, the covariant formulation with a nontrivial spin connection [51, 52]. A further extension couples the geometry directly to matter. On the curvature side, this leads to $ f(R, {\cal{T}}) $ gravity [53]; its torsional counterpart, $ f(T, {\cal{T}}) $ gravity, was introduced by Harko, Lobo, Otálora, and Saridakis [54], where $ {\cal{T}} $ is the trace of the energy–momentum tensor. The trace coupling makes the theory nonconservative, $ \nabla_\mu\Theta^{\mu}{}_{\nu}\neq0 $, which opens qualitatively new phenomenology. The $ f(T, {\cal{T}}) $ framework has since been studied in many settings: late-time cosmology and dynamical-systems analyses [55−59], gravitational baryogenesis and Big Bang nucleosynthesis constraints [60, 61], thick branes [62], black-hole thermodynamics [63], wormholes [64−66], gravastars [67], charged fluids with gravitational decoupling [68], and anisotropic strange-star models based on prescribed metric ansatze [69−71].

For isotropic compact stars obtained by direct integration of the stellar structure equations, however, the available results remain limited. Pace and Levi Said derived a working TOV model for quark stars in $ f(T, {\cal{T}}) $ gravity and solved it with the MIT bag model [72]; they later treated neutron stars in a perturbative scheme [73]. Their quark-star study, however, retained a nonzero cosmological constant, adopted the strange-quark-mass-corrected bag parameter $ \omega=0.28 $, and explored only a few negative values of the coupling through single-star mass profiles. A nonperturbative treatment with the conformal bag equation of state is still missing, as are a systematic maximum-mass scan over the full admissible range of the coupling, including its positive branch, and a direct confrontation with the modern observational bounds listed above. Since a positive β corresponds to a negative coupling in the curvature-based description, as we show in Sec. II C, this is also the region that the $ f(R,{\cal{T}}) $ stellar literature has largely left aside. This gap is addressed here.

In this paper, we derive the modified Tolman–Oppenheimer–Volkoff (TOV) equations for the linear model $ f(T, {\cal{T}})=T+\beta {\cal{T}} $ with matter Lagrangian $ {\cal{L}}_{m}=p $, solve them numerically for the conformal MIT bag equation of state, and scan the coupling over the full admissible range $ \beta\in[-10,12.25] $. Our field equations were cross-checked with independent computer-algebra software. The central result is that the maximum gravitational mass $ M_{\max}(\beta) $ is non-monotonic: it peaks at $ 2.021\, M_{\odot} $ near $ \beta\approx3.10 $ and collapses as $ \beta\to4\pi $, where the pressure equation develops a singular denominator. The $ 2\, M_{\odot} $ pulsar constraint then selects the window $ \beta\in(0,\sim 5.5) $.

The paper is organised as follows. Section II sets up the action and field equations. Section III derives the modified TOV system and states the equation of state. Section IV collects the stability criteria, and Sec. V describes the numerical method and its validation. The results are presented in Sec. VI and discussed in Sec. VII. Throughout, we use geometrised units $ G=c=1 $ and the metric signature $ (-,+,+,+) $.

Teleparallel gravity and its generalisations are formulated in Weitzenböck geometry, in which the tetrad $ e^{A}{}_{\mu} $ is the dynamical variable and the metric is a derived quantity,

where $ \eta_{AB}=\mathrm{diag}(-1,1,1,1) $. Greek indices $ \mu,\nu,\dots $ denote spacetime coordinate labels, whereas capital Latin indices $ A,B,\dots $ denote tangent-space labels. The Weitzenböck connection is curvature-free; hence, all gravitational information is encoded in the torsion tensor

where $ \omega^{A}{}_{B\mu} $ denotes the inertial spin connection. The contorsion tensor and superpotential are defined as

and the torsion scalar is

When the teleparallel Lagrangian is generalised to an arbitrary function $ f(T, {\cal{T}}) $ of the torsion scalar and the trace $ {\cal{T}}=\Theta^{\mu}{}_{\mu} $ of the energy–momentum tensor [54], the total action is

We use geometrised units, $ G=c=1 $. The matter energy–momentum tensor is defined by varying the matter action with respect to the tetrad,

For the stellar interior, we assume a perfect fluid,

and adopt the matter Lagrangian

The choice of $ {\cal{L}}_{m} $ is not merely cosmetic: it determines the variational response of the energy–momentum tensor that enters $ \delta {\cal{T}} $ and, consequently, the form of the trace-coupling term in the field equations.

We vary Eq. (6) with respect to the tetrad $ e^{A}{}_{\mu} $. The elementary variations are

The variation of the torsion scalar, integrated by parts using $ \partial_{\mu}f_{T}=f_{TT}\partial_{\mu}T +f_{T {\cal{T}}}\partial_{\mu} {\cal{T}} $, generates the purely gravitational terms

Writing $ {\cal{T}} = g^{\mu\nu}\Theta_{\mu\nu} $ and defining

Adopting $ {\cal{L}}_{m}=p $ and assuming $ \partial^{2}p/\partial g^{\mu\nu}\partial g^{\alpha\beta}=0 $, we find that the trace-coupling contribution yields the matter-coupling term

Setting the coefficient of $ \delta e^{i}{}_{\nu} $ to zero yields the $ f(T, {\cal{T}}) $ field equations,

We work in the pure-tetrad formulation, with $ \omega^{A}{}_{B\mu}=0 $. For the linear model $ f=T+\beta{\cal{T}} $ ($ \alpha=1 $ by normalisation, and $ \varphi=0 $ because it is dynamically negligible on stellar scales), consistency of the stellar equations requires

the latter reduces to $ \beta\neq4\pi $ for $ \omega=1/3 $.

The torsion scalar and the Ricci scalar differ by a boundary term, $ T=-R-B $, with $ B=(2/e)\,\partial_{\mu}(e\,T^{\nu\mu}{}_{\nu}) $. For a gravitational sector linear in T, this boundary contribution does not affect the field equations. Thus, the teleparallel field equations reproduce those of general relativity, and the two descriptions differ only in the variables in which they are written [74−76]. It is therefore natural to ask what the trace-coupled model $ f(T,{\cal{T}})=T+\beta{\cal{T}} $ adds to $ f(R,{\cal{T}})=R+2\lambda{\cal{T}} $ [53]. We address this question explicitly here.

For a perfect fluid with $ {\cal{L}}_{m}=p $, the curvature-based theory gives the effective source

and, upon imposing $ dp/d\rho=\omega $ in the contracted Bianchi identity, the hydrostatic coefficient

A comparison of Eqs. (16) and (17) with Eqs. (21)–(26) shows that the two systems coincide term by term under

The pressure-gradient coefficients satisfy $ {\cal{K}}(-2\lambda)= {\cal{K}}_{R}(\lambda) $, and the radial pressures entering Eq. (21) are identical. Moreover, the algebraic closure (24) reproduces $ e^{-B_m}=1-2m_{\rm eff}/r $ with $ m_{\rm eff}=\int 4\pi r^{2}\rho_{\rm eff}\,dr $, with both the term linear in M and the term proportional to the bag constant matching identically. The singular surfaces also map onto one another: $ \beta=4\pi $ corresponds to $ \lambda=-2\pi $, precisely where $ {\cal{K}}_{R} $ diverges.

Two consequences follow, setting the scope of the rest of this paper. First, within the linear model with $ {\cal{L}}_{m}=p $, the physical content of our results is determined by the parameter range rather than by the choice of geometric variable, and the numerical values reported below may be interpreted in either language. Second, the map reverses orientation: the positive-β region, where the maximum mass is raised and the two-solar-mass window opens, corresponds to $ \lambda<0 $. Existing $ f(R,{\cal{T}}) $ stellar studies have concentrated on $ \lambda>0 $ [41, 42, 77], and the earlier $ f(T,{\cal{T}}) $ quark-star work [72] considered only $ \beta<0 $, so the branch examined here has not been explored from either side. The correspondence is also specific to the present case. A nonlinear $ f(T,{\cal{T}}) $ leaves the boundary term in the field equations, and the two families separate [78]; the same separation occurs once a matter Lagrangian other than $ {\cal{L}}_{m}=p $ is adopted [79].

We assume a static, spherically symmetric spacetime,

where the subscript m in $ B_m $ distinguishes this metric function from the bag constant B. In the pure-tetrad formulation with a vanishing inertial spin connection $ \omega^{A}{}_{B\mu}=0 $, we employ the rotated good tetrad compatible with the metric (19) [50, 52]. For this tetrad choice, the torsion scalar is given by

where a prime denotes differentiation with respect to r. Using the rotated good tetrad ensures the consistency of the pure-tetrad field equations and avoids the spurious off-diagonal constraints that arise from the diagonal tetrad in nonlinear $ f(T) $ models [51].

Projecting Eq. (14) onto the metric (19) for $ f=T+\beta {\cal{T}} $, with a perfect-fluid source and $ {\cal{L}}_{m}=p $, we obtain

with the algebraic closure

where

and $ B_{\rm km} $ is the bag constant in geometrized units. Equation (22) follows from the contracted Bianchi identity applied to the effective source. Because the trace coupling renders the theory non-conservative for $ \beta\neq0 $, the hydrostatic equilibrium equation differs from its general-relativistic form and reduces to it only in the limit $ \beta\to0 $. In this limit, Eqs. (21)–(24) reduce to the standard Tolman–Oppenheimer–Volkoff system [80, 81].

For the conformal equation of state, $ \omega=1/3 $, the denominator of Eq. (22) simplifies to

Thus, the only singular value associated with the hydrostatic equation is $ \beta=4\pi\approx12.57 $; no regular solution exists at this value, and it is therefore excluded. The factor $ (16\pi-7\beta) $, which appears in the intermediate form of the conservation law before the equation of state is imposed, does not, by itself, produce a singularity or a sign change in the full denominator.

We model the quark-star interior using the conformal MIT bag model [10, 13],

with $ B=60\,\mathrm{MeV\,fm^{-3}} $. The stellar surface is defined by the condition $ p(R)=0 $. Outside the star, the matter sector vanishes; hence $ {\cal{T}}=0 $ and the model reduces to TEGR. The exterior geometry is therefore exactly Schwarzschild, $ e^{-B_m}|_{r>R}=1-2M/r $. Following the operational prescription adopted by Pace and Levi Said [72], we identify the terminal value $ M(R) $ of the integrated mass function (23) with the mass parameter of this exterior solution and use it to construct the mass–radius sequences.

Throughout this work, we adopt the conformal MIT bag equation of state with $ \omega=1/3 $. This choice enters the hydrostatic equation (22) and the closure (24) directly; holding it fixed therefore isolates the stellar response to β from its response to the microphysics, as appropriate for an initial scan of the coupling. The trade-off is that a self-bound conformal bag star cannot reach the masses of the heaviest observed pulsars, a limitation to which we return in Sec. VII. Stiffer descriptions are available and have been used in trace-coupled curvature theory, including a density-dependent bag parameter [43], relativistic mean-field matter [44, 45], and the unified interacting quark equation of state with its perturbative-QCD correction [82, 83]. Each of these raises the maximum mass already in general relativity and would shift the window quoted below; we consider this separately rather than incorporating it into the present scan.

Within the one-parameter equilibrium sequence, the onset of instability is identified by the first maximum of $ M(\rho_c) $ [84]. Configurations with $ dM/d\rho_c>0 $ define the candidate stable branch identified by the turning-point criterion, whereas those with $ dM/d\rho_c<0 $ beyond $ M_{\max} $ are unstable. The full radial perturbation equations in $ f(T,{\cal{T}}) $ gravity have not yet been derived. Moreover, in a non-conservative matter–geometry coupling, neither the pulsation equations nor the applicability of the general-relativistic turning-point theorem can be assumed to carry over unchanged. Therefore, the criterion is used here as a diagnostic rather than as a proof of dynamical stability, and a complete Sturm–Liouville pulsation analysis is reserved for future work. In the curvature-based theory, the corresponding calculation has been carried out, and the radial spectrum of $ R+2\lambda{\cal{T}} $ stars was found to be consistent with the turning point of the equilibrium sequence [77, 85]. Through the correspondence of Sec. II C, that result covers the linear model treated here, providing a reason to expect the diagnostic to be reliable, although it is not a substitute for deriving the pulsation equations in the torsional variables.

The conformal bag EOS yields

which is constant, independent of β, and satisfies both causality ($ c_s^2<1 $) and the Le Chatelier condition ($ c_s^2>0 $) for all $ \rho>4B $. This independence from β reflects the fact that $ c_s^2 $ is determined entirely by the equation of state, which does not depend on the gravitational coupling. The β-dependence of dynamical radial stability is a separate issue and is addressed here only diagnostically through the turning-point criterion discussed above.

The effective adiabatic index [86] is

For all physically admissible configurations ($ \rho>4B $),

which holds identically. Since $ d\Gamma/d\rho=-4B/(\rho-4B)^{2}<0 $, the central value $ \Gamma_c=\Gamma(\rho_c) $ is the minimum of Γ within the stellar interior. Thus, $ \Gamma_c>4/3 $ implies that this condition holds throughout the star for all admissible β. This provides a local matter-stability diagnostic; however, it does not, by itself, establish the absence of exponentially growing radial modes in the trace-coupled theory.

We integrate Eqs. (21)–(24) from $ r_0=10^{-4}\,\mathrm{km} $, with initial conditions $ M(r_0)=\dfrac{4\pi}{3}r_0^{3}\rho_c $ and $ p(r_0)=p(\rho_c) $, using the SCIPY RK45 integrator with relative and absolute tolerances of $ 10^{-10} $ and $ 10^{-12} $, respectively. The surface $ p(R)=0 $ is located using an ODE event. Maximum-mass configurations are identified using the turning-point method: a sign change of $ dM/d\rho_c $ from positive to negative over a 300-point density grid $\rho_c\in[200, $ $ 3000]\,\mathrm{MeV\,fm^{-3}} $ identifies the coarse turning point, which is then refined on a 400-point local grid over $ \pm100\,\mathrm{MeV\,fm^{-3}} $.

We apply a six-level validation protocol before conducting any production run.

1. GR-limit gate: At $ \beta=0 $, the solver reproduces an independent GR MIT bag solution with $\Delta M= $ $ 1.816\times10^{-13}\, M_{\odot} $.

2. Tolerance convergence: The output remains stable under a two-order-of-magnitude tightening of the ODE tolerances.

3. Density-grid convergence: the results remain stable when the central-density grid is doubled to 600 points.

4. GR benchmark: The GR maximum mass $ M_{\max}^{\rm GR}= 1.964\, M_{\odot} $ at $ B=60\,\mathrm{MeV\,fm^{-3}} $ agrees to within $ 1 $% with published MIT bag model results [11, 12].

5. Admissibility guard: The singular surface $ \beta=4\pi $ is identified and flagged before integration.

6. Cross-check: An independent evolution of $ \rho(r) $ agrees with the primary $ p(r) $ evolution to within $ \sim10^{-15} $.

The reduced stellar equations were cross-checked using independent computer-algebra software, and the numerical implementation was validated against an independent TOV integration.

The mass–radius sequences for $ \beta\in\{-10,-5,0,3,5,10\} $ are displayed in Fig. 1; the corresponding M–$ \rho_c/B $ sequences are shown in Fig. 2. Filled circles indicate the maximum-mass (turning-point) configuration on each candidate stable branch, whereas faded lines denote the unstable branch ($ dM/d\rho_c<0 $).

Figure 1.  Mass–radius relation for quark stars in $ f(T,{\cal{T}})= $$ T+\beta{\cal{T}} $ gravity using the conformal MIT bag EOS ($ B=60\,\mathrm{MeV\,fm^{-3}} $, $ \omega=1/3 $) for $ \beta\in\{-10,-5,0,3,5,10\} $. This set includes $ \beta=3 $, corresponding to the location of the $ M_{\max} $ peak identified in Fig. 5, so that the peak configuration is directly visible in the mass–radius plane. Solid lines indicate the candidate stable branch identified by the turning-point criterion ($ dM/d\rho_c>0 $), whereas faded lines indicate the unstable branch. Filled circles mark the maximum-mass configuration. The dashed horizontal line indicates the $ 2\, M_{\odot} $ observational threshold. The shaded contours show the GW170817 posterior [2] derived under hadronic-EOS priors; their overlap with quark-star sequences is shown for orientation only.

Figure 2.  Gravitational mass M as a function of the dimensionless central density $ \rho_c/B $ for the same parameter set as in Fig. 1. Solid lines indicate the candidate stable branch identified using the turning-point criterion; faded lines indicate the unstable branch beyond the turning point (filled circle).

The maximum mass varies non-monotonically with β. The asymmetry between the negative and positive branches is examined in Sec. VI C and Sec. VII. For $ \beta<0 $, the coupling suppresses $ M_{\max} $ below the GR value: $ M_{\max}= 1.549\, M_{\odot} $ at $ \beta=-10 $ and $ 1.765\, M_{\odot} $ at $ \beta=-5 $ (Table 1). The GR limit ($ \beta=0 $) yields $ M_{\max}=1.964\, M_{\odot} $ at $ R=10.71\,\mathrm{km} $, consistent with standard MIT bag results [11]. At $ \beta=3 $, the maximum mass increases to $ 2.021\, M_{\odot} $, exceeding the $ 2\, M_{\odot} $ threshold [3]. At $ \beta=5 $, it decreases to $ 1.985\, M_{\odot} $, and at $ \beta=10 $, it falls sharply to $ 1.162\, M_{\odot} $, reflecting the approach to the singular surface at $ \beta=4\pi $.

The central density at $ M_{\max} $ increases with β: $ \rho_c^{\max}/B=16.03 $ at $ \beta=-10 $, $ 19.25 $ at $ \beta=0 $, $ 21.25 $ at $ \beta=3 $, and $ 36.19 $ at $ \beta=10 $ (Table 1).

The compactness $ {\cal{C}}=M/R $ (in geometrized units, $ G=c=1 $) along the candidate stable branch is shown in Fig. 3. For all configurations, the compactness remains below the standard GR Buchdahl benchmark, $ {\cal{C}}=4/9\approx0.444 $ [87]; the peak compactness is $ {\cal{C}}=0.279 $ at $ \beta=3 $.

Figure 3.  Compactness $ {\cal{C}}=M/R $ (in geometrized units) as a function of gravitational mass M along the candidate stable branch identified by the turning-point criterion for each β. Filled circles indicate the maximum-mass configuration. All configurations remain below the standard Buchdahl compactness benchmark in GR, $ {\cal{C}}=4/9\approx0.444 $.

The surface gravitational redshift,

is shown in Fig. 4. The calculated surface redshifts also remain below the commonly quoted GR reference value of $ z_s=0.85 $ [88]. At $ M_{\max} $, the redshift ranges from $ z_s=0.354 $ ($ \beta=-10 $) to $ z_s=0.503 $ ($ \beta=3 $); the complete set of values is provided in Table 1. The enhancement of $ \Delta z_s=0.026 $ at $ \beta=3 $ relative to GR could, in principle, be observed through spectroscopic measurements of surface absorption lines [89].

Figure 4.  Surface gravitational redshift $z_s=(1-2\mathcal{C})^{-1/2}-1$ as a function of M along the candidate stable branch identified by the turning-point criterion. All configurations remain below the commonly quoted GR reference value $z_s=0.85$ [88]. See Table 1 for values at $M_{\max}$.

Figure 5 shows the results of a dense scan over $ \beta\in[-10,12.25] $ with step size $ \Delta\beta=0.25 $. The curve rises from $ M_{\max}=1.549\, M_{\odot} $ at $ \beta=-10 $, crosses the GR reference near $ \beta=0 $, reaches a maximum of $ M_{\max}=2.021\, M_{\odot} $ at $ \beta\approx3.10 $, and then decreases steeply, falling below $ 2\, M_{\odot} $ near $ \beta\approx5.5 $ and approaching zero as $ \beta\to4\pi^{-} $.

Figure 5.  Maximum gravitational mass $ M_{\max} $ as a function of the trace-coupling parameter β for quark stars in $ f(T,{\cal{T}})=T+\beta{\cal{T}} $ gravity ($ B=60\,\mathrm{MeV\,fm^{-3}} $, $ \omega=1/3 $). Each point is obtained from a refined turning-point search. The green dashed line indicates the GR reference value ($ M_{\max}^{\rm GR}=1.964\, M_{\odot} $), and the black dotted line denotes the $ 2\, M_{\odot} $ threshold. The grey dashed vertical line marks the singular surface $ \beta=4\pi\approx12.57 $. The red star indicates the peak at $ \beta\approx3.10 $, $ M_{\max}=2.021\, M_{\odot} $, which is the central result of this work.

The non-monotonic profile cannot be attributed to a sign change at $ \beta=16\pi/7 $. As shown in Eq. (27), for $ \omega=1/3 $ the denominator of Eq. (22) reduces to $ \dfrac{4}{3}(4\pi-\beta) $, whose only zero is at $ \beta=4\pi $.

reduces to the general-relativistic value $ -1/2 $ at $ \beta=0 $ and increases steadily in magnitude for $ \beta>0 $, reaching magnitudes that are $ 1.16 $ and $ 2.95 $ times the GR value at $ \beta=3 $ and $ \beta=10 $, respectively. Thus, the coupling does not enhance the maximum mass by slowing the pressure fall-off. Instead, the non-monotonicity emerges from the combined influence of the trace coupling on the pressure gradient (22), the metric potential (21), and the algebraic radial-metric closure (24). These competing effects are resolved only through the numerical solution of the full coupled system.

The trade-off among these effects can nevertheless be inferred directly from the coefficients. At $ \omega=1/3 $, the closure (24) contains $ \Sigma=-2+\beta/6\pi $ and $ \xi=\dfrac{10}{3}\beta B_{\rm km} $. A positive β therefore has two simultaneous effects: it reduces $ |\Sigma| $ below its general-relativistic value of $ 2 $, weakening the mass term in $ e^{-B_{m}} $, and adds a positive ξ contribution proportional to the bag constant. Both effects increase $ e^{-B_{m}} $. Consequently, the interior metric becomes shallower than in general relativity, and a given central density supports more mass before the turning point is reached. Opposing this behavior is $ {\cal{K}}(\beta) $, which steepens the pressure fall-off and reduces the stellar radius. For moderate positive β, the two closure terms dominate, which is why $ M_{\max} $ rises to $ 2.021\, M_{\odot} $ near $ \beta\approx3.10 $ while the radius remains close to its general-relativistic value (Table 1). Beyond that point, the balance shifts. The denominator $ 4\pi-\beta $ in (22) decreases toward zero, $ |{\cal{K}}| $ grows without bound, and the pressure drops so rapidly with radius that the surface is reached before much mass has accumulated. The star contracts, $ \rho_c^{\max}/B $ increases from $ 21.25 $ at $ \beta=3 $ to $ 36.19 $ at $ \beta=10 $, and $ M_{\max} $ collapses toward zero as $ \beta\to4\pi^{-} $. The turnover is therefore the intersection of a closure effect that saturates and a pressure-gradient effect that diverges.

We have studied quark stars in $ f(T,{\cal{T}})=T+\beta{\cal{T}} $ gravity with the conformal MIT bag EOS at $ B=60\,\mathrm{MeV\,fm^{-3}} $. The maximum mass $ M_{\max}(\beta) $ is non-monotonic in β: it peaks at $ 2.021\, M_{\odot} $ for $ \beta\approx3.10 $, exceeding the GR value by $ \Delta M=+0.057\, M_{\odot} $ ($ +2.9 $%), and is suppressed as the singular surface $ \beta=4\pi $ is approached. The conventional $ 2\, M_{\odot} $ pulsar requirement [3] is satisfied for $ \beta\in(0,\sim 5.5) $; thus, a moderate positive trace coupling is compatible with this observational threshold within a finite interval. However, the peak mass $ M_{\max}=2.021\, M_{\odot} $ does not reach the central mass estimate of PSR J0952$ - $0607, $ M\simeq2.35\, M_{\odot} $ [6]. Accommodating such a high mass would require either a stiffer quark-matter equation of state, additional microphysical effects, or an extension of the present one-parameter model; we identify this as a direction for future work rather than a result established here. All configurations on the candidate stable branch are causal, satisfy $ \Gamma>4/3 $ analytically, and remain below the standard GR compactness and surface-redshift benchmarks. The non-monotonic behaviour cannot be attributed to a sign change of the factor $ (16\pi-7\beta) $: for $ \omega=1/3 $ the denominator of Eq. (22) reduces to $ \dfrac{4}{3}(4\pi-\beta) $, and the pressure-gradient coefficient (33) increases monotonically in magnitude for $ \beta>0 $. Instead, the turnover reflects the competition among the coupling's effects on the pressure gradient, the metric potential, and the algebraic closure, and is established here numerically rather than analytically.

The correspondence established in Sec. II C places these numbers in a broader context. Under $ \beta=-2\lambda $, the peak occurs at $ \lambda\approx-1.55 $, the observationally allowed window becomes $ \lambda\in(-2.75,0) $, and the singular surface is $ \lambda=-2\pi $. Interpreted in this way, the results constitute a statement about the negative-coupling branch of $ R+2\lambda{\cal{T}} $ gravity as well as about its torsional counterpart, and this branch has not previously been mapped with a self-bound equation of state. The teleparallel formulation provides the setting in which the positive, physically favoured side of the coupling is the natural one to explore.

Several limitations should be noted. The full radial pulsation equations in $ f(T,{\cal{T}}) $ gravity have not been derived; our stability statements rest on the turning-point criterion [84] and on the analytic bound $ \Gamma>4/3 $ [86], and a complete Sturm–Liouville analysis would provide a more robust basis for them. The tidal Love number $ k_2 $ and the dimensionless deformability Λ, which are directly constrained by GW170817 [1, 2], remain to be computed in this framework. As in the original compact-star treatment of Pace and Levi Said [72], the terminal value of the integrated mass function has been identified with the gravitational mass of the star. An explicit analysis of the surface junction conditions would provide a useful independent test of this prescription; such an analysis is left for future work. Finally, more realistic quark-matter descriptions, including perturbative-QCD corrections [24−26] and colour-superconducting phases [28, 30, 31], together with a scan over the bag constant B, would broaden the parameter space and sharpen the comparison with objects such as HESS J1731$ - $347 and the GW190814 secondary [7, 8]. The unified interacting equation of state, with its repulsive perturbative-QCD term, is the most promising of these for reaching the heaviest pulsars [82, 83], and a density-dependent bag parameter offers a further route [43]. Whether the coupling window quoted here survives once the microphysics is stiffened remains an open question.

In summary, quark stars in $ f(T,{\cal{T}})=T+\beta{\cal{T}} $ gravity with the conformal MIT bag EOS display a non-monotonic maximum-mass profile bounded above by the singular surface at $ \beta=4\pi $, with a peak of $ M_{\max}=2.021\, M_{\odot} $ at $ \beta\approx3.10 $ and observational viability for $ \beta\in(0,\sim 5.5) $. These results establish $ f(T,{\cal{T}}) $ gravity as a viable framework for massive quark stars, equivalently the negative-λ branch of $ R+2\lambda{\cal{T}} $ gravity, and motivate further studies of tidal deformability, radial stability, and extended EOS families.