Energy–momentum–squared gravity: charged quark star solutions with unified interacting matter

Compact astrophysical objects—white dwarfs, neutron stars, and black holes—occupy a central place in modern astrophysics as natural laboratories for ultra-dense matter and strong-gravity phenomena. Their observable properties reveal how matter behaves at supranuclear densities and provide tests of gravitational physics in regimes that cannot be reproduced in the laboratory. Observations across the electromagnetic spectrum and via gravitational waves have therefore made the study of compact remnants a key interface among nuclear physics, quantum chromodynamics (QCD), and relativistic gravity.

Within this broader landscape, neutron stars mark the transition from degenerate electron matter to degenerate baryonic matter. However, theoretical considerations strongly suggest that, at sufficiently high densities, hadrons may dissolve into deconfined quarks. This possibility motivates the concept of quark stars or strange stars, compact configurations in which up, down, and strange quarks provide the dominant pressure support. Early work by Bodmer on "collapsed nuclei'' and by Witten on the cosmic separation of phases framed strange quark matter as a possible absolute ground state of strongly interacting matter, laying the foundation for modern quark-star scenarios [1–3]. Recent multimessenger constraints, together with QCD-informed modelling of dense matter, have provided evidence that massive neutron stars may harbour quark-matter cores, strengthening the case for deconfinement at the highest central densities [4–6]. Detailed studies of medium effects in strange quark matter and the role of interactions in strange-star configurations further illustrate how microphysical input can substantially alter compact-star phenomenology [7, 8].

To capture these effects in a controlled yet flexible manner, several interaction-based frameworks have been developed to describe quark matter at high density. A particularly versatile approach is the unified interacting quark-matter model proposed by Zhang and Mann, which incorporates perturbative QCD corrections, colour superconductivity, and finite-strange-quark-mass contributions within a thermodynamically consistent parametrization that smoothly interpolates between different quark phases [9]. This unified interacting matter (UIM) framework has been applied to various settings, including anisotropic interacting quark stars, universal relations, and comparisons with recent astrophysical observations [10–13]. Extensions to modified gravity, such as charged quark stars and extreme compact objects in regularized four-dimensional Einstein–Gauss–Bonnet gravity, further demonstrate the utility of interaction-based equations of state for exploring the parameter space of relativistic compact-star models [14]. In the present work, the UIM description provides the microphysical foundation for constructing charged quark-star configurations in energy–momentum–squared gravity.

In parallel with these microphysical developments, a broad spectrum of modified gravity theories has been proposed to address possible departures from General Relativity (GR) in regimes where curvature and matter densities become extreme. These include $ f(R) $ gravity [15, 16], $ f(R,T) $ gravity [17, 18], and symmetric teleparallel generalizations [19–22]. Within this broader landscape, Energy–Momentum–Squared Gravity (EMSG) provides a particularly economical and well-motivated framework in which quadratic contractions of the energy–momentum tensor, $ T_{\mu\nu}T^{\mu\nu} $, supplement the Einstein–Hilbert action and introduce density- and pressure-dependent corrections. These corrections become most relevant in the ultra-dense interiors of neutron and quark stars [23, 24]. In addition to modifying the effective gravitational coupling, EMSG generally relaxes the covariant conservation of the matter energy–momentum tensor. This gives rise to an extra force associated with nonminimal matter–geometry interaction, which may leave observable signatures across astrophysical and cosmological settings [24−26].

Of particular relevance to the present work are applications of EMSG to stellar structure. Neutron-star analyses have placed direct bounds on the EMSG coupling using realistic hadronic equations of state, demonstrating that the quadratic matter term can appreciably shift mass–radius relations and maximum masses relative to GR while remaining consistent with observed pulsar properties [27, 28]. Subsequent studies have extended these investigations to isotropic and anisotropic quark stars, colour-flavour-locked configurations, and interacting quark matter, showing that EMSG can support heavier and more compact stars than those predicted by standard gravity [29–34]. Charged quark stars in EMSG have also been explored, indicating that the combined effects of electric fields and nonlinear matter–geometry coupling further enhance the attainable masses and modify stability properties in ways that may be relevant to the heaviest known compact objects [35]. These developments establish EMSG as a well-motivated framework for revisiting the equilibrium and stability of quark stars and provide the theoretical backdrop for the charged, interacting quark-star models considered in this work. More recent studies have further refined the EMSG compact-object picture, including neutron-star structure and curvature diagnostics across hadronic and hadron–quark equations of state [36], universal relations linking tidal deformability, compactness, and oscillation frequencies for (proto-)neutron stars [37], anisotropic stellar interiors constructed through gravitational decoupling [38], and charged solutions in EMSG confronted with Event Horizon Telescope observations [39].

The possibility that compact stars may carry a net electric charge introduces additional physical ingredients that become particularly relevant in ultra-dense regimes. In such systems, charge-separation processes, enhanced electromagnetic pressure, and modified collapse dynamics can all influence the balance between gravity and internal stresses. The general-relativistic framework for hydrostatic equilibrium with charge was established in early analyses of relativistic charged fluid spheres and collapsing configurations. These studies showed that Coulomb repulsion and electromagnetic self-energies can delay or even prevent gravitational collapse in highly compact objects [40–42]. Subsequent studies have derived mass–radius and compactness bounds for various charged matter models, including incompressible and polytropic equations of state, and demonstrated that a modest net charge can relax classical limits such as the Buchdahl bound and give rise to quasiblack-hole-like configurations that remain regular in their interiors [43–49]. When the stellar interior is composed of deconfined quark matter, models of electrically charged strange quark stars indicate that charge separation in the outer layers and the associated electromagnetic pressure can significantly modify equilibrium conditions. These effects lead to higher maximum masses, altered mass–radius relations, and distinct oscillatory behaviour compared with neutral quark-star sequences [50, 51]. These investigations of charged neutron and quark stars therefore suggest that electric charge is a natural additional degree of freedom in the ultra-dense regime. Its inclusion may be essential for a complete description of the heaviest compact objects, especially when considered together with interaction-based quark-matter equations of state and nonlinear gravity effects such as those predicted by EMSG.

A comprehensive framework that combines charged stellar configurations, quark-matter-based compact stars, and modified-gravity theories such as EMSG provides a particularly rich setting for probing strong-field physics. On the microphysical side, unified interacting descriptions of quark matter capture the roles of QCD interactions, colour superconductivity, and the strange-quark mass in determining the stiffness of the equation of state. On the gravitational side, EMSG introduces nonlinear couplings to the energy–momentum tensor that become increasingly important at high densities. At the same time, electric charge supplies an additional source of pressure support and modifies classical compactness bounds [27, 30, 35, 50, 51]. The interplay among these ingredients can yield novel equilibrium sequences that reach the mass and compactness ranges suggested by the most massive observed compact stars, potentially addressing tensions between data and purely hadronic GR models. This unified perspective is the central motivation of the present work, which investigates charged quark stars composed of unified interacting quark matter within the EMSG framework and confronts the resulting configurations with current astrophysical constraints. While previous studies have examined neutral quark stars in EMSG [29, 30, 34] or charged quark stars with simpler equations of state [35], the present work is the first to combine the unified interacting quark-matter equation of state with electric charge and EMSG nonlinear gravity. This combination yields stellar configurations that extend significantly beyond the predictions of general relativity and existing modified gravity models in both maximum mass and compactness.

This work is organized as follows. In Sec. II, we briefly outline the EMSG framework and derive the modified field equations for the Maxwell–EMSG system, together with the corresponding stellar structure equations for static, spherically symmetric configurations. In Sec. III, we introduce the unified interacting quark-matter equation of state and the charge-density prescription used to model electrically charged quark stars. Sec. IV presents the numerical solutions for neutral and charged configurations and discusses their global properties in terms of the EMSG coupling, interaction parameter, and charge fraction. In Sec. V, we examine the static stability criterion, the adiabatic index, and the sound-speed profiles to assess dynamical stability and causality. Finally, Sec. VI summarizes the main results and outlines prospects for future work.

To investigate the equilibrium and structure of charged compact stars within the framework of EMSG, we first outline the gravitational theory and incorporate electromagnetic contributions. The Einstein-Hilbert action is generalized to include an additional term quadratic in the energy-momentum tensor:

where g denotes the determinant of the metric tensor, R is the Ricci scalar curvature, and $ T_{\mu\nu} $ is the energy–momentum tensor associated with the matter Lagrangian density $ \mathcal{L}_m $. The electromagnetic sector enters through its Lagrangian $ \mathcal{L}_e $, while the parameter α quantifies the deviation from general relativity and determines the coupling strength of the EMSG correction. In the present work, the quadratic EMSG invariant $ T_{\mu\nu}T^{\mu\nu} $ is constructed solely from the material perfect-fluid energy–momentum tensor. The electromagnetic field is excluded from the squared term and enters only through the Maxwell Lagrangian $ \mathcal{L}_e $ and its stress tensor $ \mathcal{E}_{\mu\nu} $. In particular, $ T_{\mu\nu}T^{\mu\nu} $ contains neither an $ \mathcal{E}_{\mu\nu}\mathcal{E}^{\mu\nu} $ contribution nor mixed matter–electromagnetic contractions.

Variation of this action with respect to the metric yields the modified field equations:

where $ G_{\mu\nu} $ denotes the Einstein tensor, $ \mathcal{E}_{\mu\nu} $ is the electromagnetic energy–momentum tensor, and $ \Theta_{\mu\nu} $ is a new tensor arising from the quadratic matter coupling:

where T denotes the trace of $ T_{\mu\nu} $. The matter energy–momentum tensor is related to the Lagrangian density as follows:

For a charged stellar configuration, we describe matter as a perfect fluid and include the electromagnetic contribution, following the treatment of Ray et al. [42]. The corresponding energy-momentum tensors are:

and

where ρ represents the energy density, p the pressure, $ u^\mu $ the fluid four-velocity, and $ F_{\mu\nu} $ the electromagnetic field-strength tensor. This tensor is constructed from the four-potential $ A_\mu $ as $ F_{\mu\nu} = \nabla_\mu A_\nu - \nabla_\nu A_\mu $, where $ \nabla_\mu $ denotes covariant differentiation. The Maxwell equations governing the electromagnetic field are:

where $ j^\mu = \rho_{\rm ch}u^\mu $ denotes the four-current density, and $ \rho_{\rm ch} $ is the electric charge density.

To derive the modified stellar structure equations, we employ the static, spherically symmetric metric:

where the metric functions ψ and λ depend only on the radial coordinate r. This geometry gives $\sqrt{-g}= {\rm e}^{\psi+ \lambda}r^2\sin\theta$ and $u^\mu = {\rm e}^{-\psi}\delta_0^\mu$. For a purely electric field configuration, the only nonzero component of the field strength tensor is $ F^{01} = -F^{10} $, which, through the Maxwell equation (7), leads to:

where the radial charge function is given by

Since $ \nabla^\mu G_{\mu\nu}= 0 $ identically, taking the covariant divergence of the field equations (2) demonstrates the non-conservation of the matter energy–momentum tensor in EMSG:

