Massive neutron stars as mass gap candidates: Exploring equation of state and magnetic field

The densities in the cores of the neutron stars (NSs) can reach several times that of the nuclear saturation density. The exact nature of matter at these densities is still virtually unknown. We consider a number of proposed, phenomenological relativistic mean-field equations of state to construct theoretical models of NSs. We find that, based on our selected set of models, the emergence of exotic matter at these high densities restricts the mass of NSs to $\simeq 2.2 M_\odot$. However, the presence of magnetic fields and a model anisotropy significantly increases the star's mass, placing it within the observational mass gap that separates the heaviest NSs from the lightest black holes. Therefore, we propose that gravitational wave observations, like GW190814, and other potential candidates within this mass gap, may actually represent massive, magnetized NSs.


INTRODUCTION
Neutron stars (NSs) are the end products of the stellar evolution of main-sequence stars with masses between 10 and 25M ⊙ .NSs are some of the most extreme objects in the Universe.A typical NS contains a mass on the order of the Sun squeezed into a radius of ∼ 10 km.At their cores, NSs can have densities several times that of the nuclear saturation density [1,2].
However, there is still much that is unknown about NSs -chief among them being the equation(s) of state (EOS).The strong interactions occurring within the cores of NSs are poorly constrained and, in many cases, virtually unknown.This is further complicated by the fact that, at the exceptionally high densities within NS cores, "exotic" particles like hyperons and deltas become energetically favorable [3].Notably, these exotic particles are even less constrained than nucleons themselves and lead to their own theoretical challenges (for a summary of the "hyperon puzzle" and possible solutions, see [4]).As a result, unlike the celebrated Chandrasekhar limit of white dwarfs, NSs do not have a fixed maximum mass limit (although there have been some EOS-independent estimates -see [5][6][7] for instance).NSs are the only laboratory we know of at present where matter exists in such extreme conditions.As a result, the study of NSs becomes a direct probe of the rich theoretical physics of nuclear matter at high densities.
Theoretically, high-density nuclear matter EOS fall into two main categories -microscopic and phenomenological.The microscopic EOS start from bare two-and three-nucleon interactions to reproduce the nucleon scattering data and properties of bound systems with few nucleons.We will not be focusing on these in the current work.Instead, we explore the phenomenological EOS, which are effective interaction models.A class of phenomenological EOS that has been explored extensively in recent years is the relativistic mean field (RMF) models (for a recent review on NS EOS and associated physics, see [8]).These are constructed under the quantum hadrodynamics model, where the interactions in nuclear matter are modeled at the hadronic level, with baryon-baryon interactions taking place through the exchange of mesons.These EOS are constrained in two ways -using the properties of nuclear matter at saturation density, and the observations of NSs themselves.
On obtaining an appropriate EOS, one can construct a family of NSs from it, parametrized by the central density by solving the general relativistic hydrostatic equilibrium equations first proposed by Tolman, Oppenheimer, and Volkoff (TOV) [5].Thus, every EOS maps onto a massradius (M-R) curve, which is then used to calculate the theoretical mass limit of NSs.Thus, a range of nuclear matter EOS directly results in a range of mass limits.
From an observation perspective, NSs are challenging objects.Although the first radio pulsar observation was made in 1968 [9], NS measurements still harbor a lot of uncertainties.Particularly, the radius of a NS is very difficult to measure and is poorly constrained by the available data.Indeed, it is the measurements of masses, particularly those of heavy NSs, that serve as an important discriminator among the various EOS models.
Observations of "massive" NSs are, therefore, crucial to constrain the properties of high-density nuclear matter.The heaviest observed NS (as of the time of writing) is the "black widow" pulsar PSR J0952-0607 [10,11], which has a mass M = 2.35 ± 0.17M ⊙ .On the other hand, the lowest observed black hole mass appears to be around M = 3.3M ⊙ [12] (however, see also [13]).This has led to an apparent "mass gap" that exists between the heaviest NSs and the lightest black holes.
In particular, the GW190814 observation [16] revealed the merger of two compact objects, one of which was inferred to have a mass within the range of 2.50 − 2.67M ⊙ , squarely placing it within the mass gap.Arguments proposing that this event involved two black holes have been presented in [17].However, the true nature of this object, whether it is a black hole or a NS, remains uncertain, as no electromagnetic counterpart was identified for this event.Additional observations, such as GW200210-092254 [18], have also suggested the possibility of objects residing within the mass gap.
Apart from the possibility of mass gap detections, GW observations also significantly help with constraining NS radii.A parameter that has gained relevance in recent years is the dimensionless tidal deformability (Λ).This is related to the deformation of the star under external tidal fields, such as that of a companion object in a binary.This deformation leaves an imprint on the GW signal and hence, binary mergers such as GW170817 have been used to place an observational limit on this parameter.Although the limit varies depending on the priors and models considered, it proves to be a useful discriminant between various EOS candidates as it gives a range of acceptable stiffness for a particular EOS model.Since the tidal deformability is closely linked to the radius of the NS, this ends up indirectly providing an additional observational constraint.Combining GW observations with Neutron Star Interior Composition Explorer (NICER) data gives us a radius constraint of around 11 − 13 km for a 1.4M ⊙ NS [19,20].
As mentioned previously, the heaviest NS has a mass above 2M ⊙ .Other pulsar observations such as PSR J1614-2230, M = 1.97 ± 0.04M ⊙ [21,22]; PSR J0348+0432, M = 2.01 ± 0.04M ⊙ [23]; and PSR J0740+6620 M = 2.08 ± 0.07M ⊙ [24,25] support the idea that the NS's "Chandrasekhar" limit (if it exists) could be well above this value.On the other hand, it has been shown through previous work in our group (starting with [26]) that the magnetic fields of compact stars can lead to significant enhancements in their masses, both through classical and quantum effects.
NSs are known to have significant magnetic fields.Typical surface fields fall in the range of 10 8 − 10 13 G.An important subset of NSs is the strongly magnetized "magnetars," which can have surface fields up to the order of 10 15 G. Magnetars may account for around 10% of the NS population [27,28].Emphasizing once more, the field strengths mentioned here are all for the surface values.The fields at the center of the star can end up being orders of magnitude higher than the ones quoted here.
Hence, it seems that the magnetic fields of NSs and their associated effects on the M-R relationship are not ignorable.The introduction of magnetic fields can affect the star in two ways.Classically, the magnetic field can lead to a magnetic pressure that contributes to the hydrostatic balance of the star.Further, this magnetic pressure can also lead to an overall pressure anisotropy in the NS.Another important effect the magnetic field can have is on the quantum microstates of the nuclear matter itself, i.e., modification of the EOS through Landau quantization.However, as shown previously [29], this effect is only significant for fields above 3 × 10 18 G.In the present work, we restrict ourselves to fields well below this value and consider only the classical effects of the field.
Along with the anisotropy arising from the magnetic pressure, NSs can have other sources of being anisotropic stars.In a very basic sense, the matter distribution within the star itself can lead to an inherent anisotropy.This can arise due to a variety of physical effects.For instance, convective and/or turbulent mixing in the context of two fluid models can lead to anisotropy [30].At the high densities inside NSs, other effects can contribute to anisotropy.For instance, in the superdense NS cores, the formation of pion/kaon condensates ( [31,32]) is favored, leading to pressure reduction along the radial direction, but not overall.The presence of possible superfluidity ( [33,34]) in NS cores can be an additional source of anisotropy.
On introducing a magnetic field to the star, an additional source of uncertainty arises, which is the magnetic profile within the star.The profile of the magnetic field within the compact star is unknown, which means that one can explore a number of profiles, making sure they are consistent with the Einstein-Maxwell equations.One of the most widely used profiles is a density-dependent exponential profile for the magnetic field magnitude, first proposed by Bandyopadhyay et al. [35].Although this profile can be shown to be consistent with the Einstein-Maxwell equations (see the discussion in Sec 2.4 of [36], for example), some critics have argued that simulations from LORENE indicate the field profile should be more of a polynomial fit [37][38][39].However, one should note that LORENE and hence the resulting polynomial profile are only calculated and inferred in the context of purely poloidal configurations within the star.It has been well known that purely poloidal (or purely toroidal) field profiles in stars lead to instabilities [40,41], and rather mixed field configurations ensure the most stability, e.g.[42][43][44].Nevertheless, in this work, as a first approximation, we consider two main orientations for the magnetic field introduced -radially oriented (RO) fields are fields directed toward the radial direction, while transversely oriented (TO) fields are fields directed randomly in a direction perpendicular to the radial direction.The orientation of the magnetic field plays a crucial role.As described in this paper, our approximate choice of model profiles helps to explore a series of realistic EOS in contrast to the other simulations of our group where the field profile is more accurate, but the EOS is approximate [45,46].Further, there is no consensus on whether a magnetic field always leads to an increase in the mass of a compact star [47][48][49].Thus, the study of magnetized, anisotropic NSs can help resolve multiple open issues.
However, the TOV equations hold for isotropic, spher-ical stars and are modified in the presence of anisotropy, such as that introduced by magnetic fields.We follow the approach laid out by a previous paper in our group in the same line [36], and introduce models for the magnetic field and anisotropy constructed under the assumption of spherical symmetry.This assumption of approximate spherical symmetry is valid in a number of field configurations and magnitudes, as shown previously in our group through two-dimensional simulations [45].
In the present work, we look at the theoretical possibility of massive, magnetized NSs as possible mass gap candidates.The paper is structured as follows.In Sec. 2, we set up the basic formalism and equations solved to construct our NS models.We describe the EOS used, as well as the models introduced to describe the magnetic field and anisotropy present in the star.In Sec. 3, we describe our obtained results and discuss possible implications of the same.We explore a different magnetic field profile with its caveats in Sec. 4. We end by summarizing our results and with our conclusions in Sec. 5.