which can be expressed, using Eqs. (5) and (6), as:

The explicit form of $ \Theta_{\mu\nu} $ in Eq. (3) depends on the choice of the matter Lagrangian density $ \mathcal{L}_m $. Following established approaches in the literature [52, 53], we adopt $ \mathcal{L}_m = p $ in this analysis, which yields $ \Theta_{\mu\nu} = -(\rho^2+ 4\rho p + 3p^2)u_\mu u_\nu $. Although $ \mathcal{L}_m = -\rho $ is an equally admissible choice in GR, $ \mathcal{L}_m = p $ has been shown to be the physically appropriate prescription for perfect fluids in theories with non-minimal matter-geometry coupling [52, 53]. This choice is also consistently adopted throughout the EMSG literature [27, 28, 30, 35]. Substituting this expression into Eqs. (2) and (13), the field equations and energy-momentum non-conservation take the following forms:

Evaluating these equations in the spherically symmetric geometry defined by Eq. (9), the $ 00 $ and $ 11 $ components yield:

where primes denote derivatives with respect to the radial coordinate. In the absence of electric charge, these expressions reduce to those derived in Ref. [27]. The radial component of the non-conservation equation (15) yields the pressure gradient:

To express these structural relations in a more transparent form, reminiscent of the standard Tolman–Oppenheimer–Volkoff (TOV) formulation, we introduce a gravitational mass function that quantifies the total mass–energy enclosed within a radius r. Integrating Eq. (16) gives:

which can be expressed compactly as:

where the mass function is given by:

This decomposition shows that the effective gravitational mass arises from three distinct contributions: the conventional matter energy density, represented by the first integral; the electromagnetic field energy, encoded in the charge distribution; and the modifications induced by the quadratic matter coupling in EMSG, represented by the third term. Setting $ \alpha= 0 $ recovers the standard charged perfect-fluid mass function in general relativity [45, 50], while the neutral limit ($ q=0 $) reproduces the uncharged EMSG result of Ref. [27].

Using the relation (20), the radial field equation (17) determines the gradient of the metric potential:

Combining Eqs. (11), (18), (21), and (22), we obtain the complete system of modified Tolman-Oppenheimer-Volkoff equations governing charged stellar configurations in EMSG [35]:

These differential equations describe hydrostatic equilibrium in the presence of both electromagnetic fields and EMSG modifications. For $ \alpha =0 $, the standard general relativistic TOV equations are recovered, whereas the uncharged limit ($ \rho_{\rm ch}= 0 $, $ q= 0 $) reproduces the neutral EMSG results of Ref. [27].

It is instructive to identify the physical origin of each term in this system. In the mass equation (24), the first term, $ 4\pi r^2\rho $, represents the standard contribution from the quark matter energy density; the second term, $ qq'/r $, encodes the electromagnetic field energy; and the third term, $ 4\pi\alpha r^2\rho^2(1 + 8p/\rho + 3p^2/\rho^2) $, is a purely EMSG correction arising from the quadratic $ T_{\mu\nu}T^{\mu\nu} $ coupling. Similarly, in the pressure equation (25), the factor $ (1 + 6\alpha p)^{-1} $ in the denominator and the term $ 2\alpha\rho\rho'/(1+6\alpha p) $ are exclusively of EMSG origin, whereas the term $ qq'/(4\pi r^4) $ originates from the electromagnetic sector. The remaining structure—the $ (\rho+p)\psi' $ gravitational coupling and the standard pressure gradient—is inherited directly from general relativity. When $ \alpha \to 0 $, all EMSG corrections vanish identically, and the familiar charged perfect-fluid TOV equations of general relativity are recovered [40, 50].

Numerical integration of this system requires specifying an equation of state $ p= p(\rho) $ relating pressure to energy density, together with a charge distribution profile $ \rho_{\rm ch} = \rho_{\rm ch}(\rho) $, thereby closing the system for the four unknown functions. Regularity at the stellar center requires the boundary conditions:

where $ \rho_c $ denotes the central energy density. The integration proceeds outward from the origin until the pressure vanishes, which defines the stellar surface at radius $ r_{\rm sur} $, satisfying $ p(r_{\rm sur})= 0 $.

To determine the metric potential ψ throughout the interior, an additional boundary condition is required at the surface. Taking the trace of the field equations (14) yields the Ricci scalar $ R= 8\pi(\rho- 3p)[1- 2\alpha(\rho-p)] $, which vanishes in the exterior vacuum region, where $ \rho= p =0 $. This result relies on the convention adopted above that the squared term is constructed from the fluid tensor alone. If the electromagnetic stress had been included in $ T_{\mu\nu}T^{\mu\nu} $, this term would not vanish in the charged exterior, and the external geometry would deviate from the Reissner–Nordström solution. Consequently, the external spacetime remains the Reissner–Nordström solution of general relativity, and continuity of the metric across the surface imposes:

where $ M \equiv m(r_{\rm sur}) $ and $ Q \equiv q(r_{\rm sur}) $ represent the total gravitational mass and electric charge of the configuration, respectively. We emphasize that M is a metric (gravitational) mass parameter, not a separately conserved matter rest mass. The function $ m(r) $ in Eq. (21) is determined from the radial metric component (20) and already includes the matter, electromagnetic, and effective EMSG contributions. Thus, although the material energy-momentum tensor is not separately conserved in EMSG [Eq. (13)], the contracted Bianchi identity $ \nabla^\mu G_{\mu\nu}=0 $ ensures that the total effective source is conserved and that $ m(r) $ is well defined. Because the exterior is Reissner–Nordström, $ M\equiv m(r_{\rm sur}) $ is the mass parameter of the matched asymptotically flat metric. Therefore, it governs the orbital motion of distant test bodies and binary companions, making it the appropriate quantity for comparison with pulsar mass measurements.

To investigate the equilibrium and structure of charged compact stars within the framework of EMSG, we first outline the gravitational theory and incorporate electromagnetic contributions. The Einstein-Hilbert action is generalized to include an additional term quadratic in the energy-momentum tensor:

where g denotes the determinant of the metric tensor, R is the Ricci scalar curvature, and $ T_{\mu\nu} $ is the energy–momentum tensor associated with the matter Lagrangian density $ \mathcal{L}_m $. The electromagnetic sector enters through its Lagrangian $ \mathcal{L}_e $, while the parameter α quantifies the deviation from general relativity and determines the coupling strength of the EMSG correction. In the present work, the quadratic EMSG invariant $ T_{\mu\nu}T^{\mu\nu} $ is constructed solely from the material perfect-fluid energy–momentum tensor. The electromagnetic field is excluded from the squared term and enters only through the Maxwell Lagrangian $ \mathcal{L}_e $ and its stress tensor $ \mathcal{E}_{\mu\nu} $. In particular, $ T_{\mu\nu}T^{\mu\nu} $ contains neither an $ \mathcal{E}_{\mu\nu}\mathcal{E}^{\mu\nu} $ contribution nor mixed matter–electromagnetic contractions.

Variation of this action with respect to the metric yields the modified field equations:

where $ G_{\mu\nu} $ denotes the Einstein tensor, $ \mathcal{E}_{\mu\nu} $ is the electromagnetic energy–momentum tensor, and $ \Theta_{\mu\nu} $ is a new tensor arising from the quadratic matter coupling:

where T denotes the trace of $ T_{\mu\nu} $. The matter energy–momentum tensor is related to the Lagrangian density as follows:

For a charged stellar configuration, we describe matter as a perfect fluid and include the electromagnetic contribution, following the treatment of Ray et al. [42]. The corresponding energy-momentum tensors are:

and

where ρ represents the energy density, p the pressure, $ u^\mu $ the fluid four-velocity, and $ F_{\mu\nu} $ the electromagnetic field-strength tensor. This tensor is constructed from the four-potential $ A_\mu $ as $ F_{\mu\nu} = \nabla_\mu A_\nu - \nabla_\nu A_\mu $, where $ \nabla_\mu $ denotes covariant differentiation. The Maxwell equations governing the electromagnetic field are:

where $ j^\mu = \rho_{\rm ch}u^\mu $ denotes the four-current density, and $ \rho_{\rm ch} $ is the electric charge density.

To derive the modified stellar structure equations, we employ the static, spherically symmetric metric:

where the metric functions ψ and λ depend only on the radial coordinate r. This geometry gives $\sqrt{-g}= {\rm e}^{\psi+ \lambda}r^2\sin\theta$ and $u^\mu = {\rm e}^{-\psi}\delta_0^\mu$. For a purely electric field configuration, the only nonzero component of the field strength tensor is $ F^{01} = -F^{10} $, which, through the Maxwell equation (7), leads to:

where the radial charge function is given by

Since $ \nabla^\mu G_{\mu\nu}= 0 $ identically, taking the covariant divergence of the field equations (2) demonstrates the non-conservation of the matter energy–momentum tensor in EMSG:

which can be expressed, using Eqs. (5) and (6), as:

The explicit form of $ \Theta_{\mu\nu} $ in Eq. (3) depends on the choice of the matter Lagrangian density $ \mathcal{L}_m $. Following established approaches in the literature [52, 53], we adopt $ \mathcal{L}_m = p $ in this analysis, which yields $ \Theta_{\mu\nu} = -(\rho^2+ 4\rho p + 3p^2)u_\mu u_\nu $. Although $ \mathcal{L}_m = -\rho $ is an equally admissible choice in GR, $ \mathcal{L}_m = p $ has been shown to be the physically appropriate prescription for perfect fluids in theories with non-minimal matter-geometry coupling [52, 53]. This choice is also consistently adopted throughout the EMSG literature [27, 28, 30, 35]. Substituting this expression into Eqs. (2) and (13), the field equations and energy-momentum non-conservation take the following forms:

Evaluating these equations in the spherically symmetric geometry defined by Eq. (9), the $ 00 $ and $ 11 $ components yield:

where primes denote derivatives with respect to the radial coordinate. In the absence of electric charge, these expressions reduce to those derived in Ref. [27]. The radial component of the non-conservation equation (15) yields the pressure gradient:

To express these structural relations in a more transparent form, reminiscent of the standard Tolman–Oppenheimer–Volkoff (TOV) formulation, we introduce a gravitational mass function that quantifies the total mass–energy enclosed within a radius r. Integrating Eq. (16) gives:

which can be expressed compactly as:

where the mass function is given by:

This decomposition shows that the effective gravitational mass arises from three distinct contributions: the conventional matter energy density, represented by the first integral; the electromagnetic field energy, encoded in the charge distribution; and the modifications induced by the quadratic matter coupling in EMSG, represented by the third term. Setting $ \alpha= 0 $ recovers the standard charged perfect-fluid mass function in general relativity [45, 50], while the neutral limit ($ q=0 $) reproduces the uncharged EMSG result of Ref. [27].

Using the relation (20), the radial field equation (17) determines the gradient of the metric potential:

Combining Eqs. (11), (18), (21), and (22), we obtain the complete system of modified Tolman-Oppenheimer-Volkoff equations governing charged stellar configurations in EMSG [35]:

These differential equations describe hydrostatic equilibrium in the presence of both electromagnetic fields and EMSG modifications. For $ \alpha =0 $, the standard general relativistic TOV equations are recovered, whereas the uncharged limit ($ \rho_{\rm ch}= 0 $, $ q= 0 $) reproduces the neutral EMSG results of Ref. [27].

It is instructive to identify the physical origin of each term in this system. In the mass equation (24), the first term, $ 4\pi r^2\rho $, represents the standard contribution from the quark matter energy density; the second term, $ qq'/r $, encodes the electromagnetic field energy; and the third term, $ 4\pi\alpha r^2\rho^2(1 + 8p/\rho + 3p^2/\rho^2) $, is a purely EMSG correction arising from the quadratic $ T_{\mu\nu}T^{\mu\nu} $ coupling. Similarly, in the pressure equation (25), the factor $ (1 + 6\alpha p)^{-1} $ in the denominator and the term $ 2\alpha\rho\rho'/(1+6\alpha p) $ are exclusively of EMSG origin, whereas the term $ qq'/(4\pi r^4) $ originates from the electromagnetic sector. The remaining structure—the $ (\rho+p)\psi' $ gravitational coupling and the standard pressure gradient—is inherited directly from general relativity. When $ \alpha \to 0 $, all EMSG corrections vanish identically, and the familiar charged perfect-fluid TOV equations of general relativity are recovered [40, 50].

Numerical integration of this system requires specifying an equation of state $ p= p(\rho) $ relating pressure to energy density, together with a charge distribution profile $ \rho_{\rm ch} = \rho_{\rm ch}(\rho) $, thereby closing the system for the four unknown functions. Regularity at the stellar center requires the boundary conditions:

where $ \rho_c $ denotes the central energy density. The integration proceeds outward from the origin until the pressure vanishes, which defines the stellar surface at radius $ r_{\rm sur} $, satisfying $ p(r_{\rm sur})= 0 $.

To determine the metric potential ψ throughout the interior, an additional boundary condition is required at the surface. Taking the trace of the field equations (14) yields the Ricci scalar $ R= 8\pi(\rho- 3p)[1- 2\alpha(\rho-p)] $, which vanishes in the exterior vacuum region, where $ \rho= p =0 $. This result relies on the convention adopted above that the squared term is constructed from the fluid tensor alone. If the electromagnetic stress had been included in $ T_{\mu\nu}T^{\mu\nu} $, this term would not vanish in the charged exterior, and the external geometry would deviate from the Reissner–Nordström solution. Consequently, the external spacetime remains the Reissner–Nordström solution of general relativity, and continuity of the metric across the surface imposes:

where $ M \equiv m(r_{\rm sur}) $ and $ Q \equiv q(r_{\rm sur}) $ represent the total gravitational mass and electric charge of the configuration, respectively. We emphasize that M is a metric (gravitational) mass parameter, not a separately conserved matter rest mass. The function $ m(r) $ in Eq. (21) is determined from the radial metric component (20) and already includes the matter, electromagnetic, and effective EMSG contributions. Thus, although the material energy-momentum tensor is not separately conserved in EMSG [Eq. (13)], the contracted Bianchi identity $ \nabla^\mu G_{\mu\nu}=0 $ ensures that the total effective source is conserved and that $ m(r) $ is well defined. Because the exterior is Reissner–Nordström, $ M\equiv m(r_{\rm sur}) $ is the mass parameter of the matched asymptotically flat metric. Therefore, it governs the orbital motion of distant test bodies and binary companions, making it the appropriate quantity for comparison with pulsar mass measurements.

A realistic description of deconfined quark matter is essential for modeling ultradense compact stars, particularly when such systems are studied within modified gravity frameworks. In this work, we employ the unified interacting equation of state (EoS) developed in [9], which incorporates perturbative QCD effects, color superconductivity, and a finite strange-quark mass within a single parametric framework. This construction provides a continuous and thermodynamically consistent representation of the 2SC, 2SC+s, and CFL phases and is well suited for analyzing charged quark stars in the context of EMSG theory.

The formulation begins with the thermodynamic potential [5, 6, 9]

where μ denotes the average quark chemical potential, Δ is the gap parameter for color superconductivity, $ m_s $ is the strange quark mass, and $ B_{\rm eff} $ is the effective bag constant, which encodes the nonperturbative QCD vacuum contribution [8]. The phase-dependent parameters $ (\xi_4,\xi_{2a},\xi_{2b}) $ take the following values:

Using the standard thermodynamic relations $ p=-\Omega $ and $ \rho=\Omega+\mu n_q $, we derive the unified interacting EoS by introducing an interaction parameter

where $ a_4 $ characterizes the strength of perturbative QCD corrections arising from one-gluon exchange at $ \mathcal{O}(\alpha_s^2) $. The resulting pressure–density relation can be written in the compact form

which smoothly interpolates between the noninteracting MIT bag-model limit at $ \lambda \to 0 $ and the strongly interacting regime at $ \lambda \to \infty $, in which the EoS approaches $ p=\rho - 2B_{\rm eff} $. Here, the prefactor $ 1/3 $ multiplies only the bag-model contribution $ (\rho-4B_{\rm eff}) $, whereas the interaction contribution is a separate additive term, following the parenthesization of Ref. [9]. For the positive interaction branch used in the present numerical analysis, expansion of the square root for large positive λ, with $ x= 3\pi^{2}(\rho-B_{\rm eff})/\lambda^{2} $, gives $ \sqrt{1+x}=1+x/2-x^{2}/8+\mathcal{O}(x^{3}) $. Hence, the interaction term first tends to ${2}/{3}(\rho-B_{\rm eff})$, and the full equation of state becomes

which yields $ p\to\rho-2B_{\rm eff} $ and $ {\rm d}p/{\rm d}\rho\to1 $ in the strongly interacting limit.

For numerical convenience, we adopt the rescaled variables

yielding the dimensionless representation

The parameter $ \bar{\lambda} $ effectively governs the stiffness of the EoS. Larger values correspond to stronger interactions between quarks, enhanced color superconductivity, or a reduced strange quark mass. Coupling this EoS to the modified Tolman-Oppenheimer-Volkoff equations in EMSG yields the mass-radius sequences examined later in this work. Notably, the unified parametrization avoids inconsistencies arising from the independent tuning of microscopic parameters, ensuring that the resulting stellar models remain compatible with observational constraints from massive pulsars and gravitational wave events.

A realistic description of deconfined quark matter is essential for modeling ultradense compact stars, particularly when such systems are studied within modified gravity frameworks. In this work, we employ the unified interacting equation of state (EoS) developed in [9], which incorporates perturbative QCD effects, color superconductivity, and a finite strange-quark mass within a single parametric framework. This construction provides a continuous and thermodynamically consistent representation of the 2SC, 2SC+s, and CFL phases and is well suited for analyzing charged quark stars in the context of EMSG theory.

The formulation begins with the thermodynamic potential [5, 6, 9]

where μ denotes the average quark chemical potential, Δ is the gap parameter for color superconductivity, $ m_s $ is the strange quark mass, and $ B_{\rm eff} $ is the effective bag constant, which encodes the nonperturbative QCD vacuum contribution [8]. The phase-dependent parameters $ (\xi_4,\xi_{2a},\xi_{2b}) $ take the following values:

Using the standard thermodynamic relations $ p=-\Omega $ and $ \rho=\Omega+\mu n_q $, we derive the unified interacting EoS by introducing an interaction parameter

where $ a_4 $ characterizes the strength of perturbative QCD corrections arising from one-gluon exchange at $ \mathcal{O}(\alpha_s^2) $. The resulting pressure–density relation can be written in the compact form

which smoothly interpolates between the noninteracting MIT bag-model limit at $ \lambda \to 0 $ and the strongly interacting regime at $ \lambda \to \infty $, in which the EoS approaches $ p=\rho - 2B_{\rm eff} $. Here, the prefactor $ 1/3 $ multiplies only the bag-model contribution $ (\rho-4B_{\rm eff}) $, whereas the interaction contribution is a separate additive term, following the parenthesization of Ref. [9]. For the positive interaction branch used in the present numerical analysis, expansion of the square root for large positive λ, with $ x= 3\pi^{2}(\rho-B_{\rm eff})/\lambda^{2} $, gives $ \sqrt{1+x}=1+x/2-x^{2}/8+\mathcal{O}(x^{3}) $. Hence, the interaction term first tends to ${2}/{3}(\rho-B_{\rm eff})$, and the full equation of state becomes

which yields $ p\to\rho-2B_{\rm eff} $ and $ {\rm d}p/{\rm d}\rho\to1 $ in the strongly interacting limit.

For numerical convenience, we adopt the rescaled variables

yielding the dimensionless representation

The parameter $ \bar{\lambda} $ effectively governs the stiffness of the EoS. Larger values correspond to stronger interactions between quarks, enhanced color superconductivity, or a reduced strange quark mass. Coupling this EoS to the modified Tolman-Oppenheimer-Volkoff equations in EMSG yields the mass-radius sequences examined later in this work. Notably, the unified parametrization avoids inconsistencies arising from the independent tuning of microscopic parameters, ensuring that the resulting stellar models remain compatible with observational constraints from massive pulsars and gravitational wave events.

To complement the equation of state describing ultra-dense quark matter, an explicit prescription for the electric charge density $ \rho_{\rm ch} $ is required when modeling charged stellar configurations. A commonly adopted ansatz, originally proposed in the context of general relativity by Ray et al. [42], assumes that the charge distribution follows the local energy density. This proportionality is expressed as

where the dimensionless constant β regulates the total amount of electric charge contained in the fluid. The physical motivation for this choice is straightforward: regions of higher mass-energy can support greater charge accumulation, making this relation a natural first approximation for compact stars [49]. In the ultra-dense quark star interior, this proportionality is further motivated by the fact that the local quark number density, and hence the net electric charge arising from the slight imbalance between quark flavours, scales directly with the energy density under the extreme conditions prevailing in the stellar core. Thus, $ \rho_{\rm ch} = \beta\rho $ provides a physically transparent and thermodynamically consistent first-order approximation for the charge distribution [42, 50].