FORMALISM AND SET OF EQUATIONS
To describe NSs, we solve the general relativistic hydrostatic balance equations, i.e., the TOV equations.Because of the introduction of anisotropic effects (both from the magnetic field and general matter effects), the TOV equations are modified.We follow the same modification of the TOV equations outlined in previous work by our group [36], given by ) (for TO).
(2.2) Here, m denotes the mass, ρ the density, and B the magnitude of the magnetic field at a given radius r within the star.Because of the presence of anisotropy, the pressure along the radial direction, p r is different from the pressure along the transverse direction p t .This is captured in the effective anisotropy factor ∆ defined as ∆ = p t − p r + B 2 /4π (for RO) p t − p r − B 2 /8π (for TO). (2. 3) The magnetic field thus modifies the expressions for the hydrostatic equilibrium and ∆ in different ways based on its orientation (RO or TO).

Modified Bowers-Liang model for anisotropy
To close the system of equations defined above, we need a model functional form for ∆.We use the general parametric form first introduced by Bowers and Liang [50].As done previously in our group [36], we modify the Bowers-Liang form to further include the effects of the magnetic field.∆ is then given by This model is derived keeping in mind the following key assumptions - • The anisotropic force must vanish at the center, leading to the anisotropy term vanishing quadratically at the center.
• Anisotropy varies with position inside the star.
• ∆ includes the effects due to local fluid anisotropy as well as the anisotropy due to the magnetic field (both its magnitude and orientation).
Following previous work [50,51], we restrict κ to the range [−2/3, 2/3].This is to ensure the physicality of the solution and ensure that we do not get a positive dp/dr.
We further need to supplement the set of equations (2.1, 2.2, 2.4) with the EOS and magnetic field profile to make it completely solvable.