This charge profile has been widely employed in recent investigations of charged quark stars within modified gravity frameworks, including analyses in metric $ f(R) $ gravity [54], $ f(R,T) $ gravity [55], and the regularized $ 4D $ Einstein-Gauss-Bonnet theory [14, 56]. Its use in the present EMSG context facilitates consistent comparison with these previous studies while providing a simple, physically motivated mechanism for incorporating electric charge into the stellar structure equations. We acknowledge that this proportionality is a phenomenological ansatz rather than a result derived from first-principles QCD calculations. A rigorous microscopic determination of the charge distribution in dense quark matter remains an open problem in the field. The ansatz $ \rho_{\rm ch} = \beta\rho $ represents the standard approach adopted in the charged compact star literature [14, 42, 45, 50, 54, 55], and its use here is consistent with this established practice. Exploring more sophisticated charge distributions, potentially informed by QCD calculations at finite density, would be a valuable direction for future work. Throughout this work, we fix the charge fraction at $ \beta=0.5 $, a representative intermediate value that lies well within the physically admissible range $ 0\leqslant\beta<1 $ adopted in the charged compact star literature [42, 45]. In this range, $ \beta=0 $ recovers the neutral limit, whereas $ \beta\to1 $ approaches the extremal limit ($ Q/M\to1 $). This choice yields a sizeable yet sub-extremal electric charge ($ Q/M\sim0.6 $), large enough to demonstrate the structural role of the Coulomb sector while keeping the configurations gravitationally bound. Increasing the charge fraction enhances the electric charge and the associated outward Coulomb support, thereby increasing the maximum mass and radius. This effect is quantified by the neutral ($ \beta=0 $) and charged ($ \beta=0.5 $) sequences in Tables 1 and 2. For example, at EMSG coupling $ \alpha=0 $ and $ \bar{\lambda}=0.1 $, the maximum mass increases from $ M_{\max}^{\rm n}=2.26\,M_\odot $ to $ M_{\max}^{\rm ch}=2.88\,M_\odot $ when β is raised from $ 0 $ to $ 0.5 $. A systematic scan over intermediate β values would provide a finer sensitivity map and is left for future work.

To complement the equation of state describing ultra-dense quark matter, an explicit prescription for the electric charge density $ \rho_{\rm ch} $ is required when modeling charged stellar configurations. A commonly adopted ansatz, originally proposed in the context of general relativity by Ray et al. [42], assumes that the charge distribution follows the local energy density. This proportionality is expressed as

where the dimensionless constant β regulates the total amount of electric charge contained in the fluid. The physical motivation for this choice is straightforward: regions of higher mass-energy can support greater charge accumulation, making this relation a natural first approximation for compact stars [49]. In the ultra-dense quark star interior, this proportionality is further motivated by the fact that the local quark number density, and hence the net electric charge arising from the slight imbalance between quark flavours, scales directly with the energy density under the extreme conditions prevailing in the stellar core. Thus, $ \rho_{\rm ch} = \beta\rho $ provides a physically transparent and thermodynamically consistent first-order approximation for the charge distribution [42, 50].

This charge profile has been widely employed in recent investigations of charged quark stars within modified gravity frameworks, including analyses in metric $ f(R) $ gravity [54], $ f(R,T) $ gravity [55], and the regularized $ 4D $ Einstein-Gauss-Bonnet theory [14, 56]. Its use in the present EMSG context facilitates consistent comparison with these previous studies while providing a simple, physically motivated mechanism for incorporating electric charge into the stellar structure equations. We acknowledge that this proportionality is a phenomenological ansatz rather than a result derived from first-principles QCD calculations. A rigorous microscopic determination of the charge distribution in dense quark matter remains an open problem in the field. The ansatz $ \rho_{\rm ch} = \beta\rho $ represents the standard approach adopted in the charged compact star literature [14, 42, 45, 50, 54, 55], and its use here is consistent with this established practice. Exploring more sophisticated charge distributions, potentially informed by QCD calculations at finite density, would be a valuable direction for future work. Throughout this work, we fix the charge fraction at $ \beta=0.5 $, a representative intermediate value that lies well within the physically admissible range $ 0\leqslant\beta<1 $ adopted in the charged compact star literature [42, 45]. In this range, $ \beta=0 $ recovers the neutral limit, whereas $ \beta\to1 $ approaches the extremal limit ($ Q/M\to1 $). This choice yields a sizeable yet sub-extremal electric charge ($ Q/M\sim0.6 $), large enough to demonstrate the structural role of the Coulomb sector while keeping the configurations gravitationally bound. Increasing the charge fraction enhances the electric charge and the associated outward Coulomb support, thereby increasing the maximum mass and radius. This effect is quantified by the neutral ($ \beta=0 $) and charged ($ \beta=0.5 $) sequences in Tables 1 and 2. For example, at EMSG coupling $ \alpha=0 $ and $ \bar{\lambda}=0.1 $, the maximum mass increases from $ M_{\max}^{\rm n}=2.26\,M_\odot $ to $ M_{\max}^{\rm ch}=2.88\,M_\odot $ when β is raised from $ 0 $ to $ 0.5 $. A systematic scan over intermediate β values would provide a finer sensitivity map and is left for future work.

In this section, we present numerical solutions of the modified TOV equations for charged quark stars within the EMSG framework. Using the unified interacting quark matter equation of state and a density-proportional charge profile, we examine how the coupling parameter α, the interaction parameter $ \bar{\lambda} $, and the charge fraction β affect the global stellar properties. In particular, we analyze the resulting mass-radius relations, total electric charge, compactness, and internal stability indicators that characterize the equilibrium sequences obtained in this model. The numerical integrations are performed over the following parameter ranges: the EMSG coupling strength takes the values $ \alpha \in \{-2.0,\,-1.0,\,0.0,\,+1.0,\,+2.0\} \times 10^{-38} \text{cm}^3\; \text{erg}^{-1} $, the quark matter interaction parameter spans $ \bar{\lambda} \in \{0.02,\,0.04,\,0.06,\,0.08,\,0.10\} $, and the charge fraction is fixed at $ \beta = 0.5 $ throughout, with bag constant $ B_{\rm eff} = 60\; \text{MeV fm}^{-3} $. For each combination of these parameters, we have verified that the central energy density $ \rho_c $ and central pressure $ p_c $ remain strictly positive, confirming that all computed configurations correspond to physically admissible stellar interiors without pathological behavior at the origin.

In this section, we present numerical solutions of the modified TOV equations for charged quark stars within the EMSG framework. Using the unified interacting quark matter equation of state and a density-proportional charge profile, we examine how the coupling parameter α, the interaction parameter $ \bar{\lambda} $, and the charge fraction β affect the global stellar properties. In particular, we analyze the resulting mass-radius relations, total electric charge, compactness, and internal stability indicators that characterize the equilibrium sequences obtained in this model. The numerical integrations are performed over the following parameter ranges: the EMSG coupling strength takes the values $ \alpha \in \{-2.0,\,-1.0,\,0.0,\,+1.0,\,+2.0\} \times 10^{-38} \text{cm}^3\; \text{erg}^{-1} $, the quark matter interaction parameter spans $ \bar{\lambda} \in \{0.02,\,0.04,\,0.06,\,0.08,\,0.10\} $, and the charge fraction is fixed at $ \beta = 0.5 $ throughout, with bag constant $ B_{\rm eff} = 60\; \text{MeV fm}^{-3} $. For each combination of these parameters, we have verified that the central energy density $ \rho_c $ and central pressure $ p_c $ remain strictly positive, confirming that all computed configurations correspond to physically admissible stellar interiors without pathological behavior at the origin.

Varying the EMSG coupling constant α produces clear and systematic changes in the equilibrium structure of both neutral and charged quark stars, as illustrated in Fig. 1 and summarized in Table 1. Negative values of α strengthen the effective pressure contribution generated by the nonlinear matter–geometry coupling, enabling the star to support larger gravitational masses and radii at lower central densities. This trend is apparent in the neutral sequences, but it becomes even more pronounced when electric charge is included. Coulomb repulsion amplifies the stabilizing effect of negative α, pushing the maximum mass of charged stars well beyond the general relativistic limit. For example, configurations with $ \alpha = -2\mu $ can attain nearly $ 3M_\odot $, with compactness values approaching $ M/R\simeq 0.34 $. This regime is compatible with the mass range of the secondary compact object in GW190814 [61], whose precise nature remains uncertain but which has been widely discussed as a possible quark star or heavy neutron star candidate in the literature [30, 34, 35]. A black-hole interpretation of this object cannot be excluded either. Neutral stars over the same parameter range remain bounded by $ M_{\rm max}\lesssim 2.3M_\odot $, demonstrating that although modifications to gravity alone enhance stability, electromagnetic forces provide an additional outward contribution required to achieve the highest observed stellar masses.

Figure 1.  (color online) Mass versus radius and charge versus radius relations for quark stars in energy-momentum squared gravity (EMSG), computed using the unified interacting quark matter equation of state. Each curve corresponds to a different value of the matter-geometry coupling parameter α, expressed in units of $ 10^{-38}\; \text{cm}^3\; \text{erg}^{-1} $: $ \alpha = -2 $ (cyan), $ -1 $ (green), $ 0 $ (black dashed, corresponding to general relativity), $ +1 $ (orange), and $ +2 $ (red). All sequences are generated with the interaction parameter $ \bar{\lambda} = 0.1 $, the charge fraction $ \beta=0.5 $ and bag constant $ B_{\rm eff} = 60\; \text{MeV}\; \text{fm}^{-3} $. Top left: Mass-radius curves for neutral stars. Top right: Mass-radius curves for charged stars with a charge density proportional to the energy density. Horizontal shaded bands show observational constraints from PSR J1614-2230 with $ 1.97 \pm 0.04 \,M_{\odot} $ [57], PSR J0952$ - $0607 with $ M = 2.35 \pm 0.17\,M_{\odot} $ [58], PSR J0740+6620 with $ M = 2.08^{+0.07}_{-0.07}\,M_{\odot} $ [59], PSR J0348+0432 with $ M = 2.01 \pm 0.04\,M_{\odot} $ [60], and the mass range associated with the secondary component of GW190814 [61]. Bottom left: Total electric charge as a function of radius. Bottom right: Compactness $ M/R $ as a function of gravitational mass. Black circles indicate the maximum mass point on each sequence for the corresponding parameter choice.

The influence of α is also evident in the total electric charge profiles. As shown in the bottom-left panel of Fig. 1, the charge distribution increases monotonically with radius and reaches its maximum at the stellar surface, consistent with the proportionality relation $ \rho_{\rm ch}=\beta\rho $. For negative α, the enhanced pressure support allows the star to extend to slightly larger radii at comparable central density, resulting in higher total charge values; for instance, $ Q_{\max}\approx 3.3\times10^{20}\,{\rm C} $ at $ \alpha=-2\mu $. As α increases toward positive values, both the stellar radius and the integrated charge decrease systematically, reducing the surface charge to $ Q_{\max}\approx 2.9\times10^{20}\,{\rm C} $ for $ \alpha=+2\mu $. These variations reflect the combined role of matter–geometry coupling and the charge-density prescription in shaping the star's electromagnetic content. The surface charge values obtained here, $ Q_{\max} \sim 2.9 $–$ 3.3 \times 10^{20} $ C, are consistent with those reported in recent investigations of charged compact stars in modified gravity frameworks. Studies of charged compact stars in non-minimally coupled gravity and in matter–geometry coupled theories have reported surface charges of comparable order of magnitude [62, 63], as have analyses of charged quark stars in $ f(R) $ and $ f(R,T) $ gravity [54, 55] and in regularized four-dimensional Einstein–Gauss–Bonnet theory [56]. The agreement across these independent frameworks indicates that the charge magnitudes obtained in the present work are comparable to those adopted in the broad class of charged compact-star models explored in the recent literature. We emphasize, however, that a net charge of this order is an idealization: a surface charge $ Q\sim3\times10^{20}\; {\rm C} $ on a star of radius $ R\sim13\; {\rm km} $ produces a surface field $ E_{\rm surf}\sim 1.7\times 10^{22}\; {\rm V\,m^{-1}} $, exceeding the Schwinger critical field $ E_{\rm c}\simeq1.3\times10^{18}\; {\rm V\,m^{-1}} $ for electron–positron pair production by roughly four orders of magnitude. In a realistic quantum vacuum, such fields would drive prolific pair creation and rapid discharge, so charges of this magnitude are not expected to be sustained. The maximal charges adopted here should therefore be understood as bracketing the largest possible influence of electric charge on the equilibrium structure within EMSG, rather than as necessarily sustainable, globally unscreened charges in isolated compact stars. A realistic star may carry a much smaller net charge, and a quantitative assessment of how these trends persist would require a dedicated lower-charge or discharge-regulated analysis beyond the scope of the present equilibrium treatment.

Across the full range of α, all computed stellar sequences remain consistent with current astrophysical constraints, including the mass and radius measurements from NICER for PSR J1614$ - $2230, J0348$ + $0432, J0740$ + $6620, and J0952$ - $0607, as well as the bounds inferred from multimessenger gravitational-wave observations. Collectively, these results highlight the complementary influence of modified gravity and moderate electric charge in governing the structure of ultra-dense matter and emphasize the potential of observational data to distinguish general relativity from alternative theories such as EMSG in the strong-field regime.

Varying the EMSG coupling constant α produces clear and systematic changes in the equilibrium structure of both neutral and charged quark stars, as illustrated in Fig. 1 and summarized in Table 1. Negative values of α strengthen the effective pressure contribution generated by the nonlinear matter–geometry coupling, enabling the star to support larger gravitational masses and radii at lower central densities. This trend is apparent in the neutral sequences, but it becomes even more pronounced when electric charge is included. Coulomb repulsion amplifies the stabilizing effect of negative α, pushing the maximum mass of charged stars well beyond the general relativistic limit. For example, configurations with $ \alpha = -2\mu $ can attain nearly $ 3M_\odot $, with compactness values approaching $ M/R\simeq 0.34 $. This regime is compatible with the mass range of the secondary compact object in GW190814 [61], whose precise nature remains uncertain but which has been widely discussed as a possible quark star or heavy neutron star candidate in the literature [30, 34, 35]. A black-hole interpretation of this object cannot be excluded either. Neutral stars over the same parameter range remain bounded by $ M_{\rm max}\lesssim 2.3M_\odot $, demonstrating that although modifications to gravity alone enhance stability, electromagnetic forces provide an additional outward contribution required to achieve the highest observed stellar masses.

Figure 1.  (color online) Mass versus radius and charge versus radius relations for quark stars in energy-momentum squared gravity (EMSG), computed using the unified interacting quark matter equation of state. Each curve corresponds to a different value of the matter-geometry coupling parameter α, expressed in units of $ 10^{-38}\; \text{cm}^3\; \text{erg}^{-1} $: $ \alpha = -2 $ (cyan), $ -1 $ (green), $ 0 $ (black dashed, corresponding to general relativity), $ +1 $ (orange), and $ +2 $ (red). All sequences are generated with the interaction parameter $ \bar{\lambda} = 0.1 $, the charge fraction $ \beta=0.5 $ and bag constant $ B_{\rm eff} = 60\; \text{MeV}\; \text{fm}^{-3} $. Top left: Mass-radius curves for neutral stars. Top right: Mass-radius curves for charged stars with a charge density proportional to the energy density. Horizontal shaded bands show observational constraints from PSR J1614-2230 with $ 1.97 \pm 0.04 \,M_{\odot} $ [57], PSR J0952$ - $0607 with $ M = 2.35 \pm 0.17\,M_{\odot} $ [58], PSR J0740+6620 with $ M = 2.08^{+0.07}_{-0.07}\,M_{\odot} $ [59], PSR J0348+0432 with $ M = 2.01 \pm 0.04\,M_{\odot} $ [60], and the mass range associated with the secondary component of GW190814 [61]. Bottom left: Total electric charge as a function of radius. Bottom right: Compactness $ M/R $ as a function of gravitational mass. Black circles indicate the maximum mass point on each sequence for the corresponding parameter choice.

The influence of α is also evident in the total electric charge profiles. As shown in the bottom-left panel of Fig. 1, the charge distribution increases monotonically with radius and reaches its maximum at the stellar surface, consistent with the proportionality relation $ \rho_{\rm ch}=\beta\rho $. For negative α, the enhanced pressure support allows the star to extend to slightly larger radii at comparable central density, resulting in higher total charge values; for instance, $ Q_{\max}\approx 3.3\times10^{20}\,{\rm C} $ at $ \alpha=-2\mu $. As α increases toward positive values, both the stellar radius and the integrated charge decrease systematically, reducing the surface charge to $ Q_{\max}\approx 2.9\times10^{20}\,{\rm C} $ for $ \alpha=+2\mu $. These variations reflect the combined role of matter–geometry coupling and the charge-density prescription in shaping the star's electromagnetic content. The surface charge values obtained here, $ Q_{\max} \sim 2.9 $–$ 3.3 \times 10^{20} $ C, are consistent with those reported in recent investigations of charged compact stars in modified gravity frameworks. Studies of charged compact stars in non-minimally coupled gravity and in matter–geometry coupled theories have reported surface charges of comparable order of magnitude [62, 63], as have analyses of charged quark stars in $ f(R) $ and $ f(R,T) $ gravity [54, 55] and in regularized four-dimensional Einstein–Gauss–Bonnet theory [56]. The agreement across these independent frameworks indicates that the charge magnitudes obtained in the present work are comparable to those adopted in the broad class of charged compact-star models explored in the recent literature. We emphasize, however, that a net charge of this order is an idealization: a surface charge $ Q\sim3\times10^{20}\; {\rm C} $ on a star of radius $ R\sim13\; {\rm km} $ produces a surface field $ E_{\rm surf}\sim 1.7\times 10^{22}\; {\rm V\,m^{-1}} $, exceeding the Schwinger critical field $ E_{\rm c}\simeq1.3\times10^{18}\; {\rm V\,m^{-1}} $ for electron–positron pair production by roughly four orders of magnitude. In a realistic quantum vacuum, such fields would drive prolific pair creation and rapid discharge, so charges of this magnitude are not expected to be sustained. The maximal charges adopted here should therefore be understood as bracketing the largest possible influence of electric charge on the equilibrium structure within EMSG, rather than as necessarily sustainable, globally unscreened charges in isolated compact stars. A realistic star may carry a much smaller net charge, and a quantitative assessment of how these trends persist would require a dedicated lower-charge or discharge-regulated analysis beyond the scope of the present equilibrium treatment.

Across the full range of α, all computed stellar sequences remain consistent with current astrophysical constraints, including the mass and radius measurements from NICER for PSR J1614$ - $2230, J0348$ + $0432, J0740$ + $6620, and J0952$ - $0607, as well as the bounds inferred from multimessenger gravitational-wave observations. Collectively, these results highlight the complementary influence of modified gravity and moderate electric charge in governing the structure of ultra-dense matter and emphasize the potential of observational data to distinguish general relativity from alternative theories such as EMSG in the strong-field regime.

The interaction parameter $ \bar{\lambda} $ of the unified interacting quark matter equation of state plays a central role in determining the stiffness of dense matter and, consequently, the global properties of quark stars. The trends associated with variations in $ \bar{\lambda} $ are shown in Fig. 2 and summarized in Table 2. As shown in Table 2, the central pressure $ p_c $ at the maximum-mass configuration decreases from $ 241.17\; \text{MeV fm}^{-3} $ at $ \bar{\lambda}=0.02 $ to $ 220.39\; \text{MeV fm}^{-3} $ at $ \bar{\lambda}=0.10 $. This decrease reflects the progressive stiffening of the equation of state and the corresponding reduction in the central pressure required to support the maximum-mass star as the interaction parameter increases. Smaller values of $ \bar{\lambda} $ correspond to a softer equation of state, producing stars that are less massive and more compact. In contrast, larger values yield progressively stiffer matter capable of supporting higher gravitational masses and slightly larger radii. This behavior is clearly evident in the neutral sequences: the maximum mass increases from $ M_{\rm max}^n=2.08\,M_\odot $ at $ \bar{\lambda}=0.02 $ to $ 2.25\,M_\odot $ at $ \bar{\lambda}=0.10 $, accompanied by a systematic decrease in the central energy density at the stability limit.

Figure 2.  (color online) Mass versus radius and charge versus radius relations for quark star configurations in energy-momentum squared gravity, obtained by varying the interaction parameter $ \bar{\lambda} $ of the unified interacting quark matter equation of state. All models employ $ B_{\text{eff}} = 60\; \text{MeV}\; \text{fm}^{-3} $, a fixed matter-geometry coupling $ \alpha = 1.0\mu $ with $ \mu = 10^{-38}\; \text{cm}^3\; \text{erg}^{-1} $, and the charge fraction $ \beta=0.5 $. Five representative values of $ \bar{\lambda} $ are shown: $ 0.02 $ (cyan), $ 0.04 $ (green), $ 0.06 $ (yellow), $ 0.08 $ (orange), and $ 0.10 $ (purple). Upper left panel: Mass-radius sequences for electrically neutral stars. Upper right panel: Mass-radius sequences for charged stars. Horizontal shaded bands indicate observational constraints from PSR J1614-2230 [57], J0348+0432 [60], J0740+6620 [59], J0952-0607 [58], and the mass range associated with the secondary component of GW190814 [61]. Lower left panel: Total electric charge as a function of radius for the charged configurations. Lower right panel: Compactness $ M/R $ plotted against gravitational mass. Black circles mark the maximum mass configuration for each choice of $ \bar{\lambda} $.

The impact of $ \bar{\lambda} $ becomes more pronounced for charged stars. As this parameter increases, the additional matter stiffness enhances the support provided by the Coulomb field, allowing charged configurations to reach larger masses within the same gravitational framework. For instance, the maximum mass increases from $ M_{\rm max}^{\rm ch}=2.62\,M_\odot $ at $ \bar{\lambda}=0.02 $ to $ 2.84\,M_\odot $ at $ \bar{\lambda}=0.10 $, while the corresponding stellar radius expands from approximately $ 12.0 $ km to $ 12.8 $ km. These variations reflect the combined influence of the quark matter interaction strength and the electric charge contributions, both of which stiffen the stellar interior and suppress excessive compaction.

The total charge profiles shown in the lower-left panel of Fig. 2 further support this trend. Higher values of $ \bar{\lambda} $ produce configurations that sustain larger integrated charge and, correspondingly, higher surface charge, with $ Q_{\max} $ increasing monotonically with $ \bar{\lambda} $, from $ 2.74\times10^{20}\; \text{C} $ at $ \bar{\lambda}=0.02 $ to $ 2.98\times10^{20}\; \text{C} $ at $ \bar{\lambda}=0.10 $. This behavior parallels the changes observed in the mass-radius relations and underscores the sensitivity of the Coulomb sector to the stiffness of quark matter.