Magnetic field profile
We introduce a density-dependent magnetic field in the star [35], given by This profile gives us the magnitude of the field as a function of the density, and hence, the radius within the star.Here, B s corresponds to the surface field of the star, B 0 and ρ 0 control the field at the center, and η and γ are model parameters that control how the field decays from center to surface.
Throughout this work, we have chosen B s to be 10 15 G.However, the results we obtain are found to be largely independent of this parameter, as long as B s is not comparable to B 0 .

Equations of state
To describe the matter present within the star, we use a selection of phenomenological EOS, constructed using the RMF approach.Here, the matter is modeled at the hadron level, with the interactions between the baryons modeled using meson fields as mediators.The meson field strengths are then set to their mean values as per the RMF approximation.Three such meson fields are included in this work -the scalar meson σ, describing attraction between baryons; the vector meson ω, describing repulsion; and the isovector meson ρ, explaining isospin asymmetric interactions.
RMF EOS are constrained in two ways -they must be able to reproduce the observed NS properties, and they must reproduce the properties of symmetric nuclear matter (SNM) at saturation density (n 0 ).
At the high densities present in NS cores, the presence of exotic particles is energetically favorable.Exotic particles are those that do not exist in stable form under terrestrial conditions.In the present work, we have considered two such classes of exotic particles -hyperons (particles with at least one strange baryon in quark content) and ∆ particles (nonstrange baryons of spin 3/2).Thus, we explore pure nucleonic (npeµ) EOS along with hyperon and ∆ admixed (npeµ − Y ∆) EOS.
The general Lagrangian of the relativistic model in the mean field approximation is given by ( The baryon content (B) here comprises the nucleons The Lagrangian for interacting baryons is given by where g σB (n), g ωB (n), and g ρB (n) are the meson-baryon coupling constants, m B is the mass of the baryon, n is the baryon number density, τ = (τ 1 , τ 2 , τ 3 ) is the Pauli isospin matrix, and γ µ are the Dirac matrices.The coupling constants are set to reproduce the properties of SNM at n 0 within their experimental bounds.The inclusion of hyperonic matter has been done by including meson-hyperon couplings based on the SU(3) ESC08 model [52] and the inclusion of ∆ particles done by including a near-universal meson-∆ coupling : x i∆ = g i∆ /g iN = 1.2.This value of x i∆ is in accordance with the consideration of meson-∆ coupling based on group theory as explored in recent work [53].Here, g represents the coupling constants; the subscript ∆ indicates the ∆ particles, N the nucleons, and i represents the mesons mentioned above.
The Lagrangians for the leptons (∈ {e − , µ − }) and the mesons are given by (2.8) Additional nonlinear self-interaction terms (dependent on σ 3 and σ 4 ) are often introduced to ensure that the empirical values of properties like nuclear incompressibility (K 0 ) and effective nuclear mass (m * /m N ) are reproduced correctly.
Apart from introducing nonlinear terms to the Lagrangian, another approach that ensures the proper reproduction of all saturation properties is to make the coupling constants density dependent.The elimination of the need for nonlinear terms leads to density-dependent RMF parametrizations (DDRMF), introduced as where i ∈ {σ, ω, ρ}, x = n/n 0 .f (x) gives the form of the density dependence, typically described using an ansatz of the form for i ∈ {σ, ω}, and Thus, now, the parameters that must fit to the SNM properties at n 0 are the nine density-dependent parameters: a σ , b σ , c σ , d σ , a ω , b ω , c ω , d ω , and a ρ , along with the meson-baryon couplings at saturation density: g σB (n 0 ), g ωB (n 0 ), and g ρB (n 0 ).
Another density-dependent modification to the standard RMF scheme comes from the recently constrained slope of the symmetry energy at n 0 (L 0 ).Here, the isovector meson-baryon coupling constant alone is modified to take a density-dependent form as (2.13) Following previous work [54], we refer to this as the RMFL approach.The value of L 0 is then fixed based on the coefficient a ρ along with the standard meson-baryon coupling constants.
From the Lagrangian constructed using the RMF/DDRMF/RMFL approach, one can then obtain the baryon and meson field equations by applying the Euler-Lagrange method.As we are constructing RMF models, the meson fields are set to their mean values (values in the ground state).
Along with the meson mean field equations, we impose charge neutrality and baryon number conservation, leading to a system of five coupled non-linear equations to be solved.Additionally, for the RMFL and DDRMF models, baryon chemical equilibrium has to be imposed as the density-dependent coupling leads to a rearrangement energy contribution (for the detailed systems of equations, see [54]).The solution of this system of equation gives us the meson mean fields, ω, σ, and ρ, along with the Fermi momenta of the baryons under consideration [54].
We choose a few RMF EOS -GM1L [55], SWL [54], DD2 [56], DD-ME1 [57], DD-ME2 [58], and DDMEX [59]-that best satisfy both the constraints from SNM properties and the kind of astrophysical observations we seek to explain in this work.The first two EOS are constructed using the RMFL approach whereas the rest are DDRMF EOS.The SNM properties at n 0 for this list of EOS are shown in Table 1.Properties tabulated here are -energy per nucleon (E/N ), symmetry energy (J), slope of the symmetry energy (L 0 ), the curvature of the symmetry energy (J ′′ ), the incompressibility (K) and the effective mass (m * /m).
FIG. 1: Different RMF EOS: Here, the upper branch denotes the npeµ EOS, while the lower branch represents the hyperon-∆ admixed npeµ − Y ∆ EOS.
Figure 1 shows that the pure nucleonic EOS are much stiffer (i.e., have a higher pressure at a given energy density) when compared to the cases with exotic particles.This is consistent with past literature.Indeed, the softening of the EOS by hyperonic/exotic matter and the necessity to reconcile with massive NS observations has led to the hyperon puzzle, with different resolutions being proposed, including introducing further repulsive interactions.

Tidal deformability
We can further constrain the NS EOS using tidal deformability limits from GW observations.In the presence of an external gravitational field (ϵ ij ), a star develops a quadrupole moment (Q ij ) such that Q ij = −λϵ ij , where λ is the tidal deformability of the star.
Theoretically, one can link λ to the dimensionless second Love number k 2 , arising from gravitational multipole expansion, as where C=M/R is the compactness of the star.
We compute the tidal love number k 2 by solving the static, linearized perturbation equation arising from the external tidal field.This corresponds to a metric αβ + h αβ , where h αβ is the linearized, perturbation metric.By expanding h αβ into spherical harmonics, Y m l (θ, ϕ), and restricting to the l = 2, static, even-parity perturbations in the Regge-Wheeler gauge, we can write [60,61] where H 0 , H 2 , and K are all radial functions determined by the perturbed Einstein equations.Expanding the perturbed stress-energy tensor and subsequently inserting the fluid and metric perturbations in the linearized Einstein equations, we obtain H 0 = H 2 ≡ H and K ′ = Hν ′ + H ′ .Further subtracting the equation δG θ θ + δG ϕ ϕ = 16πδp from the tt component of the perturbed Einstein equations, we obtain a differential equation for H as Here, H(r) is a radial function arising from the static, linearized perturbations of the Einstein equations.As we are solving for k 2 , we restrict ourselves to the l = 2, static, even-parity perturbations of the perturbation metric.The other quantities are A = dp t /dp r , c 2 s = dp r /dρ (the speed of sound squared), e λ = [1 − 2m/r] −1 , and ν ′ = 2e λ (m + 4πp r r 3 )/r 2 .For isotropic stars, A = 1.
On obtaining H(r) by solving Eq. (2.15) simultaneously with our system of equations, one can compute the tidal Love number k 2 as In recent years, due to GW observations, observational constraints have been placed on the dimensionless tidal deformability.From the binary NS merger event GW170817, the dimensionless tidal deformability at 1.4M ⊙ has been constrained to be Λ 1.4 < 800 [14] and Λ 1.4 < 580 [15].The two different limits are due to different models.We have enforced both in our subsequent analyses.

RESULTS AND DISCUSSION
As mentioned previously, this work investigates the theoretical possibility of NSs with masses high enough to fall in the observational mass gap.

Isotropic cases: Pure EOS effect
To start with, we look at the completely isotropic case -magnetic field and anisotropy parameter κ are both set to zero.This gives us an idea of the role that pure EOS effects play in the M-R relationship.M-R curves for the npeµ and npeµ − Y ∆ EOS are displayed in Fig. 2. As expected, the stiffer npeµ EOS give higher masses, with the maximum mass limit (M max ) reaching as high as 2.57M ⊙ for the DDMEX EOS.The npeµ − Y ∆ EOS, being softer, give a mass limit M max of about 2.25M ⊙ (for the same DDMEX EOS).The results for the isotropic cases are shown in Table 2.However, on computing the tidal deformability (Fig. 3) of these cases, it is clear that the nucleon-only cases all violate at least one of the observational upper bounds, while the hyperon-∆ admixed cases do not.These results are consistent with NS properties obtained in previous work in this line using the same EOS [54,62,63].Nevertheless, as we will see in further results, the introduction of anisotropy and/or magnetic field helps to change the value of tidal deformability and even bring the nucleon-only cases below observational thresholds.However, based on energy considerations mentioned previously and the tidal deformability constraints, it appears that npeµ − Y ∆ cases are better favored overall.For the rest of this analysis, we consider only the npeµ − Y ∆ admixed cases.
FIG. 3: Tidal deformability Λ as a function of M for npeµ EOS (dotted) and npeµ − Y ∆ EOS (solid).The various curves are the same as Fig. 2.  The physical parameters of isotropic NSs, both by inclusion of exotic particles, and for pure nucleonic matter.
As seen from the isotropic results, pure EOS effects are not sufficient to bring the stable NS to mass gap levels (> 2.5M ⊙ ).As a result, we turn to the additional physics of magnetic fields and/or anisotropy.

Introducing magnetic field with fixed anisotropy parameter, κ
We first introduce magnetic field following Eq.(2.5) in the presence of a fixed anisotropy parameter, κ = 0.5.We start with a profile described by η = 0.2, γ = 2.The general shape of this profile within a star is given in Fig. 4. We vary B 0 to take values 1.2 × 10 18 and 0.9 × 10 18 G in the transverse direction and 0.6 × 10 18 and 0.9 × 10 18 G in the radial direction.B s is kept fixed at 10 15 G.As mentioned previously, varying B s (in a range not comparable to B 0 ) does not have an effect on the results.The trend obtained is similar to previous work done in this line [36].Introducing anisotropy enhances the mass of the star, even in the absence of magnetic fields.For the latter, it is the matter anisotropy that is ultimately contributing to the enhanced mass of the star.However, it is important to note here that in general anisotropy is not wholly independent of the magnetic field and its effect.On introducing a TO field to an already anisotropic (by matter) star, its mass further increases, whereas a RO field tends to decrease the mass.In these cases, the anisotropy is not wholly a matter effect, but instead has contributions from the magnetic field introduced.This leads to the anisotropy being a consequence and/or an extension of the magnetic effect on the NS's structure and properties.Thus, different M-R curves are obtained on varying the magnetic field.The results for the two stiffest EOS (DD-ME2 and DDMEX) are shown in Figs. 5 and 6. Results for all the EOS are given in Table 3.
It seems that the highly magnetized TO stars are most promising in terms of mass gap candidates.As seen from Table 3, all the EOS give stars of maximum mass > 2.5M ⊙ for the TO field with B 0 = 1.2 × 10 18 G.The mass of the NS is increased up to even 2.8M ⊙ for the DDMEX EOS.However, on computing the tidal deformability for these cases (Figs. 7 and 8), we see that   these highly magnetized cases fail to meet either of the observational bounds set by GW observations.Similar issues arise when we perform a stability analysis [42], as all the B 0 = 1.2×10 18 G cases give high magnetic energy (E mag ) to gravitational energy (E grav ) ratios of order 0.1 (not explicitly shown here).Can we then rule out this much high fields at the cores of NSs completely?The answer is no.It is important to note that the results discussed so far are profile specifici.e., specific to a particular η and γ.We next try a second profile: η = 0.01, γ = 2.The general shape of this profile within the star is shown in Fig. 9. On comparing with the previous profile, we see that the magnetic field here has a much more gradual fall in the star, leading to a lower Lorentz force near the surface of the star.Hence, we expect that this profile will lead to lower masses, however, with better prospects for stability.As seen from the results, we require central fields (B c ) of around 5 × 10 17 G or higher for the magnetic field to significantly affect the mass of the star.For comparable values of B c in both chosen profiles, we see that the second profile leads to lower masses.For DDMEX, for instance, the second profile gives a maximum mass of 2.69M ⊙ as opposed to 2.84M ⊙ in the first profile.However, E mag /E grav is now reduced from nearly 0.2 for the first profile to 0.05 for the second.Similar trends are seen for the other EOS used.Thus, this second profile gives us much more stable stars overall.
We further look at the tidal deformability computed for the second profile.The results are shown in Fig. 11.We see that the lower η leads to much lower Λ 1.4 .In fact, most of the values of Λ 1.4 are clustered around 580.This means that these stars are better candidates to satisfy even the stricter observational limit on Λ 1.4 .
Although there was no measurable tidal deformations in the GW190814 signal, there have been tidal deformability constraints from this event, obtained from reexamining the analysis of GW170817 under the assumption of GW190814 having a massive NS.In particular, the spectral EOS distribution from GW170817 was taken and each EOS was reweighted by the probability that its maximum mass falls in the mass of the secondary object detected in the GW190814 [16].This favors stiffer EOS with the updated tidal deformability constraint being 458 < Λ 1.4 < 889.These updated limits are additionally shown in Fig. 11.We see that all the cases, including the highly magnetized ones, are consistent with all the tidal deformability constraints laid out so far.
Exploring just these two profile cases, we thus see a wide variety of behavior possible by experimenting with the model parameters of the magnetic field.Since the true form of the magnetic field within the star is essentially unknown, it is important to experiment with different profiles and check the results.In particular, we cannot rule out the possibility of high fields existing within NSs based on results from one profile alone.
We see that higher values of κ lead to higher masses being supported by the star, even at zero field.From the results tabulated in Table 5, we see that for a given field, a higher κ leads to a lower E mag /E grav (and lower B c ) along with a higher mass.Higher κ also reduces the tidal deformability (Fig. 14).EOS that violate the observational bounds on Λ in the case of isotropic cases can be seen to be consistent with the same bounds when we introduce anisotropy to the star.This dependence of Λ on the degree of anisotropy present in the star could be crucial in determining the observational bounds on κ.Indeed, previous work [61] discussed how κ could be eventually constrained using GW observations.However, as of now, there are not any observational bounds on κ.Nevertheless, the presence of anisotropy within a NS is well motivated, both by the possible presence of a magnetic field and by high-density effects, such as superfluidity.Previous studies also explored the possibility of the mass gap object in GW190814 being an anisotropic NS [64].From top to bottom, the curves for each EOS correspond to κ = 0.6, 0.5, 0.4, 0.3, 0.2, and 0.1 sequentially.

Universal relations
As mentioned throughout this paper, the primary challenge arising in the theoretical study of NSs is the unknown EOS.However, unlike the M-R curves explored so far, there are "universal relations" between NS parameters that are independent/insensitive to the underlying EOS.Such relations are found to exist between the moment of inertia (I) of the star, stellar compactness (C), and the tidal deformability (Love).This leads to the B0(G) κ Bc(10   They can also be a source of EOS-independent tests of the fundamental physics of general relativity [65,66].These universal relations were first introduced/studied in the context of isotropic NSs constructed under the slow rotation approximation.Recent publications [61,67,68] have extended this to anisotropic stars as well.We now examine if this universality can be extended to the magnetic, anisotropic NSs we have discussed so far in this paper.
In Figs.
15, 16 and 17, the C − I, C − Λ and I − Λ relations are shown for the DDMEX EOS by varying κ considering its different values as 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6.As seen from the figures, the C −I relations are unaffected by the presence of anisotropy in the star.However, the C −Λ and I −Λ relations show slight deviation once anisotropy is introduced, particularly at the low Λ regime.This slight anisotropic deviation was also seen in previous work [67].To find the extent of the deviation, we fit the relations using the least squares fitting.
where c n represents the best-fit parameters.The calculated values of parameters a n , b n , and c n , along with the standard error (square root of variance) on each parameter are given in Table 6.To quantify the goodness of fit, we calculate the root-mean-squared error (RMSE), given by the square root of the mean of the difference between the predicted and actual results squared.The corresponding RMSE for the varying κ cases is given in Table 7.
To examine the effect of magnetic field on the universal relations, we look at the same three relations -C − I, C − Λ, and I − Λ, but by fixing κ to be 0.5.Along with the nonmagnetized cases, we also include B 0 = 2 × 10 18 and B 0 = 5 × 10 18 G cases.The results for the DDMEX EOS are shown in Figs. 18, 19, and 20.We see that the universal relations are unaffected by the presence of the magnetic field.
The fitting is done as explained previously.The bestfit parameters, along with the errors on the same, are given in Table 8.The RMSE is also listed in Table 9.       the fits are shown in Table 11.on the universality of these relations, more precisely on C − Λ and I − Λ relations.Varying κ makes the RMSE of the fit roughly 2-4 times that of the case with fixed κ.However, varying magnetic field while keeping κ fixed gives, indeed, a better universal fit.As explained in previous work [61], the rigorous exploration of the effect of anisotropy on the I-Love-C relations could lead to possible observational constraints on κ and the degree of anisotropy in the star in general.

COMPARISON BETWEEN POLYNOMIAL AND EXPONENTIAL PROFILES
All the above explorations are based on a particular law of the density-dependent magnetic field profile, which offers an exponential variation.One may question if the density-dependent magnetic field profile used in this work (Eq.2.5) is inconsistent with the Maxwell equations [37].Actually, this should not be an issue.Under the formalism of approximate spherical symmetry (which is quite a good approximation, particularly, for toroidal fields) described in this paper, the density-dependent profile only gives us the magnitude of the magnetic field at a particular density (and hence, radius) within the star.The orientation is determined by whether we consider "RO" or "TO" fields.As we have shown in the previous work in the same line, we can always make the magnitude obtained from the density-dependent profile consistent with the Maxwell equations, approximately in the regime of spherical symmetry, for both of the orientations considered in this work.A detailed proof of this is shown in Sec.2.4 of [36], under similar formalism and assumptions.
One of the alternative profiles is a polynomial profile, with magnetic field expressed as a function of the baryon chemical potential (µ B ).The profile is obtained by fitting the two-dimensional results from the numerical code LORENE using quadratic polynomials [38,39], given by We introduce magnetic fields to our star using the above profile by varying µ with values 2 × 10 31 , 5 × 10 31 , 10 32 , and 2×10 32 Am 2 .An important caveat, however, here is that LORENE can only generate poloidal magnetic fields.In the current formalism, this translates to an "RO" field.We set the anisotropy κ = 0.1 throughout.
The M-R relations for the four magnetized cases mentioned above, along with a nonmagnetized (but still anisotropic) case for the DDMEX EOS are shown in Fig. 24.The corresponding values of the tidal deformability are shown in Fig. 25.As this profile leads to pure poloidal/RO fields, we see that under the current formalism, it leads to an overall decrease in the mass of the NS.
It appears that the magnetic field effect is extremely limited in the case of this profile as opposed to the previous profile.From the results listed in  see that the NS properties are not all that different from the nonmagnetized cases even for high magnitudes of the central (maximum) magnetic field and E mag /E grav .The trend of variation of the magnetic field within the star for each of the µ cases is shown in Fig. 26.We see that this profile always gives, at most, an order-of-magnitude difference between the center and the surface fields.This means that a maximum magnetic field of order 10 18 G within the star corresponds to a surface field of similar order or, at the maximum, reduced up to 10 17 G.This is much higher than the present measurements for the surface magnetic field of NSs.Additionally, such a high field throughout the star leads to high E mag /E grav (Table 12), which might destabilize the star.On the other hand, if we aim to achieve the surface field ∼ 10  Indeed, this 1 order-of-magnitude variation of the magnetic field from the center to surface is consistent with previous work, done using this profile [38,39,72].It is also consistent with the numerical simulation results for the poloidal fields from our group using the XNS code [45,46].However, it has been long known [40,41] that purely poloidal or purely toroidal field configurations lead to unstable stars.The most stable field configuration in NSs is a mixed field configuration consisting of both poloidal and toroidal fields; perhaps toroidally dominated [73], also see [44,74].Hence, it would be incorrect to rule out the magnetic field's effect on the star based on the results from the purely poloidal polynomial profile.In fact, XNS code shows that toroidal fields can have 2 orders of magnitude variation between the maximum field and the outer crust field (exactly at the surface it exhibits zero or very low field by the very geometry).Hence, the chosen field profile with the exponential variation, as used in all the previous sections, is likely to be a natural situation with the right combination of toroidal and poloidal fields within the star.The toroidal field varies from its maximum to minimum (close to the surface) by 2 orders of magnitude and the corresponding poloidal field varies by one order.If the core field is toroidally dominated with 2 orders of magnitude higher than its poloidal counterpart, then for a maximum toroidal field of 10 18 G, the surface poloidal field could be 10 15 G (exactly at the surface, the toroidal field vanishes).In such a star, the poloidal field hardly affects the stellar structure.This further justifies the approximate spherical choice of our model, when a significant deviation from the spherical symmetry is only by the poloidal field.
Thus, although the polynomial profile is constructed in order to be consistent with the Einstein-Maxwell equations, it does not account for the instabilities that arise from purely poloidal configurations as well as possible instability due to high E mag /E gav .Indeed, all the profilebased approaches, whether they be the exponential profile or the polynomial profile, are ultimately approximations that have limitations in predicting the NS physics.For a first-order study into the NS's mass limits, such as the one done in this work, it appears that the exponential profile is more than adequate.For more sophisticated studies of the magnetic field profiles and the precise effect of its geometry, one must turn to two-dimensional Einstein-Maxwell solvers such as XNS, as done already by some of us [45,46].However, at present, the XNS code cannot deal with realistic EOS such as the ones studied in this paper, it can deal only with the polytropic form.Thus the approximate studies, as reported in the present paper, become important to gain insights into the EOS dependence of NS properties.

CONCLUSIONS
In this work, we investigate the possibility of massive NSs being mass gap candidates to explain observations such as GW190814.We first examine the pure EOS effect.Because of the appearance of exotic particles such as hyperons, the EOS is softened considerably as compared to the pure nucleonic case.The maximum mass obtained for isotropic stars constructed from our chosen set of npeµ−Y ∆ EOS was only ≃ 2.2M ⊙ .To truly bring our theoretical models to mass gap ranges, we next investigate the effects of magnetic field and a model anisotropy in the NS.Crucially, both magnetic field and anisotropy are well-established physical effects that show up in NSs.The observations of magnetars, for example, support NSs having surface fields as high as 10 15 G. Anisotropy is also expected to be present, both as a consequence of the magnetic field introduced, and as an independent effect arising from high-density effects such as superfluidity.Introducing anisotropy enhances the mass of the star, even in the absence of the magnetic field.Depending on the orientation of the magnetic field introduced (RO or TO), the maximum mass can be either enhanced or decreased from the nonmagnetized value.Only the toroidal field leads to a reasonable increment in the mass of a NS.Although we introduce these two nonspherical effects to the system, we construct our NS models in approximate spherical symmetry.This is a reasonable approximation, particularly for toroidally dominated stars, as shown previously by our group using the two-dimensional numerical simulation code XNS [45].
Nevertheless, many authors (e.g.[75][76][77]) tackle the EOS based restriction to the NS mass through the introduction of additional effects such as rotation in the star.There have also been arguments based on GW170817 observations combined with universal relations that the maximum mass of non-rotating, nonmagnetized NSs is within 2.3M ⊙ [78]; see also [79].However, the hyperonization/softening issue may be confronted through the adjustment of the coupling constant as well [80], which though appears to be an extreme assumption [17].Moreover, the possibility of the ligher component of GW190814 to be a quark star is not ruled out [81].Therefore, while the alternative choices may lead to slightly varied outcomes, the fundamental findings of the paper persist in their robustness.
We further use Λ and E mag /E grav to ensure the physicality of our results.We find that high magnetic fields from certain profiles are ruled out under these constraints.However, with the right choice of profile, NSs of masses in the range 2.5 − 2.67M ⊙ are demonstrated to be possible, while satisfying all stability criteria.On examining the role of anisotropy in the star, we find that increasing the anisotropy parameter κ leads to a decrease in the tidal deformability.This is due to the enhanced mass, as the tidal deformability is inversely proportional to the compactness.This means that EOS that were previously ruled out based on isotropic studies of the Λ 1.4 constraint are still viable candidates if we consider their anisotropic pressure modification.The sensitivity of Λ on κ is extremely important as it may help to ultimately establish observational bounds on this parameter.
We finally examine the EOS independence of the I-Love-C universal relations.We see that of the three effects -EOS, magnetic field, and anisotropy -it is the anisotropy that has a maximal effect of varying the universal fit between the parameters.Nevertheless, the deviation from universality is always with RMSE of the or-der of 1%, which is still much smaller than experimental errors on, say, the moment of inertia I.
In summary, NSs could possibly be mass-gap candidates.However, the EOS alone is unlikely to yield a massive NS with a mass exceeding 2.5M ⊙ .The introduction of anisotropy, arising from either matter, magnetic fields, or a combination of both, appears to be a crucial factor in achieving this.However, the magnetic field geometry must be carefully controlled to prevent instability and comply with the constraint imposed by tidal deformability.
With the appropriate magnetic field geometry and strength, it becomes feasible to have a massive, stable NS that fits within the mass-gap while satisfying the tidal deformability constraint.It is important to note that we have not considered the effects of rotation.Building on the findings of previous work [17], an additional increase in mass will occur, which should facilitate the interpretation of compact stellar objects within the mass gap as NSs rather than low-mass black holes.This aspect will be a subject of investigation in future research.

FIG. 5 :
FIG. 5: M-R curves for varying B 0 and orientation in case of the profile corresponding to η = 0.2, γ = 2 for the DD-ME2 EOS.Anisotropic parameter κ is set to 0.5 throughout.

FIG. 7 :FIG. 8 :
FIG. 7: Tidal deformability Λ as a function of M for varying B 0 in case of the profile corresponding to η = 0.2, γ = 2 for the DD-ME2 EOS.Anisotropic parameter κ is set to 0.5 throughout.

FIG. 10 :
FIG. 10: M-R curves for varying B 0 in case of the field profile with η = 0.01, γ = 2 for the DDMEX EOS.Anisotropic parameter, κ, is set to 0.5 throughout.Dotted lines indicate RO fields, dashed lines indicate TO fields.

FIG. 11 :
FIG. 11: Tidal deformability Λ as a function of M for varying B 0 in case of profile corresponding to η = 0.01, γ = 2 for the DDMEX EOS.Anisotropic parameter κ is set to 0.5 throughout.The red, dotted lines show the limits on Λ 1.4 from GW190814.The various curves are the same as Fig. 10.

FIG. 14 :
FIG. 14: Change in Λ − M relation of the DDMEX EOS due to varying anisotropy parameter κ.B 0 has been fixed to 0 G throughout.

FIG. 15 : 4 0 2 0 3 0
FIG. 15: C − I relation for nonmagnetic NSs constructed under DDMEX EOS with varying κ.I is normalized in units of the cube of the NS's mass (M 3 ).

FIG. 18 :
FIG. 18: C − I relation for DDMEX EOS with fixed κ = 0.5.Two values of B 0 are chosen along with the nonmagnetized case.
) where µ B is the baryon chemical potential and µ is the dipole magnetic moment, respectively, expressed in MeV and A/m 2 , to obtain B in units of the electron's critical field B c = 4.414 × 10 13 G.We use coefficients, a, b, c, obtained from the numerical simulation of a 2.2M ⊙ star in previous work[38,39], which are a = −0.769G 2 /Am 2 , b = 1.20 × 10 3 G 2 /Am 2 MeV, c = −3.46× 10 −7 G 2 /Am 2 MeV 2 .As implied, this profile, in principle, is valid only for 2.2M ⊙ NS.Other NSs should have different a, b, c.

TABLE 12 :FIG. 26 :
FIG. 26: Variation of the magnetic field within the star for varying µ based on the polynomial profile with the DDMEX EOS.

TABLE 1 :
Properties of SNM at n 0 for the EOS used in this work.

TABLE 3 :
M max for the different EOS under varying B 0 and orientations.Anisotropy parameter κ = 0.5 throughout.The field profile is for η = 0.2, γ = 2. κ by itself leads to higher masses, which are further enhanced (reduced) by TO (RO) fields.

TABLE 5 :
Results for the DDMEX EOS under varying κ for three

TABLE 6 :
The best-fit parameters and calculated errors for the universal relations for NSs constructed by DDMEX EOS with varying degrees of anisotropy.

TABLE 7 :
RMSE for the universal relations for NSs constructed by DDMEX EOS with varying degrees of anisotropy.

TABLE 8 :
The best-fit parameters and calculated errors for the universal relations for NSs constructed by DDMEX EOS with fixed κ = 0.5 and varying B 0 .

TABLE 9 :
RMSE for the universal relations of NSs constructed by DDMEX EOS with fixed κ = 0.5 and varying B 0 .

Table 12 , we
15−16G, then the corresponding central field hardly affects the stellar structure.Therefore, this profile does not seem to be a viable profile, from either stability or observability.