Importantly, all stellar sequences across the considered range of $ \bar{\lambda} $ remain compatible with current observational constraints. Neutral configurations comfortably satisfy the mass measurements from NICER and radio pulsars, while the charged models with moderate $ \bar{\lambda} $ values extend into the mass range associated with the secondary component, under a possible compact-star interpretation, in GW190814. Taken together, these results highlight the role of quark matter interactions in shaping the macroscopic properties of quark stars and demonstrate how varying $ \bar{\lambda} $ provides a complementary mechanism—alongside the EMSG coupling α and electric charge fraction β—for accessing the broader parameter space of ultra-dense compact stars in modified gravity.

To facilitate independent reproduction of the models presented here, we highlight two representative maximum-mass configurations from Table 1. For the strongest negative EMSG coupling, $ \alpha = -2.0\,\mu $, with $ \bar{\lambda} = 0.1 $, $ \beta = 0.5 $, and $ B_{\rm eff} = 60\; \text{MeV fm}^{-3} $, the maximum-mass charged star has a central density $ \rho_c = 968\; \text{MeV fm}^{-3} $, central pressure $ p_c = 277.77\; \text{MeV fm}^{-3} $, gravitational mass $ M = 2.98\,M_{\odot} $, radius $ R = 12.91\; \text{km} $, and total surface charge $ Q = 3.30\times10^{20}\; \text{C} $. The corresponding general-relativistic configuration ($ \alpha = 0 $) has $ \rho_c = 844\; \text{MeV fm}^{-3} $, $ p_c = 233.88\; \text{MeV fm}^{-3} $, $ M = 2.88\,M_{\odot} $, $ R = 12.84\; \text{km} $, and $ Q = 3.07\times10^{20}\; \text{C} $. Together, these two models bracket the range of EMSG effects examined in this work and provide sufficient information for independent numerical verification.

The interaction parameter $ \bar{\lambda} $ of the unified interacting quark matter equation of state plays a central role in determining the stiffness of dense matter and, consequently, the global properties of quark stars. The trends associated with variations in $ \bar{\lambda} $ are shown in Fig. 2 and summarized in Table 2. As shown in Table 2, the central pressure $ p_c $ at the maximum-mass configuration decreases from $ 241.17\; \text{MeV fm}^{-3} $ at $ \bar{\lambda}=0.02 $ to $ 220.39\; \text{MeV fm}^{-3} $ at $ \bar{\lambda}=0.10 $. This decrease reflects the progressive stiffening of the equation of state and the corresponding reduction in the central pressure required to support the maximum-mass star as the interaction parameter increases. Smaller values of $ \bar{\lambda} $ correspond to a softer equation of state, producing stars that are less massive and more compact. In contrast, larger values yield progressively stiffer matter capable of supporting higher gravitational masses and slightly larger radii. This behavior is clearly evident in the neutral sequences: the maximum mass increases from $ M_{\rm max}^n=2.08\,M_\odot $ at $ \bar{\lambda}=0.02 $ to $ 2.25\,M_\odot $ at $ \bar{\lambda}=0.10 $, accompanied by a systematic decrease in the central energy density at the stability limit.

Figure 2.  (color online) Mass versus radius and charge versus radius relations for quark star configurations in energy-momentum squared gravity, obtained by varying the interaction parameter $ \bar{\lambda} $ of the unified interacting quark matter equation of state. All models employ $ B_{\text{eff}} = 60\; \text{MeV}\; \text{fm}^{-3} $, a fixed matter-geometry coupling $ \alpha = 1.0\mu $ with $ \mu = 10^{-38}\; \text{cm}^3\; \text{erg}^{-1} $, and the charge fraction $ \beta=0.5 $. Five representative values of $ \bar{\lambda} $ are shown: $ 0.02 $ (cyan), $ 0.04 $ (green), $ 0.06 $ (yellow), $ 0.08 $ (orange), and $ 0.10 $ (purple). Upper left panel: Mass-radius sequences for electrically neutral stars. Upper right panel: Mass-radius sequences for charged stars. Horizontal shaded bands indicate observational constraints from PSR J1614-2230 [57], J0348+0432 [60], J0740+6620 [59], J0952-0607 [58], and the mass range associated with the secondary component of GW190814 [61]. Lower left panel: Total electric charge as a function of radius for the charged configurations. Lower right panel: Compactness $ M/R $ plotted against gravitational mass. Black circles mark the maximum mass configuration for each choice of $ \bar{\lambda} $.

The impact of $ \bar{\lambda} $ becomes more pronounced for charged stars. As this parameter increases, the additional matter stiffness enhances the support provided by the Coulomb field, allowing charged configurations to reach larger masses within the same gravitational framework. For instance, the maximum mass increases from $ M_{\rm max}^{\rm ch}=2.62\,M_\odot $ at $ \bar{\lambda}=0.02 $ to $ 2.84\,M_\odot $ at $ \bar{\lambda}=0.10 $, while the corresponding stellar radius expands from approximately $ 12.0 $ km to $ 12.8 $ km. These variations reflect the combined influence of the quark matter interaction strength and the electric charge contributions, both of which stiffen the stellar interior and suppress excessive compaction.

The total charge profiles shown in the lower-left panel of Fig. 2 further support this trend. Higher values of $ \bar{\lambda} $ produce configurations that sustain larger integrated charge and, correspondingly, higher surface charge, with $ Q_{\max} $ increasing monotonically with $ \bar{\lambda} $, from $ 2.74\times10^{20}\; \text{C} $ at $ \bar{\lambda}=0.02 $ to $ 2.98\times10^{20}\; \text{C} $ at $ \bar{\lambda}=0.10 $. This behavior parallels the changes observed in the mass-radius relations and underscores the sensitivity of the Coulomb sector to the stiffness of quark matter.

Importantly, all stellar sequences across the considered range of $ \bar{\lambda} $ remain compatible with current observational constraints. Neutral configurations comfortably satisfy the mass measurements from NICER and radio pulsars, while the charged models with moderate $ \bar{\lambda} $ values extend into the mass range associated with the secondary component, under a possible compact-star interpretation, in GW190814. Taken together, these results highlight the role of quark matter interactions in shaping the macroscopic properties of quark stars and demonstrate how varying $ \bar{\lambda} $ provides a complementary mechanism—alongside the EMSG coupling α and electric charge fraction β—for accessing the broader parameter space of ultra-dense compact stars in modified gravity.

To facilitate independent reproduction of the models presented here, we highlight two representative maximum-mass configurations from Table 1. For the strongest negative EMSG coupling, $ \alpha = -2.0\,\mu $, with $ \bar{\lambda} = 0.1 $, $ \beta = 0.5 $, and $ B_{\rm eff} = 60\; \text{MeV fm}^{-3} $, the maximum-mass charged star has a central density $ \rho_c = 968\; \text{MeV fm}^{-3} $, central pressure $ p_c = 277.77\; \text{MeV fm}^{-3} $, gravitational mass $ M = 2.98\,M_{\odot} $, radius $ R = 12.91\; \text{km} $, and total surface charge $ Q = 3.30\times10^{20}\; \text{C} $. The corresponding general-relativistic configuration ($ \alpha = 0 $) has $ \rho_c = 844\; \text{MeV fm}^{-3} $, $ p_c = 233.88\; \text{MeV fm}^{-3} $, $ M = 2.88\,M_{\odot} $, $ R = 12.84\; \text{km} $, and $ Q = 3.07\times10^{20}\; \text{C} $. Together, these two models bracket the range of EMSG effects examined in this work and provide sufficient information for independent numerical verification.

Figure 3 shows the radial profile of the charge density $ \rho_{\rm ch}(r) = \beta\,\rho(r) $ obtained from the numerical integration of the modified TOV equations. Each curve is evaluated at the central density corresponding to the maximum-mass configuration for that parameter set, as listed in Tables 1 and 2. In both panels, $ \rho_{\rm ch} $ reaches its peak value at the stellar centre, $ \rho_{\rm ch}(0) = \beta\,\rho_c $, where $ \rho_c $ differs for each curve according to the maximum-mass entry in the corresponding table. It then decreases monotonically outward while remaining strictly positive throughout the interior. The profile is finite and smooth at $ r = 0 $, consistent with the regularity condition imposed by the boundary condition $ q(0) = 0 $. At the stellar surface, $ \rho_{\rm ch} $ does not vanish but approaches a nonzero value determined by the finite surface density of quark matter. This is a characteristic feature of the bag-model equation of state, for which the pressure vanishes at a nonzero energy density set by $ B_{\rm eff} $. In the left panel, more negative values of α correspond to higher maximum-mass central densities and therefore produce larger central charge densities, while also extending the star to slightly larger radii. In the right panel, smaller values of $ \bar{\lambda} $ are associated with softer equations of state and higher maximum-mass central densities, yielding larger central charge densities but shorter stellar radii. In all cases, the charge density remains positive definite, confirming that the adopted ansatz $ \rho_{\rm ch} = \beta\rho $ with $ \beta \gt 0 $ introduces a physically consistent, positive-definite charge distribution throughout the stellar interior.

Figure 3.  (color online) Gravitational mass M as a function of the central energy density $ \rho_c $ for charged quark star configurations in EMSG. The parameter spaces used in the left and right panels are identical to those adopted in Fig. 1 and Fig. 2, respectively. Left panel: Mass–central-density curves for different values of the EMSG coupling parameter α. Right panel: Corresponding sequences obtained by varying the quark-matter interaction parameter $ \bar{\lambda} $. Black dots denote the maximum-mass configurations along each sequence.

Figure 3 shows the radial profile of the charge density $ \rho_{\rm ch}(r) = \beta\,\rho(r) $ obtained from the numerical integration of the modified TOV equations. Each curve is evaluated at the central density corresponding to the maximum-mass configuration for that parameter set, as listed in Tables 1 and 2. In both panels, $ \rho_{\rm ch} $ reaches its peak value at the stellar centre, $ \rho_{\rm ch}(0) = \beta\,\rho_c $, where $ \rho_c $ differs for each curve according to the maximum-mass entry in the corresponding table. It then decreases monotonically outward while remaining strictly positive throughout the interior. The profile is finite and smooth at $ r = 0 $, consistent with the regularity condition imposed by the boundary condition $ q(0) = 0 $. At the stellar surface, $ \rho_{\rm ch} $ does not vanish but approaches a nonzero value determined by the finite surface density of quark matter. This is a characteristic feature of the bag-model equation of state, for which the pressure vanishes at a nonzero energy density set by $ B_{\rm eff} $. In the left panel, more negative values of α correspond to higher maximum-mass central densities and therefore produce larger central charge densities, while also extending the star to slightly larger radii. In the right panel, smaller values of $ \bar{\lambda} $ are associated with softer equations of state and higher maximum-mass central densities, yielding larger central charge densities but shorter stellar radii. In all cases, the charge density remains positive definite, confirming that the adopted ansatz $ \rho_{\rm ch} = \beta\rho $ with $ \beta \gt 0 $ introduces a physically consistent, positive-definite charge distribution throughout the stellar interior.

Figure 3.  (color online) Gravitational mass M as a function of the central energy density $ \rho_c $ for charged quark star configurations in EMSG. The parameter spaces used in the left and right panels are identical to those adopted in Fig. 1 and Fig. 2, respectively. Left panel: Mass–central-density curves for different values of the EMSG coupling parameter α. Right panel: Corresponding sequences obtained by varying the quark-matter interaction parameter $ \bar{\lambda} $. Black dots denote the maximum-mass configurations along each sequence.

The equilibrium stability of charged quark stars in EMSG is closely linked to the topology of the mass–central-density relation, as shown in Fig. 4. For each parameter set, the maximum mass marks the onset of radial instability, consistent with the classical stability argument developed by Harrison et al. [64] and further established by Zeldovich and Novikov [65]. Along the $ M(\rho_c) $ curve, models on the rising branch (${\rm d}M/{\rm d}\rho_c \gt 0$) are dynamically stable, whereas those on the descending branch (${\rm d}M/{\rm d}\rho_c \lt 0$) correspond to unstable configurations.

Figure 4.  (color online) Radial profiles of the charge density $ \rho_{\rm ch}(r) = \beta\,\rho(r) $ for charged quark star configurations in EMSG, each evaluated at the central density of the maximum-mass star for the corresponding parameter set. Left panel: Dependence of $ \rho_{\rm ch}(r) $ on the EMSG coupling parameter α, with central densities taken from Table 1 ($ \bar{\lambda}=0.1 $, $ \beta=0.5 $, $ B_{\rm eff}=60\; \text{MeV fm}^{-3} $). Right panel: Corresponding profiles for varying interaction parameter $ \bar{\lambda} $, with central densities from Table 2 ($ \alpha=+1.0\,\mu $, $ \beta=0.5 $). In all cases, $ \rho_{\rm ch} $ is finite and smooth at the stellar centre, remains strictly positive throughout the interior, and decreases monotonically to a nonzero surface value characteristic of quark matter with a finite bag constant.

Variations in the EMSG coupling parameter α and the matter interaction strength $ \bar{\lambda} $ systematically shift the location and extent of the stable sequence. Strongly negative α values or enhanced interaction parameters yield higher maximum masses and extend the stable region toward lower central densities; conversely, positive α values or weaker interaction strengths reduce the stability domain. This behavior reflects both the influence of nonlinear gravity and the underlying quark matter microphysics, reaffirming the utility of the turning-point criterion, as indicated in Refs. [64, 65], for diagnosing stellar stability in highly compact, strongly interacting regimes.

The equilibrium stability of charged quark stars in EMSG is closely linked to the topology of the mass–central-density relation, as shown in Fig. 4. For each parameter set, the maximum mass marks the onset of radial instability, consistent with the classical stability argument developed by Harrison et al. [64] and further established by Zeldovich and Novikov [65]. Along the $ M(\rho_c) $ curve, models on the rising branch (${\rm d}M/{\rm d}\rho_c \gt 0$) are dynamically stable, whereas those on the descending branch (${\rm d}M/{\rm d}\rho_c \lt 0$) correspond to unstable configurations.

Figure 4.  (color online) Radial profiles of the charge density $ \rho_{\rm ch}(r) = \beta\,\rho(r) $ for charged quark star configurations in EMSG, each evaluated at the central density of the maximum-mass star for the corresponding parameter set. Left panel: Dependence of $ \rho_{\rm ch}(r) $ on the EMSG coupling parameter α, with central densities taken from Table 1 ($ \bar{\lambda}=0.1 $, $ \beta=0.5 $, $ B_{\rm eff}=60\; \text{MeV fm}^{-3} $). Right panel: Corresponding profiles for varying interaction parameter $ \bar{\lambda} $, with central densities from Table 2 ($ \alpha=+1.0\,\mu $, $ \beta=0.5 $). In all cases, $ \rho_{\rm ch} $ is finite and smooth at the stellar centre, remains strictly positive throughout the interior, and decreases monotonically to a nonzero surface value characteristic of quark matter with a finite bag constant.

Variations in the EMSG coupling parameter α and the matter interaction strength $ \bar{\lambda} $ systematically shift the location and extent of the stable sequence. Strongly negative α values or enhanced interaction parameters yield higher maximum masses and extend the stable region toward lower central densities; conversely, positive α values or weaker interaction strengths reduce the stability domain. This behavior reflects both the influence of nonlinear gravity and the underlying quark matter microphysics, reaffirming the utility of the turning-point criterion, as indicated in Refs. [64, 65], for diagnosing stellar stability in highly compact, strongly interacting regimes.

The adiabatic index,

The adiabatic index plays a fundamental role in diagnosing the dynamical stability of relativistic stars against small radial perturbations. According to the classical analysis by Chandrasekhar [66], a stellar configuration becomes dynamically unstable when the effective stiffness of matter falls below the critical threshold $ \gamma \lt 4/3 $, a limit characteristic of relativistically supported compact objects. Later refinements have emphasized that the precise stability boundary depends sensitively on the underlying microphysical equation of state and the strong-gravity environment [67]. In modified gravity theories in which energy-momentum conservation may be violated, additional corrections can also affect the radial stability condition [68].

Figure 5 presents the radial behavior of the adiabatic index for charged quark star configurations in the EMSG framework, obtained using the same parameter sets employed in Figs. 1 and 2. In all cases, the profiles increase monotonically from the core toward the stellar surface and consistently remain above the relativistic instability threshold $ \gamma = 4/3 $. This trend reflects the progressive stiffening of quark matter at lower densities and indicates that the equilibrium sequences satisfy the necessary condition $ \gamma \gt 4/3 $ for stability against radial perturbations.

Figure 5.  (color online) Radial profiles of the adiabatic index $ \gamma(r) $ for charged quark star configurations in energy-momentum squared gravity. The left panel shows how γ varies as a function of radius for representative values of the EMSG coupling parameter α, based on the parameter space adopted in Fig. 1. In the right panel, the adiabatic index is plotted as a function of the quark matter interaction parameter $ \bar{\lambda} $, with parameter choices matching those in Fig. 2. In all cases, γ rises steeply near the stellar surface and consistently remains above the relativistic instability limit ($ \gamma = 4/3 $, dash-dotted line), {so that the explored models satisfy this necessary condition for stability against radial perturbations even in regimes of strong gravity and quark interactions.

Variations in the EMSG coupling parameter α or in the quark-matter interaction strength $ \bar{\lambda} $ primarily affect the outer layers of the star, with more negative α or larger $ \bar{\lambda} $ yielding slightly higher values of γ near the surface. Nevertheless, no parameter choice within the explored domain produces regions with $ \gamma \lt 4/3 $, so the necessary condition $ \gamma \gt 4/3 $ for radial stability is satisfied throughout the stellar interiors. This behavior is fully consistent with earlier analyses of quark stars in the EMSG framework [34, 69]. It demonstrates that the nonlinear matter-geometry coupling does not introduce destabilizing effects for the parameter ranges considered here. We emphasize that all stability diagnostics presented here—the adiabatic index, the sound speed, and the turning-point criterion—are evaluated for the total effective fluid, which incorporates both the quark matter contributions and the EMSG gravity correction terms, following the approach adopted in Refs. [70, 71].

The adiabatic index,

The adiabatic index plays a fundamental role in diagnosing the dynamical stability of relativistic stars against small radial perturbations. According to the classical analysis by Chandrasekhar [66], a stellar configuration becomes dynamically unstable when the effective stiffness of matter falls below the critical threshold $ \gamma \lt 4/3 $, a limit characteristic of relativistically supported compact objects. Later refinements have emphasized that the precise stability boundary depends sensitively on the underlying microphysical equation of state and the strong-gravity environment [67]. In modified gravity theories in which energy-momentum conservation may be violated, additional corrections can also affect the radial stability condition [68].

Figure 5 presents the radial behavior of the adiabatic index for charged quark star configurations in the EMSG framework, obtained using the same parameter sets employed in Figs. 1 and 2. In all cases, the profiles increase monotonically from the core toward the stellar surface and consistently remain above the relativistic instability threshold $ \gamma = 4/3 $. This trend reflects the progressive stiffening of quark matter at lower densities and indicates that the equilibrium sequences satisfy the necessary condition $ \gamma \gt 4/3 $ for stability against radial perturbations.

Figure 5.  (color online) Radial profiles of the adiabatic index $ \gamma(r) $ for charged quark star configurations in energy-momentum squared gravity. The left panel shows how γ varies as a function of radius for representative values of the EMSG coupling parameter α, based on the parameter space adopted in Fig. 1. In the right panel, the adiabatic index is plotted as a function of the quark matter interaction parameter $ \bar{\lambda} $, with parameter choices matching those in Fig. 2. In all cases, γ rises steeply near the stellar surface and consistently remains above the relativistic instability limit ($ \gamma = 4/3 $, dash-dotted line), {so that the explored models satisfy this necessary condition for stability against radial perturbations even in regimes of strong gravity and quark interactions.

Variations in the EMSG coupling parameter α or in the quark-matter interaction strength $ \bar{\lambda} $ primarily affect the outer layers of the star, with more negative α or larger $ \bar{\lambda} $ yielding slightly higher values of γ near the surface. Nevertheless, no parameter choice within the explored domain produces regions with $ \gamma \lt 4/3 $, so the necessary condition $ \gamma \gt 4/3 $ for radial stability is satisfied throughout the stellar interiors. This behavior is fully consistent with earlier analyses of quark stars in the EMSG framework [34, 69]. It demonstrates that the nonlinear matter-geometry coupling does not introduce destabilizing effects for the parameter ranges considered here. We emphasize that all stability diagnostics presented here—the adiabatic index, the sound speed, and the turning-point criterion—are evaluated for the total effective fluid, which incorporates both the quark matter contributions and the EMSG gravity correction terms, following the approach adopted in Refs. [70, 71].

The squared sound speed,

The squared sound speed provides a direct measure of the stiffness of the equation of state and must satisfy the causality requirement $ v_s^{2} \leqslant 1 $ throughout the stellar interior. Figure 6 illustrates the radial behavior of $ v_s^{2} $ for the charged quark star models considered in this work.

Figure 6.  (color online) Radial profiles of the squared sound speed $ v_s^{2} = dp/d\rho $ for charged quark star configurations in EMSG. Left panel: Dependence of $ v_s^{2}(r) $ on different values of the matter-geometry coupling parameter α, using the same parameter space as in Fig. 1. Right panel: Corresponding sound-speed profiles for varying interaction strength $ \bar{\lambda} $, with parameter choices matching those of Fig. 2. In all cases, $ v_s^{2} $ increases gradually from the stellar core toward the surface and remains well below the causal limit $ v_s^{2} = 1 $, confirming that the explored models satisfy the causality condition throughout the interior.

Across all parameter sets, the sound speed increases slightly from the stellar core toward the outer layers and remains consistently below the causal limit. Variations in the EMSG coupling parameter α (left panel) produce only minimal changes in the profiles, indicating that nonlinear matter-geometry effects do not significantly modify the microphysical stiffness of matter. Changes in the interaction parameter $ \bar{\lambda} $ (right panel) generate a slightly wider spread in the curves, with larger $ \bar{\lambda} $ corresponding to a stiffer response while remaining safely within the subluminal regime.

These results demonstrate that all configurations examined—regardless of α, $ \bar{\lambda} $, or the presence of electric charge—satisfy the causality condition, confirming the physical viability of the charged quark star solutions obtained within the EMSG framework. We emphasize that the stability assessment presented in this work is based exclusively on static criteria: the turning-point condition applied to the $ M(\rho_c) $ relation, the adiabatic index condition $ \gamma \gt 4/3 $, and the causality requirement $ v_s^2 \leqslant 1 $. A rigorous determination of dynamical stability would require a full radial oscillation (eigenfrequency) analysis—solving the radial pulsation equations for the normal-mode spectrum $ \omega_n^2 $—that consistently incorporates the EMSG gravity corrections, the electromagnetic field, and the quark matter interactions. Such an analysis represents a technically demanding extension that lies outside the scope of the present study and is left for future work [68].

The squared sound speed,

The squared sound speed provides a direct measure of the stiffness of the equation of state and must satisfy the causality requirement $ v_s^{2} \leqslant 1 $ throughout the stellar interior. Figure 6 illustrates the radial behavior of $ v_s^{2} $ for the charged quark star models considered in this work.

Figure 6.  (color online) Radial profiles of the squared sound speed $ v_s^{2} = dp/d\rho $ for charged quark star configurations in EMSG. Left panel: Dependence of $ v_s^{2}(r) $ on different values of the matter-geometry coupling parameter α, using the same parameter space as in Fig. 1. Right panel: Corresponding sound-speed profiles for varying interaction strength $ \bar{\lambda} $, with parameter choices matching those of Fig. 2. In all cases, $ v_s^{2} $ increases gradually from the stellar core toward the surface and remains well below the causal limit $ v_s^{2} = 1 $, confirming that the explored models satisfy the causality condition throughout the interior.

Across all parameter sets, the sound speed increases slightly from the stellar core toward the outer layers and remains consistently below the causal limit. Variations in the EMSG coupling parameter α (left panel) produce only minimal changes in the profiles, indicating that nonlinear matter-geometry effects do not significantly modify the microphysical stiffness of matter. Changes in the interaction parameter $ \bar{\lambda} $ (right panel) generate a slightly wider spread in the curves, with larger $ \bar{\lambda} $ corresponding to a stiffer response while remaining safely within the subluminal regime.

These results demonstrate that all configurations examined—regardless of α, $ \bar{\lambda} $, or the presence of electric charge—satisfy the causality condition, confirming the physical viability of the charged quark star solutions obtained within the EMSG framework. We emphasize that the stability assessment presented in this work is based exclusively on static criteria: the turning-point condition applied to the $ M(\rho_c) $ relation, the adiabatic index condition $ \gamma \gt 4/3 $, and the causality requirement $ v_s^2 \leqslant 1 $. A rigorous determination of dynamical stability would require a full radial oscillation (eigenfrequency) analysis—solving the radial pulsation equations for the normal-mode spectrum $ \omega_n^2 $—that consistently incorporates the EMSG gravity corrections, the electromagnetic field, and the quark matter interactions. Such an analysis represents a technically demanding extension that lies outside the scope of the present study and is left for future work [68].

In this work, we conducted a detailed investigation of electrically charged quark stars within the framework of energy-momentum squared gravity (EMSG), employing a unified interacting quark-matter equation of state that incorporates perturbative QCD effects, color superconductivity, and finite strange-quark-mass contributions. By extending the Maxwell-EMSG field equations to include nonlinear matter-geometry couplings and a density-dependent charge profile, we derived the corresponding stellar-structure equations and systematically analyzed the model-parameter dependence of equilibrium configurations. Our study demonstrates that the combined effects of electric charge and EMSG corrections can significantly influence the global properties of ultra-dense compact stars.

The mass-radius relations show that negative values of the EMSG coupling parameter α substantially enhance the pressure support generated by the quadratic matter term, enabling both neutral and charged stars to sustain larger gravitational masses at lower central densities. When a moderate electric charge is included, this stabilizing effect becomes even more pronounced, allowing maximum masses approaching $ 3M_{\odot} $ with compactness values near $ M/R \simeq 0.34 $. These results place charged quark stars in EMSG within the observationally relevant regime of massive pulsars such as PSR J0952-0607 {and, under a possible compact-star interpretation, the mass range associated with the secondary component of GW190814}. Variations in the interaction parameter $ \bar{\lambda} $ further broaden the accessible mass range, confirming that quark-matter microphysics and nonlinear gravitational effects jointly determine the structure of dense stellar objects.

We also examined the stability of the resulting configurations using several complementary diagnostics. The turning-point criterion applied to the $ M(\rho_{c}) $ relation indicates that the onset of instability occurs at the maximum-mass point for each sequence, consistent with the classical analyses of Chandrasekhar and subsequent refinements. Moreover, the radial profiles of the adiabatic index $ \gamma(r) $ remain above the relativistic instability threshold $ \gamma = 4/3 $ throughout the interior, regardless of the choice of α or $ \bar{\lambda} $. This behaviour is further supported by the sound-speed analysis, which shows that $ v_{s}^{2} $ remains well below the causal limit across the entire stellar radius. These complementary indicators — the turning-point criterion, $ \gamma \gt 4/3 $, and $ v_s^2 \leqslant 1 $ — are necessary conditions for stability and are satisfied throughout the explored parameter space, indicating physically viable charged quark-star configurations within the EMSG framework; a rigorous assessment of dynamical stability through a full radial oscillation (eigenfrequency) analysis is left for future work.

The broader implications of our results emphasize the importance of nonlinear matter-geometry couplings in shaping the behaviour of strongly interacting systems. The ability of EMSG to support high-mass compact stars without violating causality or stability constraints provides a compelling extension of general relativity in the strong-field regime. At the same time, the observed sensitivity of stellar properties to the interaction strength $ \bar{\lambda} $ underscores the need for consistent, phenomenologically motivated quark-matter models when interpreting astrophysical observations. The unified interacting EoS employed here provides a versatile framework for studying matter under extreme conditions while remaining consistent with current multimessenger constraints.

Several avenues for future research naturally emerge from the present investigation. A detailed analysis of tidal deformability, particularly in the context of binary mergers observed by LIGO–Virgo–KAGRA, would help constrain the viable parameter ranges for α, β, and $ \bar{\lambda} $. Extensions to slowly or rapidly rotating configurations, universal relations for moment-of-inertia compactness, and the potential role of anisotropic pressures would also provide deeper insight into the phenomenology of quark stars in EMSG. Additionally, alternative charge distributions or more general matter Lagrangians could further clarify the interplay between nonlinear gravitational effects and electromagnetic contributions. These efforts would help determine whether EMSG leaves distinctive observational signatures that could distinguish it from general relativity using current and forthcoming astrophysical measurements.

In this work, we conducted a detailed investigation of electrically charged quark stars within the framework of energy-momentum squared gravity (EMSG), employing a unified interacting quark-matter equation of state that incorporates perturbative QCD effects, color superconductivity, and finite strange-quark-mass contributions. By extending the Maxwell-EMSG field equations to include nonlinear matter-geometry couplings and a density-dependent charge profile, we derived the corresponding stellar-structure equations and systematically analyzed the model-parameter dependence of equilibrium configurations. Our study demonstrates that the combined effects of electric charge and EMSG corrections can significantly influence the global properties of ultra-dense compact stars.

The mass-radius relations show that negative values of the EMSG coupling parameter α substantially enhance the pressure support generated by the quadratic matter term, enabling both neutral and charged stars to sustain larger gravitational masses at lower central densities. When a moderate electric charge is included, this stabilizing effect becomes even more pronounced, allowing maximum masses approaching $ 3M_{\odot} $ with compactness values near $ M/R \simeq 0.34 $. These results place charged quark stars in EMSG within the observationally relevant regime of massive pulsars such as PSR J0952-0607 {and, under a possible compact-star interpretation, the mass range associated with the secondary component of GW190814}. Variations in the interaction parameter $ \bar{\lambda} $ further broaden the accessible mass range, confirming that quark-matter microphysics and nonlinear gravitational effects jointly determine the structure of dense stellar objects.

We also examined the stability of the resulting configurations using several complementary diagnostics. The turning-point criterion applied to the $ M(\rho_{c}) $ relation indicates that the onset of instability occurs at the maximum-mass point for each sequence, consistent with the classical analyses of Chandrasekhar and subsequent refinements. Moreover, the radial profiles of the adiabatic index $ \gamma(r) $ remain above the relativistic instability threshold $ \gamma = 4/3 $ throughout the interior, regardless of the choice of α or $ \bar{\lambda} $. This behaviour is further supported by the sound-speed analysis, which shows that $ v_{s}^{2} $ remains well below the causal limit across the entire stellar radius. These complementary indicators — the turning-point criterion, $ \gamma \gt 4/3 $, and $ v_s^2 \leqslant 1 $ — are necessary conditions for stability and are satisfied throughout the explored parameter space, indicating physically viable charged quark-star configurations within the EMSG framework; a rigorous assessment of dynamical stability through a full radial oscillation (eigenfrequency) analysis is left for future work.

The broader implications of our results emphasize the importance of nonlinear matter-geometry couplings in shaping the behaviour of strongly interacting systems. The ability of EMSG to support high-mass compact stars without violating causality or stability constraints provides a compelling extension of general relativity in the strong-field regime. At the same time, the observed sensitivity of stellar properties to the interaction strength $ \bar{\lambda} $ underscores the need for consistent, phenomenologically motivated quark-matter models when interpreting astrophysical observations. The unified interacting EoS employed here provides a versatile framework for studying matter under extreme conditions while remaining consistent with current multimessenger constraints.

Several avenues for future research naturally emerge from the present investigation. A detailed analysis of tidal deformability, particularly in the context of binary mergers observed by LIGO–Virgo–KAGRA, would help constrain the viable parameter ranges for α, β, and $ \bar{\lambda} $. Extensions to slowly or rapidly rotating configurations, universal relations for moment-of-inertia compactness, and the potential role of anisotropic pressures would also provide deeper insight into the phenomenology of quark stars in EMSG. Additionally, alternative charge distributions or more general matter Lagrangians could further clarify the interplay between nonlinear gravitational effects and electromagnetic contributions. These efforts would help determine whether EMSG leaves distinctive observational signatures that could distinguish it from general relativity using current and forthcoming astrophysical measurements.

The authors sincerely thank the anonymous reviewers for their valuable comments and constructive suggestions, which have significantly improved the quality and clarity of the manuscript.

The authors sincerely thank the anonymous reviewers for their valuable comments and constructive suggestions, which have significantly improved the quality and clarity of the manuscript.