# Neutron Stars in Scalar Torsion Theories with Nonminimal Coupling

Youcef Kehal<EMAIL_ADDRESS>Khireddine Nouicer  khnouicer@univ-
jijel.dz Laboratory of Theoretical Physics, Departement of Physics, University
Mohamed Seddik Ben Yahia, BP 98, Ouled Aissa, Jijel 18000, Algeria  Hamza
Boumaza<EMAIL_ADDRESS>Laboratoire de Physique des Particules et
Physique Statistique (LPPPS),
Ecole Normale Supérieure-Kouba, B.P. 92, Vieux Kouba, 16050 Algiers, Algeria

###### Abstract

The aim of this paper is to investigate the existence and the structure of
static and slowly rotating neutron stars (NSs) in scalar tensor theories
within the framework of teleparallel gravity, namely scalar torsion theories.
We focus on a particular type of theory $F(T,\phi,X)$ which includes a
canonical scalar field $\phi$ that is nonminimally coupled to the torsion
scalar $T$ through a coupling parameter $\xi$ plus a potential of the form
$V(\phi)=-\mu^{2}\phi^{2}/2+\lambda\phi^{4}/4$. We derive the hydrostatic
equilibrium equations in static case and solve them numerically for both
interior and exterior regions using the appropriate boundary conditions near
the center and far from the star. We plot the radial profiles of the metric
functions and the scalar field as well as extracting the mass-radius diagram
of the star employing a set of different realistic equations of state (EoSs).
Our findings show a high degree of compatibility with the observational
constraints of the GW170817 event and indicate a maximum mass of
$2.37M_{\odot}$ for BSk21 with $\xi=0.25$ which is in close proximity to the
recent discovery of the most massive NS ever observed, PSR J0952-0607. We end
by extending our analysis to include slow rotation regime and determine the
relation between the star’s moment of inertia and its mass, then we
demonstrate the insensitivity of the dimensionless moment of inertia and
normalized moment of inertia to the EoSs.

##  I Introduction

Beside its mathematical elegance when describing geometry of spacetime,
general relativity (GR) has also been successfully tested in weak
gravitational regime, it fits properly with observational data inside the
solar and stellar system scales, as it was able to accurately predict a large
number of astrophysical phenomena will2014confrontation . Nonetheless, the
theory seems to be struggling to cope with cosmological facts since that
Supernovae Ia data (SN Ia) riess1998observational ; perlmutter1999measurements
and cosmic microwave background (CMB) measurements spergel2003first  suggest
that the universe is currently passing through an accelerated expansion phase
which can not be occurred without introducing an exotic form of energy with
negative pressure, called dark energy (DE) or by replacing GR with a more
complete theory of gravity. Also, on larger scales the existence of
hypothetical dark matter (DM) is necessary to fully explain the rotation
curves of galaxies and structure formation. All these considerations, along
with the challenge of unifying Einstein’s theory with quantum mechanics, have
spurred researchers to look for alternative theories of gravity.

Efforts in this direction have gained significant attention notably after the
first detection of gravitational waves (GWs) GW150914, made by the LIGO-VIRGO
collaboration. The observation was the result of the merger of two binary
black holes with $29M_{\odot}$ and $36M_{\odot}$ into a black hole of mass
$62M_{\odot}$ abbott2016gw150914  which radiated a $3.0_{-0.5}^{+0.5}$ solar
masses of energy abbott2016observation  or soon later from GW signals GW170817
emitted by the coalescing of nearly equal mass binary neutron stars inspiral
(NSs) with combined mass sways around the order of
$2.74_{-0.04}^{+0.04}M_{\odot}$ abbott2017gw170817  followed by its
electromagnetic counterpart GW170817A detected by the Fermi GRB monitor
LIGOScientific:2017zic  spawning the establishment of multi-messenger GW
astronomy metzger2017welcome . Events of this kind can provide additional
insights into their sources and allow for exploration of physics in regions of
high density. For instance, compact objects like NS, which are originally
formed from the collapsed stellar relics, where the gravitational field is
extremely strong and the inner matter is compressed to densities of almost
$2.8\times 10^{14}g/cm^{3}$ with a mass around $M\sim 1.4-2M_{\odot}$ and
typical radius lies between $10$ km and $13$ km olmo2020stellar  are
considered as an ideal astrophysical laboratories to test modified theories of
gravity in strong field regimes.

In this context, one straightforward way to modify gravity is to substitute
the Ricci scalar $R$ in Einstein-Hilbert action with a general function of it
$f(R)$ sotiriou2010f . These theories have been extensively studied for the
purpose of explaining the accelerated expansion of the universe and the
inflationary epoch de2010f , as well as for examining the structure of NS in
spherical symmetric and static background for R-squared gravity and
Starobinsky models using perturbative cooney2010neutron ;
arapouglu2011constraints ; orellana2013structure ; astashenok2013further ;
resco2016neutron  or exact and self-consistent ganguly2014neutron ;
yazadjiev2014non ; capozziello2016mass ; feng2018equation ;
astashenok2017realistic  analyses. Furthermore, the inclusion of an extra
scalar degree of freedom (DoF) to the usual two tensor modes of GR releases
wide options to modify gravity through the so-called scalar tensor theories
(STT) quiros2019selected . The main reason for adding scalar fields to the
Lagrangian lies in evading Lovelock theorem assumptions whereupon screening
mechanisms must be intercede to circumvent the local experiments. As a result,
various classes of STT have been studied in-depth to explore their
phenomenological implications on extreme relativistic objects
babichev2010relativistic ; sotani2012slowly ; minamitsuji2016relativistic . In
particular, the existence and the proprieties of NS solutions were treated in
the context of Brans-Dicke theories kase2019neutron  using massive and self-
interacting potentials for different coupling values between the scalar field
and the gravitational sector while the works of arapouglu2019neutron ;
arapouglu2020neutron ; odintsov2022neutron  are devoted to study the effects
of a non-minimally coupled scalar field in Jordan frame on the mass-radius
profile of the star employing multiple potential choices.

Additionally, a larger class of STT known as Horndeski theories
kobayashi2019horndeski  along with their generalization customary dubbed as
DHOST whose Lagrangian is characterized by a second-order derivatives of the
scalar field and immunity against Ostrogradski instabilities
langlois2016degenerate  have analyzed the NS configuration cisterna2016slowly
; boumaza2022neutron  in the Hartle-Thorne slow-rotation approximation
hartle1967slowly . This procedure, when combined with the consideration of
massive scalar fields offers more substantial deviation from GR since the
finite range of the scalar field determined by its Compton wave length and the
observations outside the star beyond this scale can not set tight constraints
on the parameters of massive STT as opposed to the massless case
yazadjiev2016slowly ; staykov2018static . In addition to their standard
properties, NSs have been found to possess new features in certain theories.
These features include the acquisition of non-trivial scalar hairs, a
discovery first made by the Damour and Esposito-Farese (DEF) model
PhysRevLett.70.2220  and thoroughly discussed in salgado1998spontaneous ;
pani2014slowly ; yagi2021neutron ; hu2021scalarized ; ventagli2021neutron ,
the manifestation of spontaneous vectorization in Proca theories
kase2020neutron  and the possible existence of new types of NSs with
topological charges doneva2020topological . And we also note that a lot of
works sought after the exploration of the universal relations between the
normalized moments of inertia $I$ and the stellar compactness of the star $C$
breu2016maximum ; staykov2016moment ; popchev2019moment . In order to broaden
the picture, and stimulated with some astrophysical standpoints such as the
recycled pulsars and the highly magnetized stars where the rotational effect
becomes considerable, the study of NS dynamics could be extended even further
to incorporate the rapid rotating regime doneva2013rapidly ;
yazadjiev2015rapidly .

Along the same line, to describe gravity in cosmological and astrophysical
systems one can shift the perspective by choosing to work in context of non-
Riemannian teleparallel geometry bahamonde2021teleparallel . The teleparallel
equivalent of GR (TEGR) can be identified as a gauge theory for the
translation group in which the gravitational interaction is not geometrized
but it rather incarnates itself in torsion $T$ aldrovandi2012teleparallel ;
krvsvsak2019teleparallel . By mimicking $f(R)$ models, the TEGR Lagrangian can
be promoted to $f(T)$ cai2016f  or enrich it with scalar fields $f(T,\phi)$
geng2011teleparallel ; horvat2015nonminimally  and more general constructions
$f(T,\phi,X,Y)$ that comprises two other scalar quantities: the kinetic energy
term $X$ and the derivative coupling term $Y$, namely scalar torsion theories
hohmann2018scalar1 ; hohmann2018scalar2 ; hohmann2018scalar3 . Although, the
modified teleparallel theories are commonly known in literature from being not
invariant under local Lorentz transformations li2011f , this issue was cured
by the restoration of the spin connection as a dynamical variable alongside
the tetrads, this is the covariant version of TEGR krvsvsak2016covariant .
Therefore, it was alluring to display how much $f(T)$ theories can influence
the configuration of stars in static and spherically symmetric spacetime
boehmer2011existence ; lin2022realistic ; fortes2022solving  or adopt the
scalar torsion theories to describe slowly rotating NS boumaza2021slowly .

The primary focus of this manuscript is to investigate the existence of NS
solutions in a particular subfamily of teleparallel scalar tensor theories,
denoted as $F(T,\phi,X)$ whose form incorporates a nonminimal coupling
function between the scalar field $\phi$ and the torsion scalar $T$ with a
coupling parameter $\xi$. We analyze the ramifications of including a scalar
potential of the form $V(\phi)=-\mu^{2}\phi^{2}/2+\lambda\phi^{4}/4$ on the
structure of a spherically symmetric and static NS, this potential was mainly
examined alongside a coupling function in the context of the symmetron
screening mechanism hinterbichler2010screening . By deriving the modified
Tolman-Oppenheimer-Volkoff (TOV) equations we show that the exterior solution
is a Schwarzschild-de Sitter type where the effective cosmological constant is
induced by considering a non-vanishing asymptotic potential
$V(\phi_{\infty})$. Our goal is to extract the mass-radius relations and
compare our results with the available observations data of GWs using four
different realistic equations of state (EoSs): SLy, BSk19, BSk20 and BSk21.
Specifically, The GW170817 event has reported that a static NS with a mass of
approximately $1.6M_{\odot}$ must have a radius greater than or equal to
$10.68$ km while the maximum mass must have a radius larger than $9.6$ km.
These findings are critical in constraining the EoS for the matter at
extremely high densities, as well as the parameters of any given theory. Our
numerical simulations have been carried out to investigate the compatibility
of these constraints with our model and have yielded a maximum mass of
$2.37M_{\odot}$ for BSk21 with $\xi=0.25$, this places us in the immediate
vicinity of the groundbreaking discovery of PSR J0952-0607 romani2022psr
which represents the most massive and fastest-spinning NS discovered so far.
Throughout our investigation, we aim to shed light on the inertia-mass
relations by extending the study into slow rotation approximation. For this
purpose, we adopt a Weitzenböck regular tetrad to generate an axial symmetric
spacetime then we finish the paper by demonstrating the EoS-independence of
the dimensionless and the normalized moments of inertia.

The layout of this paper is organized as follows. In section (II) we set up
the stage of our work by presenting the general formalism and the field
equations of the scalar torsion theories $F(T,\phi,X)$ then we briefly discuss
the notions of axial symmetry in teleparallel gravity. In section (III), we
derive the general modified TOV equations for static case as well as the
initial conditions and the asymptotic behavior of the metric functions and the
scalar field. The (IV) is devoted to display our numerical results and the
mass-radius graphs of the star. Finally in (V) we explore the Hartle-Thorne
formalism of slowly rotating NS and illustrate the relation between the moment
of inertia and the mass of the star.

##  II Scalar Torsion Theories

###  II.1 The formalism

This subsection aims to present briefly the main ingredients of teleparallel
theory of gravity in which the dynamical variables are usually given by a
tetrad fields $e^{a}=e^{a}{}_{\mu}dx^{\mu}$ and a flat spin connection
$\omega^{a}{}_{b}=\omega^{a}{}_{b\mu}dx^{\mu}$. The tetrads are a set of four
linearly independent vectors that form an orthonormal basis in the tangent
space at each point of the spacetime manifold $(\mathcal{M},g_{\mu\nu})$ while
the spin connections are one-forms that represent an additional degree of
freedom (DoF) related to inertial effects. Note that throughout this article,
the Greek letters $(\mu,\nu,\rho...=0,1,2,3)$ denote the external spacetime
coordinates whereas the Latin letters $(a,b,c...=0,1,2,3)$ denote the internal
Lorentz coordinates. The metric tensor $g_{\mu\nu}$ can be constructed via the
relation

$g_{\mu\nu}=\eta_{ab}e^{a}{}_{\mu}e^{b}{}_{\nu}\,,$ (1)

where $\eta_{ab}=\mathrm{diag}(-1,1,1,1)$ is the Minkowski metric associated
with the tangent space. The tetrad components are normalized by the
orthogonality conditions $e^{a}{}_{\nu}e_{a}{}^{\mu}=\delta^{\mu}_{\nu}$ and
$e^{a}{}_{\mu}e_{b}{}^{\mu}=\delta^{a}_{b}$. Under a local Lorentz
transformations (LLTs) $\Lambda^{a}{}_{b}(x)$ the metric tensor (1) remains
invariant and for a fully covariant approach both the tetrads and spin
connections must transform simultaneously bahamonde2021teleparallel

$e^{\prime a}{}_{\mu}=\Lambda^{a}{}_{b}\,e^{b}{}_{\mu}\,,\quad\omega^{\prime
a}{}_{b\mu}=\Lambda^{a}{}_{c}\,\omega^{c}{}_{d\mu}(\Lambda^{-1})^{d}{}_{b}+\Lambda^{a}{}_{c}\,\partial_{\mu}(\Lambda^{-1})^{c}{}_{b}\,.$
(2)

The role of the spin connection is to compensate inertial effects
krvsvsak2019teleparallel  and its existence is necessary to avoid the breaking
of LLTs and make the whole theory frame-independent krvsvsak2016covariant .
Locally it is convenient to work within a purely inertial gauge under which
the spin connection vanishes identically, this is often stated as the
Weitzenböck gauge

$\omega^{a}{}_{b\mu}=\Lambda^{a}{}_{c}\,\partial_{\mu}\Lambda_{b}{}^{c}\,.$
(3)

In contrast to the torsion-less Levi-Civita connection
${\overset{\circ}{\Gamma}}\vphantom{\Gamma}^{\rho}{}_{\mu\nu}$ used in GR
where the gravitational interaction is geometrized by curvature
aldrovandi2012teleparallel , the formulation of teleparallel geometry adopt
the so-called Weitzenböck connection $\Gamma^{\rho}{}_{\mu\nu}$ which is the
unique linear affine connection that satisfies metric compatibility condition,
possesses torsion and yields to a vanishing curvature thus the field strength
of gravity is represented by the torsion tensor

$T^{\rho}{}_{\mu\nu}=-2\Gamma^{\rho}{}_{[\mu\nu]}=e_{a}{}^{\rho}\left(\partial_{\mu}e^{a}{}_{\nu}-\partial_{\nu}e^{a}{}_{\mu}+\omega^{a}{}_{b\mu}e^{b}{}_{\nu}-\omega^{a}{}_{b\nu}e^{b}{}_{\mu}\right)\,,$
(4)

noteworthy to mention that the torsion tensor is invariant under LLTs (2) and
diffeomorphisms. The tensor that relates the Levi-Civita and Weitzenböck
connections namely the contortion tensor is given by Ricci theorem

$\begin{split}K^{\rho}{}_{\mu\nu}&=\Gamma^{\rho}{}_{\mu\nu}-\overset{\circ}{\Gamma}\vphantom{\Gamma}^{\rho}{}_{\mu\nu}\\\
&=\frac{1}{2}\left(T_{\mu\ \nu}^{\ \rho}+T_{\nu\ \mu}^{\ \rho}-T_{\
\mu\nu}^{\rho}\right)\,.\end{split}$ (5)

It is also useful to define the superpotential tensor as follows

$S_{\rho}{}^{\mu\nu}=K^{\mu\nu}{}_{\rho}-\delta_{\rho}^{\mu}T_{\sigma}{}^{\sigma\nu}+\delta_{\rho}^{\nu}T_{\sigma}{}^{\sigma\mu}=-S_{\rho}{}^{\nu\mu}\,,$
(6)

therefore, the torsion scalar can be written in compact way

$T=\frac{1}{2}S_{\rho}{}^{\mu\nu}T^{\rho}{}_{\mu\nu}\,.$ (7)

In fact, the Lagrangian of TEGR is established from a linear contribution of
the torsion scalar $T$, this latter differs from the Ricci scalar
$\overset{\circ}{R}$ only by a total divergence term $B$ such that

$\begin{split}\overset{\circ}{R}&=-T+B\\\
&=-T+\frac{2}{e}\partial_{\mu}\left(eT_{\lambda}{}^{\lambda\mu}\right)\,,\end{split}$
(8)

where $e=\sqrt{-g}$ is the determinant of the tetrads $e^{a}{}_{\mu}$. The
above relation (8) shows that the Lagrangians of TEGR and GR are equivalent to
each other up to a surface term which ensures that their corresponding field
equations and thus all the phenomenological predictions at astrophysical and
cosmological levels are the same. This equivalence is represented by the
following identity

$G_{\alpha\beta}=\overset{\circ}{\nabla}\vphantom{\nabla}_{\rho}S_{\beta\alpha}{}^{\rho}-\left(K^{\rho}{}_{\beta\nu}+T^{\rho}{}_{\beta\nu}\right)S_{\rho\alpha}{}^{\nu}+\frac{1}{2}g_{\alpha\beta}T\,,$
(9)

where $G_{\alpha\beta}$ is the Einstein tensor and the Ricci theorem (5) can
be expressed in terms of spin connections as
$K^{b}{}_{a\mu}=\omega^{b}{}_{a\mu}-\overset{\circ}{\omega}\vphantom{\omega}^{b}{}_{a\mu}$
while the Levi-Civita covariant derivative is defined as

$\overset{\circ}{\nabla}\vphantom{\nabla}_{\nu}S_{a}{}^{\mu\nu}=e^{-1}\partial_{\nu}\left(eS_{a}{}^{\mu\nu}\right)-{\overset{\circ}{\omega}\vphantom{\omega}^{b}{}_{a\nu}}S_{b}{}^{\mu\nu}\,.$
(10)

Now, we come to the bottom line of this subsection. After defining the torsion
scalar (7) one can also generalize the TEGR Lagrangian by adding scalar fields
$\phi$ and its kinetic energy term
$X=-\left(1/2\right)g^{\mu\nu}\phi_{,\mu}\phi_{,\nu}$. So the total action
reads as

$\mathcal{S}=\int
d^{4}xe\left[-\frac{1}{2}F\left(T,\phi,X\right)+\mathcal{L}_{m}\right]\,,$
(11)

where $F$ is a function of all its arguments and $\mathcal{L}_{m}$ refers to
the matter Lagrangian density, we suppose here that the matter fields is
minimally coupled to the scalar field but not allowed to couple directly with
the spin connection hohmann2018scalar1 . The energy-momentum tensor is defined
as

$\Theta_{\mu\nu}=-\frac{2}{e}\frac{\delta\left(e\mathcal{L}_{m}\right)}{\delta
g^{\mu\nu}}=e^{a}{}_{\mu}\left(\frac{1}{e}\frac{\delta\left(e\mathcal{L}_{m}\right)}{\delta
e^{a}{}_{\lambda}}\right)g_{\nu\alpha}=e^{a}{}_{\mu}\Theta_{a}{}^{\lambda}g_{\nu\alpha}\,.$
(12)

By ignoring boundary terms, the variations of the action (11) with respect to
the tetrad fields yields to

$-\frac{1}{2}Fe_{a}{}^{\mu}+\overset{\circ}{\nabla}\left({F_{T}S_{a}{}^{\mu\nu}}\right)-F_{T}\left(T^{\sigma}{}_{a\nu}+K^{\sigma}{}_{a\nu}\right)S_{\sigma}{}^{\mu\nu}+\frac{1}{2}F_{X}g^{\mu\nu}e_{a}{}^{\rho}\phi_{,\nu}\phi_{,\rho}=\Theta_{a}{}^{\mu}\,.$
(13)

Contracting with $e^{a}{}_{\beta}g_{\mu\alpha}$ and using the identity (9) we
can write the tetrad field equations in covariant form

${G}_{\alpha\beta}F_{T}+\frac{1}{2}g_{\alpha\beta}\left(F-F_{T}T\right)+S_{\beta\alpha}{}^{\rho}{}\overset{\circ}{\nabla}_{\rho}F_{T}+\frac{1}{2}F_{X}\phi_{,\alpha}\phi_{,\beta}=\Theta_{\alpha\beta}\,.$
(14)

Taking variation of (11) with respect to the scalar field $\phi$ leads to the
Klein-Gordon equation

$g^{\mu\nu}\overset{\circ}{\nabla}_{\mu}\left(F_{X}\phi_{,\nu}\right)+F_{\phi}=0\,,$
(15)

where $F_{A}\equiv\partial F/\partial A$ such that $A=T,X,\phi$. Finally, the
antisymmetric part of (14) coincides with the variations of the action (11)
with respect to the spin connection hohmann2018scalar2

$\left(S_{\alpha\beta}{}^{\rho}-S_{\beta\alpha}{}^{\rho}\right)\overset{\circ}{\nabla}_{\rho}F_{T}=0\,.$
(16)

###  II.2 Axial Symmetric Spacetime

After introducing the fundamental concepts of teleparallel gravity and present
the field equations of the scalar-torsion theory $F(T,\phi,X)$, our goal here
is to provide a quick glimpse into the notions of symmetries within the
teleparallel geometric framework hohmann2019modified  and reach the
Weitzenböck tetrad that reproduce an axial symmetric spacetime
bahamonde2021exploring  i.e a non-diagonal metric tensor with one cross term.
The teleparallel geometry denoted as $(\mathcal{M},\bm{e},\bm{\omega})$ is
invariant under the action $\Phi$ of a Lie group $G$ on a spacetime manifold
$\mathcal{M}$ if and only if the metric and the affine connection are
invariant under the diffeomorphisms
$\Phi_{u}:\mathcal{M}\rightarrow\mathcal{M}$ such that $u\in G$. For
convenience, when discussing teleparallel symmetries we focus solely on the
infinitesimal symmetries where the diffeomorphisms $\Phi_{u}$ are generated by
a set of vector fields $X_{\varsigma}$ that preserves the geometric structure
of the spacetime $\mathcal{M}$ which in turns require that the Lie derivative
of the metric tensor as well as the teleparallel affine connection vanishes in
the direction of $X_{\varsigma}$

$\mathcal{L}_{X_{\varsigma}}g_{\mu\nu}=0,\quad\mathcal{L}_{X_{\varsigma}}\Gamma^{\lambda}{}_{\mu\nu}=0\,.$
(17)

A teleparallel geometry is symmetric only if there exists a local Lie algebra
homomorphism $\bm{\lambda}$ mapping the symmetry group into the Lorentz group
$SO(1,3)$ defined by hohmann2019modified

$\lambda_{\varsigma}(x)=\frac{d}{dt}\Lambda_{\varsigma(t)}(x)\,,$ (18)

where $\Lambda:G\times\mathcal{M}\rightarrow SO(1,3)$ is the corresponding
local Lie group homomorphism. The requirement of invariance i.e the vanishing
Lie derivatives (17) implies that the Lie derivative of the tetrads and the
spin connections must vanish

$\mathcal{L}_{X_{\varsigma}}e^{a}{}_{\mu}=-\lambda^{a}{}_{b}e^{b}{}_{\mu}\,,\quad\mathcal{L}_{X_{\varsigma}}\omega^{a}{}_{b\mu}=\mathcal{D}_{\mu}\lambda^{a}{}_{b}=\partial_{\mu}\lambda^{a}{}_{b}+\omega^{a}{}_{c\mu}\lambda^{c}{}_{b}-\omega^{c}{}_{b\mu}\lambda^{a}{}_{c}\,,$
(19)

these are the teleparallel Killing equations. Clearly, it is much easier to
solve them in Weitzenböck gauge by picking a particular choice of
$\Lambda^{a}{}_{b}$ where the spin connections vanish. Doing so, the second
condition of (19) implies that the local Lie algebra homomorphism $\lambda$
are spacetime independent and hence the symmetric tetrads satisfying the first
condition of (19) are related by a global Lorentz transformation
hohmann2019modified . It can be stated that the symmetries of the tetrad-spin
connection pair guarantee that the torsion tensor also satisfies these
symmetries. and we require that the symmetry conditions are also fulfilled by
the scalar fields bahamonde2021exploring

$\mathcal{L}_{X_{\varsigma}}T^{\rho}{}_{\mu\nu}=0\,,\quad\mathcal{L}_{X_{\varsigma}}\phi=0\,.$
(20)

In spherical coordinates $(t,r,\theta,\varphi)$ the single Killing vector
field that forms the $so(2)$ algebra and generates the $SO(2)$ group is
$X_{z}=\partial_{\varphi}$, this lies in generating a rotation $R(\varphi)\in
SO(2)$ around the polar axis. Within the Weitzenböck gauge, one needs to
determine the Lie algebra homomorphism $\lambda:g\rightarrow so(1,3)$
providing that $\lambda$ is constant over the spacetime. Obviously, multiple
choices can be made but we restrict ourselves to a nontrivial canonical one
that yields to the regular branch of the Weitzenböck tetrads whose components
take the form hohmann2019modified

$e^{0}{}_{\mu}=C^{0}{}_{\mu}\,,\quad
e^{1}{}_{\mu}=C^{1}{}_{\mu}\cos\varphi-C^{2}{}_{\mu}\sin\varphi\,,\quad
e^{2}{}_{\mu}=C^{1}{}_{\mu}\sin\varphi+C^{2}{}_{\mu}\cos\varphi\,,\quad
e^{3}{}_{\mu}=C^{3}{}_{\mu}\,,$ (21)

this branch possesses a proper limit to spherically symmetric tetrads where
the functions $C^{a}{}_{\mu}$ depend on the coordinates $(r,\theta)$ after
assuming stationarity. In general, substituting the tetrad (21) in (1) will
produce a non-circular metric with more than just one off-diagonal component
so to obtain an axial symmetric spacetime with only one cross term, we impose
a particular coordinate gauge to get rid of all the off-diagonal components of
the metric tensor except for the $dtd\varphi$ component bahamonde2021exploring

$\displaystyle
C^{0}{}_{1}=C^{1}{}_{0}=C^{0}{}_{2}=C^{2}{}_{0}=C^{1}{}_{3}=C^{2}{}_{1}=C^{2}{}_{2}=C^{3}{}_{0}=C^{3}{}_{3}=0\,$
$\displaystyle\qquad\qquad\qquad C^{2}{}_{3}=-C^{3}{}_{2}\,,\quad
C^{3}{}_{1}=-\frac{C^{1}{}_{1}C^{1}{}_{2}}{C^{2}{}_{3}}\,,$ (22)

and we are only left with $05$ functions to determine. The components of the
axial symmetric metric for the reduced tetrad are given by
bahamonde2021exploring

$\displaystyle
g_{tt}=-\left(C^{0}{}_{0}\right)^{2},g_{rr}=\left(C^{1}{}_{1}\right)^{2}\left[\left(\frac{C^{1}{}_{2}}{C^{2}{}_{3}}\right)^{2}+1\right],g_{\theta\theta}=\left(C^{1}{}_{2}\right)^{2}+\left(C^{2}{}_{3}\right)^{2}$
$\displaystyle\quad\quad\quad\quad\quad
g_{\varphi\varphi}=\left(C^{2}{}_{3}\right)^{2}-\left(C^{0}{}_{3}\right)^{2},g_{t\varphi}=-C^{0}{}_{0}C^{0}{}_{3}\,.$
(23)

We pointed out that the presence of the free function $C^{0}{}_{3}$ induces a
rotation around the polar axis while in the case of static and spherically
symmetric spacetime $C^{0}{}_{3}$ is expected to vanish.

##  III Static and Spherically Symmetric Background

###  III.1 Modified Tolman-Oppenheimer-Volkoff equations

We are looking for solutions of the field equations assuming a static and
spherically symmetric background metric with two independent functions of
radial coordinate $f$ and $g$, that is

$ds^{2}=-f(r)dt^{2}+h(r)dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2})\,,$ (24)

for the matter sector which is encapsulated in the energy-momentum tensor, we
shall also model the stellar matter inside the neutron star as a perfect fluid

$\Theta_{\mu\nu}=(\rho+P)u_{\mu}u_{\nu}+Pg_{\mu\nu}\,,$ (25)

where $\rho=\rho(r)$ and $P=P(r)$ represent the matter density and the
pressure respectively, while $u^{\mu}$ is the fluid’s four velocity. Recall
that the fluid is at rest so the normalization condition $u^{\mu}u_{\mu}=-1$
implies that $u^{\mu}=(f(r)^{-1/2},0,0,0)$. From the contracted Bianchi
identity $\nabla^{\mu}\Theta_{\mu\nu}=0$, we obtain the continuity equation

$P^{\prime}+\frac{f^{\prime}}{2f}(\rho+P)=0\,.$ (26)

Hereafter, the prime denotes a derivative with respect to $r$. Among an
infinity possible tetrads that can recover the metric (24), the simplest
tetrad is the diagonal gauge one

$\tilde{e}^{a}{}_{\mu}=\mbox{diag}(\sqrt{f(r)},\sqrt{h(r)},r,r\mbox{sin}\theta)\,.$
(27)

The literature commonly identifies this tetrad as a "bad" choice for
spherically symmetric solutions in $f(T)$ gravity tamanini2012good  as it
gives rise to extra off-diagonal component of the field equations which in
turn restricts the functional form of $f(T)$ and reduces the theory to TEGR
case boehmer2011existence . However, it was shown in krvsvsak2016covariant
that this issue is an outcome of the old pure-tetrad formulation of
teleparallel gravity which assumes a vanishing spin connection, since that the
latter enters the TEGR action as a total divergence term so any arbitrary
choice of the spin connection will not affect the TEGR field equations and
leave them invariant under LLTs, but the ignorance of the role of spin
connections causes a delicate problem when we promote the action and consider
a non-linear general function of $T$. The $f(T)$ gravity in standard
formulation suffers clearly from Lorentz invariance violation. The solution of
this problem lies in following a covariant version of $f(T)$
krvsvsak2016covariant  and calculate the associated spin connections for each
tetrads using the switch-off gravity technique presented in krvsvsak2015spin .
The non-vanishing components of the spin connection for the diagonal tetrad
(27) are given by

$\omega_{\ \hat{2}\theta}^{\hat{1}}=-1\,,\quad\omega_{\
\hat{3}\phi}^{\hat{1}}=-\sin\theta\,,\quad\omega_{\
\hat{3}\phi}^{\hat{2}}=-\cos\theta\,.$ (28)

Through an ansatz choice based on symmetry principles hohmann2019modified  and
discussed above in subsection (II.2), one can also reproduce the metric (24)
by using a rotated tetrad $e^{a}{}_{\mu}$ related to the diagonal one (27) by
a Lorentz transformation matrix (2) such that

$C^{0}{}_{0}=\sqrt{f(r)}\,,\quad C^{1}{}_{1}=\sqrt{h(r)}\sin\theta\,,\quad
C^{1}{}_{2}=r\cos\theta\,,\quad C^{2}{}_{3}=r\sin\theta\,,\quad
C^{0}{}_{3}=0\,,$ (29)

inserting the functions (29) in (21) taking into account the constraint
relations (II.2), the non diagonal tetrad takes the form

$e^{a}{}_{\mu}=\left(\begin{array}[]{cccc}\sqrt{f(r)}&0&0&0\\\
0&\sqrt{h(r)}\sin(\theta)\cos(\varphi)&r\cos(\theta)\cos(\varphi)&-r\sin(\theta)\sin(\varphi)\\\
0&\sqrt{h(r)}\sin(\theta)\sin(\varphi)&r\cos(\theta)\sin(\varphi)&r\sin(\theta)\cos(\varphi)\\\
0&\sqrt{h(r)}\cos(\theta)&-r\sin(\theta)&0\end{array}\right)\,,$ (30)

all components of the spin connection associated with this tetrad are null and
the antisymmetric field equation (16) is satisfied. So substituting the tetrad
(30) and the vanishing spin connection or the diagonal tetrad (27) with their
associated spin connections (28) to calculate the torsion scalar (7) or the
boundary term in (8) will lead to the same result

$\displaystyle T$ $\displaystyle=$
$\displaystyle\frac{2\left(\sqrt{h}-1\right)\left[rf^{\prime}+f\left(1-\sqrt{h}\right)\right]}{r^{2}fh}\,,$
(31) $\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{4f^{2}\left(2h^{3/2}+rh^{\prime}-2h\right)+r^{2}f^{\prime
2}h+rf\left[(4f^{\prime}h^{3/2}+rf^{\prime}h^{\prime}-2h\left(rf^{\prime\prime}+4f^{\prime}\right)\right]}{2r^{2}f^{2}h^{2}}\,.$
(32)

Using any tetrad with its corresponding spin connection and the above torsion
scalar (31) in (14) we finally obtain the $"tt"$, $"rr"$ and $\theta\theta$
components of the field equations in general form

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle-F+\left[\frac{2}{rh}\left(\left(\sqrt{h}-1\right)\frac{f^{\prime}}{f}+\frac{h^{\prime}}{h}\right)+\frac{4}{r^{2}h}\left(\sqrt{h}-1\right)\right]F_{T}+\frac{4}{rh}\left(\sqrt{h}-1\right)F^{\prime}_{T}\,,$
(33) $\displaystyle P$ $\displaystyle=$ $\displaystyle
F+\left[-\frac{2}{rh}\left(\sqrt{h}-2\right)\frac{f^{\prime}}{f}-\frac{4}{r^{2}h}\left(\sqrt{h}-1\right)\right]F_{T}+\frac{\phi^{\prime
2}}{h}F_{X}\,,$ (34) $\displaystyle P$ $\displaystyle=$ $\displaystyle
F+\left[\frac{2\left(\sqrt{h}-1\right)^{2}}{r^{2}h}-\frac{f^{\prime}}{f}\left(\frac{rhf^{\prime}+f\left(-6h+4h^{3/2}+rh^{\prime}\right)}{2rfh^{2}}\right)-\frac{h^{\prime}}{rh^{2}}+\frac{f^{\prime\prime}}{fh}\right]F_{T}$
(35)
$\displaystyle-\left[\frac{2f\left(\sqrt{h}-1\right)-rf^{\prime}}{rfh}\right]F^{\prime}_{T}\,.$

Similarly, solving the scalar field equation (15) yields to

$\phi^{\prime\prime}=\left[-\frac{2}{r}+\frac{1}{2}\left(\frac{h^{\prime}}{h}-\frac{f^{\prime}}{f}\right)-\frac{1}{2}\frac{F^{\prime}_{X}}{F_{X}}\right]\phi^{\prime}-\frac{hF_{\phi}}{F_{X}}\,.$
(36)

Throughout this paper, we propose to study a specific type of gravitational
Lagrangian whose form can be written as

$\displaystyle F(T,\phi,X)$ $\displaystyle=F(\phi)T-X+V(\phi)\,,$ (37)
$\displaystyle=\frac{\kappa}{2}T+\frac{\xi}{2}T\phi^{2}+\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi+V(\phi)\,,$

this action is written in Jordan frame with the nonminimal coupling function
$F(\phi)=1+(\xi/\kappa)\phi^{2}$ such that $\xi$ is the coupling parameter and
$\kappa=1/8\pi G$ . The second term represents a kinetic energy term preceded
by a minus sign to avoid the appearance of ghosts and tachyonic instabilities
while the last one is a general potential of the scalar field. We also note
that the action (37) belongs to a more general class of theories known as
scalar torsion theories hohmann2018scalar1 ; hohmann2018scalar2 ;
hohmann2018scalar3  and the functional form is similar to the non-minimal
quintessence model in TEGR geng2011teleparallel . Clearly, the GR action can
be recovered when setting $\xi=0$. Using our model (37) the general solutions
of the field equations (33) and (34) reduce to

$\displaystyle\frac{f^{\prime}}{f}$ $\displaystyle=$
$\displaystyle\frac{h-1}{r}+\frac{2r^{2}h(P-V)+r^{2}\phi^{\prime
2}}{2r\left(\kappa+\xi\phi^{2}\right)}\,,$ (38)
$\displaystyle\frac{h^{\prime}}{h}$ $\displaystyle=$
$\displaystyle\frac{1-h}{r}+\frac{2r^{2}h(\rho+V)-8\xi\left(\sqrt{h}-1\right)r\phi\phi^{\prime}+r^{2}\phi^{\prime
2}}{2r\left(\kappa+\xi\phi^{2}\right)}\,.$ (39)

As expected, in the absence of scalar fields and a vanishing non-minimal
coupling parameter $\xi=0$ these equations simply converges toward the general
relativistic TOV equations

$\displaystyle\frac{f^{\prime}}{f}$ $\displaystyle=$
$\displaystyle\frac{h-1}{r}+\frac{rhP}{\kappa}\,,$ (40)
$\displaystyle\frac{h^{\prime}}{h}$ $\displaystyle=$
$\displaystyle\frac{1-h}{r}+\frac{rh\rho}{\kappa}\,.$ (41)

Expressing the torsion scalar (31) in terms of (38) leads to

$T=(\sqrt{h}-1)\left[\frac{2(h-\sqrt{h})}{r^{2}h}+\frac{2h\left(P-V\right)+\phi^{\prime
2}}{h\left(\kappa+\xi\phi^{2}\right)}\right]\,,$ (42)

this shows that only the pressure, the potential and the kinetic energy
contribute to $T$ while the contribution of matter density and the coupling
parameter is transparent. Now we replace the reduced equations (38) and (39)
in (36) to obtain

$\displaystyle\phi^{\prime\prime}$ $\displaystyle=$
$\displaystyle\left[\frac{r^{2}h\left(\rho-P+2V\right)}{2r\left(\kappa+\xi\phi^{2}\right)}-\frac{1+h}{r}\right]\phi^{\prime}+\xi
hT\phi+hV_{\phi}\,,$ (43)

where $V_{\phi}\equiv\partial V/\partial\phi,$ and the continuity equation
(26) becomes

$P^{\prime}=-\left(\rho+P\right)\left(\frac{h-1}{2r}+\frac{2rh\left(P-V\right)+r\phi^{\prime
2}}{4\left(\kappa+\xi\phi^{2}\right)}\right)\,.$ (44)

Together (38),(39),(43) and (44) describe the interior structure of the star
and represent the modified TOV equations of hydrostatic equilibrium in static
and spherically symmetric background for a canonical scalar field $\phi$ non-
minimally coupled to torsion scalar $T$ through the coupling parameter $\xi$
with a general potential $V(\phi)$. Before solving the modified TOV equations
numerically and plot the mass-radius diagram we will need to extract the
boundary conditions inside and outside the stars.

###  III.2 Boundary Conditions

Let us now discuss the initial conditions of the star near the center $r=0$
and its behavior at spatial infinity for our model $F(T,\phi,X)$ (37). Around
the center of NSs we need to express the metric functions, the scalar field
and pressure in Taylor series form

$f(r)=f_{c}+\sum_{n=2}^{\infty}f_{n}r^{n}\,,\quad
h(r)=h_{c}+\sum_{n=2}^{\infty}h_{n}r^{n}\,,\quad\phi(r)=\phi_{c}+\sum_{n=2}^{\infty}\phi_{n}r^{n}\,,\quad
P(r)=P_{c}+\sum_{n=2}^{\infty}P_{n}r^{n}\,,$ (45)

where $f_{c},f_{n},h_{c},h_{n},\phi_{c},\phi_{n},P_{c},P_{n}$ are constants
and the lowercase index $"c"$ stands for center. To avoid singularities and
ensure the geometrical regularity of the star at the center $r=0$ we must
impose the following boundary conditions

$f^{\prime}(r=0)=0\,,\quad
h^{\prime}(r=0)=0\,,\quad\phi^{\prime}(r=0)=0\,,\quad P^{\prime}(r=0)=0\,.$
(46)

Likewise, we can expand the potential function $V(\phi)$ and its derivative as
well

$V(\phi)=V(\phi_{c})+\sum_{n=1}^{\infty}\frac{1}{n!}\frac{d^{n}V}{d\phi^{n}}\biggl{|}_{\phi=\phi_{c}}(\phi-\phi_{c})^{n}\,.$
(47)

By replacing the Eqs. (45) and ignoring higher orders than
$\mathcal{O}(r^{2})$ in the background equations (33)-(34) and (15)-(26) we
obtain the boundary conditions around the center

$\displaystyle f(r)$ $\displaystyle=$ $\displaystyle
f_{c}\left[1+\frac{\left(\rho_{c}+3P_{c}-2V(\phi_{c})\right)}{6\left(\kappa+\xi\phi_{c}^{2}\right)}r^{2}\right]\,,$
(48) $\displaystyle h(r)$ $\displaystyle=$ $\displaystyle
h_{c}+\frac{\left(\rho_{c}+V(\phi_{c})\right)}{3\left(\kappa+\xi\phi_{c}^{2}\right)}r^{2}\,,$
(49) $\displaystyle\phi(r)$ $\displaystyle=$
$\displaystyle\phi_{c}+\frac{1}{6}V_{\phi}(\phi_{c})r^{2}\,,$ (50)
$\displaystyle P(r)$ $\displaystyle=$ $\displaystyle
P_{c}-\frac{\left(\rho+P\right)\left(\rho_{c}+3P_{c}-2V(\phi_{c})\right)}{12\left(\kappa+\xi\phi^{2}\right)}r^{2}\,,$
(51)

we remark that the influence of $V(\phi)$ and the coupling parameter $\xi$
appears explicitly on the metric functions and pressure whereas the derivative
of $V(\phi)$ modifies the scalar field $\phi$ kase2019neutron . It is
convenient to define the Arnowitt-Deser-Misner (ADM) mass of the star through
the mass function $\mathcal{M}(r)$ at spatial infinity as $M\equiv
4\pi\kappa\left(1-h^{-1}\right)r|_{r\rightarrow\infty}$. However in the
presence of scalar fields we calculate the total mass through the following
integral formula salgado1998spontaneous

$M=\lim_{r\rightarrow\infty}m(r)=4\pi\int_{0}^{\infty}r^{2}E(r)dr\,,$ (52)

such that

$E(r)=\frac{1}{\kappa+\xi\phi^{2}}\left[\rho+V(\phi)+\frac{\phi^{\prime
2}}{2h}+\frac{4\xi}{rh}\left(1-\sqrt{h}\right)\phi\phi^{\prime}\right]\,.$
(53)

Performing a similar analysis to the one presented in equations (45) and using
the above mass function (52) we obtain

$M(r)=M_{c}+\frac{2\pi\kappa\left(\rho_{c}+V\left(\phi_{c}\right)\right)}{\kappa+\xi\phi_{c}^{2}}r^{3}\,.$
(54)

We note from (54) that the NS mass depends on $r^{3}$ which means that
$h_{c}=h(r=0)=1$ and thus one can easily conclude from (42) that the torsion
scalar vanishes at the center of the star. Also it is always possible to shift
$f_{c}$ appeared in (48) into an arbitrary constant and we typically choose
$f_{c}=1$. The central value of the scalar field $\phi_{c}$ will be determined
by a shooting method to satisfy the boundary condition of the asymptotic
scalar field $\phi(r_{\infty})=\phi_{\infty}$ kase2019neutron . The internal
solution of the NS will be integrated later starting from the center $r=0$ up
to the radius of the star $r_{s}$ which is characterized by $P(r_{s})=0$, this
will indicate that outside the star both the pressure and the density smoothly
goes to zero. For the external solution, we start from the star surface
$r_{s}$ to spatial infinity $r_{\infty}$ where we usually require that
$h(r_{\infty})=1$, this condition can be trivially understood if we used the
ADM mass formula and assumed a constant mass $M$ while the time
reparametrization of (24) allows us to normalize $f(r_{\infty})$
kase2019neutron . Outside the star, the requirement of asymptotic flatness is
a crucial consideration in several works, as seen in kase2019neutron ;
arapouglu2019neutron ; arapouglu2020neutron . It lies in imposing a vanishing
potential $V(\phi_{\infty})=0$ when the scalar field goes over to a constant
asymptotic value $\phi_{\infty}$ but this is not strictly necessary from a
physical standpoint as it is a natural consequence of postulating the
existence of isolated NSs at spatial infinity sotiriou2012black . In other
words, the imposition of a vanishing potential $V(\phi_{\infty})=0$ serves to
ensure that the gravitational field of the NS has a finite range, thus
limiting its influence on the behavior of other objects. However, throughout
our paper we will assume a priori that the NS is surrounded by a cosmological
background and we suppose an asymptotically non-vanishing potential
$V(\phi_{\infty})\neq 0$ such that the value of the scalar field at spatial
infinity rather satisfies the following condition

$\Box\phi=\frac{dV_{\text{eff}}}{d\phi}\biggl{|}_{\phi=\phi_{\infty}}=0\,,$
(55)

this equation could be easily deduced from plugging (37) in (15) where the
effective potential is given by

$V_{\text{eff}}(\phi)=\frac{1}{2}\xi T\phi^{2}+V(\phi)\,.$ (56)

Before closing this subsection, let us discuss the behavior of the metric
functions for external configurations $P=\rho=0$ with a constant asymptotical
value of the scalar field $\phi_{\infty}$. The equations (38) and (39) reduce
to

$\displaystyle\frac{f^{\prime}}{f}$ $\displaystyle=$
$\displaystyle-\frac{\kappa+\xi\phi_{\infty}^{2}-h\left(\kappa+\xi\phi_{\infty}^{2}-r^{2}V(\phi_{\infty})\right)}{r\left(\kappa+\xi\phi_{\infty}^{2}\right)}\,,$
(57) $\displaystyle\frac{h^{\prime}}{h}$ $\displaystyle=$
$\displaystyle\frac{\kappa+\xi\phi_{\infty}^{2}-h\left(\kappa+\xi\phi_{\infty}^{2}-r^{2}V(\phi_{\infty})\right)}{r\left(\kappa+\xi\phi_{\infty}^{2}\right)}\,,$
(58)

whose solutions are simply given by

$f(r)=h(r)^{-1}=1+\frac{c_{1}}{r}-\frac{\Lambda_{\text{eff}}}{3}r^{2}\,,\quad\text{with}\quad\Lambda_{\text{eff}}=\frac{V(\phi_{\infty})}{\kappa+\xi\phi_{\infty}^{2}}\,.$
(59)

These solutions represent a Schwarzschild de Sitter spacetime (SdS) where the
constant of integration $c_{1}$ is actually the mass of the star $\mathcal{M}$
measured by an observer at infinity and $\Lambda_{\text{eff}}$ is the
effective cosmological constant induced by considering a non-vanishing
asymptotic potential. Clearly, in theories with a constant value of the
asymptotic scalar field that satisfies a vanishing potential at spatial
infinity $V(\phi_{\infty})=0$ we recover the Schwarzschild limit. On the other
hand, the assumptions of the sub-horizon limit $\Lambda_{\text{eff}}\hskip
2.84526ptr^{2}\ll 1$ and the weak field approximation $\mathcal{M}/r\ll 1$
babichev2016relativistic  allow us to expand the metric functions as
$f(r)=1+\delta_{f}(r)$ and $h(r)=1+\delta_{h}(r)$ (such that $\delta_{f}(r)\ll
1$ and $\delta_{h}(r)\ll 1$). By following the same line we perturb the
asymptotically scalar field as $\phi(r)=\phi_{\infty}+\delta_{\phi}(r)$ so the
scalar field equation (15) becomes

$\delta^{\prime\prime}_{\phi}+\frac{2}{r}\delta^{\prime}_{\phi}-V^{\prime\prime}(\phi_{\infty})\delta_{\phi}-V^{\prime}(\phi_{\infty})=0\,,$
(60)

whose solution has the particularly simple form

$\delta_{\phi}(r)=k_{1}\frac{e^{-r\sqrt{V^{\prime\prime}(\phi_{\infty})}}}{r}+k_{2}\frac{e^{r\sqrt{V^{\prime\prime}(\phi_{\infty})}}}{r\sqrt{V^{\prime\prime}(\phi_{\infty})}}-\frac{V^{\prime}(\phi_{\infty})}{V^{\prime\prime}(\phi_{\infty})}\,,$
(61)

the solution involves a mass term $\sqrt{V^{\prime\prime}(\phi_{\infty})}$
where $k_{1}$ and $k_{2}$ are constants of integration but the requirement of
a finite asymptotic scalar field implies that $k_{2}$ must vanish. It is
interesting to point out after demonstrating that the torsion scalar vanishes
outside the star one can deduce from equations (55) and (56) that the first
derivative of the potential also vanishes $V^{\prime}(\phi_{\infty})=0$
leaving out only the Yukawa-type term in (61) which ultimately follows that
the perturbed asymptotic scalar field can be expressed as

$\phi(r)=\phi_{\infty}+k_{1}\frac{e^{-r\sqrt{V^{\prime\prime}(\phi_{\infty})}}}{r}\,.$
(62)

##  IV Neutron Stars Profiles and The Mass-Radius Graphs

The crucial key for solving the system of four equations (38)-(39) and
(43)-(44) numerically is to relate the pressure and the energy density of the
matter distribution through an equation of state (EoS) that fits properly with
observational data. The outer crust and the inner core of a neutron star can
be described by a unified EoS where the pressure and the energy density are
parametrized by two logarithmic quantities
$\zeta\equiv\log_{10}\left(P/dyn.cm^{-2}\right)$ and
$\xi\equiv\log_{10}\left(\rho/g.cm^{-3}\right)$. Our numerical analysis is
performed using a set of four different realistic EoSs: SLy, BSk19, BSk20 and
BSk21 and their analytical representations are given in haensel2004analytical
; potekhin2013analytical . For numerical purposes it is much easier to deal
with squared roots of the metric function $f(r)\rightarrow
f(r)^{2}\,,h(r)\rightarrow h(r)^{2}$ and use a new dimensionless parameter $s$
instead of the radial coordinate $r$ such that kase2019neutron

$s=\text{ln}(\frac{r}{r_{0}})\,,\qquad
r_{0}=\left(\frac{8\pi\kappa}{\rho_{0}}\right)^{1/2}=89.664\,\,km\,,\qquad\rho_{0}=m_{n}n_{0}=1.6749\times
10^{14}\,\,g/cm^{-3}\,,$ (63)

where $m_{n}$ is the mass of a neutron and $n_{0}$ is the typical number of
the star. Also we introduce the following dimensionless variables for the
pressure and density in terms of the new dimensionless variable $s$

$P_{m}(r)\rightarrow
c^{2}\rho_{0}P(s)\,\,,\rho_{m}(r)\rightarrow\rho_{0}\rho(s)\,.$ (64)

In this paper, we will focus mainly on exploring the effects of a particular
potential on the NS stellar configurations whose explicit form is quartic and
comprises two positive constants $\mu$ and $\lambda$

$V(\phi)=-\frac{1}{2}\mu^{2}\phi^{2}+\frac{\lambda}{4}\phi^{4}\,,$ (65)

this potential was primary investigated alongside a coupling function in the
context of symmetron screening mechanism. The effective potential in this case
becomes

$V_{\text{eff}}(\phi)=-\frac{1}{2}\mu_{\text{eff}}^{2}\hskip
1.42262pt\phi^{2}+\frac{\lambda}{4}\phi^{4}\,,$ (66)

where $\mu_{\text{eff}}^{2}=\mu^{2}-\xi T$. Outside the star it is obvious
that $T=0$ so the critical points of the effective potential are
$\phi_{\text{cr}}=0,\pm\mu/\sqrt{\lambda}$, the first one is unstable whereas
the others are stable points. In arapouglu2019neutron ; arapouglu2020neutron
the choice of taking whether stable or unstable critical points depends on the
asymptotic potential condition $V(\phi_{\text{cr}})=0$. However, as we
mentioned earlier that this condition is not necessary and throughout our work
we will consider a non-trivial stable point
$\phi_{\text{cr}}=\mu/\sqrt{\lambda}$ that obeys (55) outside the star but do
not satisfy $V(\phi_{\text{cr}})=0$. For the potential (65) we rescale the
scalar field and the parameters $\mu$ and $\lambda$ as

$\phi\rightarrow\phi\sqrt{\kappa}\,,\quad\mu\rightarrow\mu\sqrt{\frac{\rho_{0}}{\kappa}}\,,\quad\lambda\rightarrow\lambda\frac{\rho_{0}}{\kappa^{2}}\,.$
(67)

The initial conditions (48)-(51) for the potential (65) near the center $r=0$
are

$\displaystyle f(r)$ $\displaystyle=$ $\displaystyle
1+\frac{\left(2\rho_{c}+6P_{c}+2\mu^{2}\phi_{c}^{2}-\lambda\phi_{c}^{4}\right)}{12\left(\kappa+\xi\phi_{c}^{2}\right)}r^{2}\,,$
(68) $\displaystyle h(r)$ $\displaystyle=$ $\displaystyle
1+\frac{\left(4\rho_{c}-2\mu^{2}\phi_{c}^{2}+\lambda\phi_{c}^{4}\right)}{12\left(\kappa+\xi\phi_{c}^{2}\right)}r^{2}\,,$
(69) $\displaystyle\phi(r)$ $\displaystyle=$
$\displaystyle\phi_{c}-\frac{1}{6}\left(\mu^{2}\phi_{c}-\lambda\phi_{c}^{3}\right)r^{2}\,,$
(70) $\displaystyle P(r)$ $\displaystyle=$ $\displaystyle
P_{c}-\frac{\left(\rho_{c}+P_{c}\right)\left(2\rho_{c}+6P_{c}+2\mu^{2}\phi_{c}^{2}-\lambda\phi_{c}^{4}\right)}{24\left(\kappa+\xi\phi_{c}^{2}\right)}r^{2}\,.$
(71)

Similarly the equation (54) reduces to

$M(r)=M_{c}+\frac{\pi\kappa\left(4\rho_{c}-2\mu^{2}\phi_{c}^{2}+\lambda\phi_{c}^{4}\right)}{2\left(\kappa+\xi\phi_{c}^{2}\right)}r^{3}\,.$
(72)

Figure 1: The behavior of the radial profiles of the metric functions $f$
(solid lines) and $h$ (dashed lines) and the scalar field $\phi(r)$ for SLy
using four different coupling parameters. The vertical lines are locations of
the star radius $r_{s}$

With the initial conditions for the system of equations (38)-(39) and
(43)-(44) in hand, we can now proceed to solve them numerically. The explicit
form of the initial conditions as well as the modified TOV equations in terms
of the new variable $s$ can be found in the appendix (A). To integrate the
system of equations, we must consider both interior and exterior regions of
the star. For the interior region, we start from the center of the star at
approximately $s_{0}=-10$ which is small enough to consider it as the center
$r=0$ and integrate outwards to the surface of the star at $r_{s}$, where the
pressure vanishes $P(r_{s})=0$. For the exterior region, we integrate from the
star surface $r_{s}$ to infinity ($s=15$) keeping in mind that both the
pressure and density vanish outside the star. Throughout the numerical
process, we employ a double shooting method to calculate the central value of
the scalar field $\phi_{c}$ by setting the scalar field to reach its critical
value $\phi_{cr}=\mu/\sqrt{\lambda}$ at asymptotical infinity, where the
potential parameters take the values $\mu=15.04\times 10^{-8}r_{0}^{-1}$ and
$\lambda=2.512\times 10^{-14}\rho_{0}^{-1}r_{0}^{-2}$
hinterbichler2010screening ; de2021neutron . Also we match the internal and
external solutions by imposing the continuity conditions of the metric
functions $f(r_{\text{in}})=f(r_{\text{out}})$ and
$h(r_{\text{in}})=h(r_{\text{out}})$ and the continuity of the scalar field
$\psi(r_{\text{in}})=\psi(r_{\text{out}})$ such that
$\psi(s)=\phi^{\prime}(s)$. The behavior of the metric functions $f$ and $h$,
as well as the scalar field $\phi$ inside and outside the star, are shown in
Fig. (1) for a central density $\rho_{c}=6.638\rho_{0}$ using a sample EoS
(SLy) with four different values of the coupling parameter $\xi$. The
numerical results for the radial profiles of $f$ and $h$ shows that outside
the object, the spacetime becomes almost Minkowskian, verifying that our
previous assumptions of sub-horizon limit and weak field approximation are
valid and that asymptotic flatness can be achieved, even if we did not assume
$V(\phi_{\infty})=0$ a priori. Outside the star, the scalar field is given by
(62), this expression demonstrates that the exponential decay term is
controlled by a mass term and modulated by the inverse radial distance factor
$1/r$ , which causes the scalar field to approach its asymptotic value
$\phi_{\infty}=\mu/\sqrt{\lambda}$ or equivalently
$\phi/\phi_{\infty}\rightarrow 1$ as $r$ goes to infinity.

Figure 2: The behavior of the torsion and the boundary term as well as the
Ricci scalar as functions of $s$ inside and outside the star for the SLy EoS
using four different values of central density.

On the other hand, it is well known that the gravitational interaction in
teleparallel gravity incarnates itself in torsion scalar $T$ which differs
from the Ricci scalar $R$ only by a boundary term $B$. Therefore, in Fig. (2)
we illustrate the behavior of both torsion scalar and boundary term as
functions of the radial coordinate inside and outside the star. The excellent
agreement between our theoretical result (42) and the numerical simulation can
be explained by the fact that the torsion scalar is zero at the center of the
star ($r=0$), as expected from the value of $h_{c}=1$. While for the boundary
term, it starts from a positive value for lower central densities and a
negative value for higher central densities, thus the initial value of $B$
decreases as we increase the central density $\rho_{c}$. It is noteworthy to
mention that inside the matter distribution of the star, both the torsion
scalar and the boundary term curves increase. However, the boundary curve is
always higher than the torsion curve, indicating that the contribution of the
boundary term is more significant than that of torsion in the teleparallel
framework. As we move towards the star radius, both curves drop rapidly until
they reach almost zero outside the star.

Figure 3: The Mass-radius relation for the potential (65) with
$\mu=15.04\times 10^{-8}r_{0}^{-1}$ and $\lambda=2.512\times
10^{-14}\rho_{0}^{-1}r_{0}^{-2}$ for four realistic EoSs and different
coupling parameters $\xi$. For the GR solutions $\xi=0$ (black line),
$\xi=0.05$ (Blue), $\xi=0.25$ (Green), $\xi=0.75$ (Red) and $\xi=1.5$
(Magenta).

To gain a deeper understanding of how any gravitational theory differs from
GR, it is essential to examine the relation between the mass and the radius of
the NS. In our model (37) with the potential (65) we have found that the mass-
radius relation of NS deviates significantly from the GR counterpart. The
relation (52) with (53) were used to calculate the mass-radius relation for
four different realistic EoSs and multiple coupling parameters $\xi$. The
results presented in Fig. (3) shows excellent agreement with the constraints
imposed by the GW170817 event. Specifically, the radius of a static NS of mass
$1.6M_{\odot}$ and the maximum mass has to be larger than $10.68$ km and $9.6$
km, respectively. Additionally, our findings indicate that the EoSs: SLy,
BSk20 and BSk21 can support NSs with masses exceeding $2M_{\odot}$.
Remarkably, the BSk21 yielded a maximum mass of $2.37M_{\odot}$ for a coupling
parameter value of $\xi=0.25$, which is strikingly close to the measured
pulsar mass $2.35\pm 0.17M_{\odot}$ of PSR J0952-0607, the heaviest and
fastest known NS to date that has recently been detected romani2022psr . All
the data of the maximum mass as well as the and their corresponding radii are
collected in the following table.

Equation of State | SLy | BSk19 | BSk20 | BSk21
---|---|---|---|---
$\text{The coupling parameter}\,\c{c}\xi$ | $0.2452$ | $0.2444$ | $0.2208$ | $0.2505$
$M_{max}(M_{\odot})$ | $2.13975$ | $1.94308$ | $2.25290$ | $2.37809$
$\text{The corresponding radii for}\,M_{max}\,\,\text{(km)}$ | $11.0368$ | $10.0519$ | $11.1168$ | $12.2028$
$\text{The corresponding radii for}\,1.6M_{\odot}\,\text{(km)}$ | $11.7985$ | $10.7013$ | $11.8681$ | $12.5041$
Table 1: The maximum mass and their corresponding radii

Interestingly, in Fig. (4) it was observed that increasing the coupling
parameter value from $\xi=0$ which represents the GR solutions to $\xi=0.25$
led to an increase in the maximum mass of NS. However, when the coupling
parameter value was increased beyond a specific threshold $\xi=0.25$ the
maximum mass of the NS decreased. These results suggest that an optimal range
of the coupling parameter exists, where the maximum mass of the NS is highest.

Figure 4: Evolution of the maximum mass of the NS as well as its minimum
radius in terms of the coupling parameter $\xi$

##  V Slowly Rotating Neutron Stars

The next step is to consider a slowly rotating neutron stars. We extend the
metric (24) to a stationary axially symmetric spacetime hartle1967slowly

$ds^{2}=-f(r)dt^{2}+h(r)dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2})-2\omega(r,\theta)r^{2}\sin^{2}\theta
dtd\varphi+\mathcal{O}(\Omega^{2})\,,$ (73)

the regular tetrad in the Weitzenböck gauge that corresponds to the above
axial metric (73) can be determined by satisfying both conditions (II.2) and
(29) providing that
$C^{0}{}_{3}=r^{2}f(r)^{-1/2}\sin^{2}(\theta)\,\omega(r,\theta)$

$e^{a}{}_{\mu}=\left(\begin{array}[]{cccc}\sqrt{f(r)}&0&0&\frac{r^{2}\sin^{2}\left(\theta\right)}{\sqrt{f(r)}}\omega(r,\theta)\\\
0&\sqrt{h(r)}\sin(\theta)\cos(\phi)&r\cos(\theta)\cos(\phi)&-r\sin(\theta)\sin(\phi)\\\
0&\sqrt{h(r)}\sin(\theta)\sin(\phi)&r\cos(\theta)\sin(\phi)&r\sin(\theta)\cos(\phi)\\\
0&\sqrt{h(r)}\cos(\theta)&-r\sin(\theta)&0\end{array}\right)\,,$ (74)

since this tetrad by default satisfies the symmetries discussed in the
subsection II.2, it automatically solves the antisymmetric field equations
(16) and we are only left with the symmetric tetrad field equations (14). Here
$\omega(r,\theta)$ is the frame dragging function, it represents the rate of
rotation at first order of the angular velocity $\Omega$. The latter is
assumed to rotate uniformly and sufficiently small to not disturb the other
geometrical and physical properties of the star, like the gravitational field,
density and pressure. This constraint is encoded in the condition $\Omega
r_{s}\ll c$ where $r_{s}$ is the radius of the star and $c$ is the speed of
light in vacuum.

In the same vein, to model the matter field of a slowly rotating stars we
consider a perfect fluid similar to (25) but with a uniformly rotating angular
velocity $\Omega$ such that

$u^{t}=\left[f(r)-\Omega
r^{2}sin^{2}\theta(\Omega-2\omega(r,\theta))\right]^{-1/2}\,,\quad
u^{r}=u^{\theta}=0\,,\quad u^{\varphi}=\Omega u^{t}\,,$ (75)

when we ignore terms of $\mathcal{O}(\Omega^{2})$ the $(t\varphi)$ component
of the stress energy tensor of the fluid (25) becomes

$\Theta_{\varphi}{}^{t}=\frac{r^{2}sin^{2}\theta}{f}(\Omega-\omega)(\rho+P)+\mathcal{O}(\Omega^{2})\,.$
(76)

Inserting the corresponding tetrads frame (74) in the tetrad field equations
(14), the only surviving equation to first order in $\Omega$ is

$\frac{F_{T}}{r^{4}}\partial_{r}\left(\frac{r^{4}\bar{\omega}_{r}}{\sqrt{fh}}\right)+\frac{F_{T}}{r^{2}}\sqrt{\frac{h}{f}}\frac{1}{\sin^{3}\theta}\partial_{\theta}\left(\sin^{3}\theta\,\bar{\omega}_{\theta}\right)-\frac{1}{\sqrt{fh}}\left[2\sqrt{h}\left(\Omega-\bar{\omega}\right)-r\bar{\omega}_{r}\right]\frac{F^{\prime}_{T}}{r}-2\sqrt{\frac{h}{f}}\left(\rho+P\right)\bar{\omega}=0\,,$
(77)

where the function $\bar{\omega}(r,\theta)=\Omega-\omega(r,\theta)$ is the
coordinate angular velocity of the fluid seen by a freely falling observer
hartle1967slowly . Now we solve the above equation by decomposing
$\bar{\omega}$ in terms of the Legendre polynomial $P_{l}$

$\bar{\omega}(r,\theta)=\sum_{l=1}^{\infty}\bar{\omega}_{l}(r)\left(-\frac{1}{\sin\theta}\frac{dP_{l}}{d\theta}\right)\,,$
(78)

we can finally write the radial differential equation for
$\bar{\omega}_{l}(r)$ as

$\bar{\omega}^{\prime\prime}_{l}+\left[\frac{4}{r}-\frac{1}{2}\left(\frac{f^{\prime}}{f}+\frac{h^{\prime}}{h}\right)+\frac{F^{\prime}_{T}}{F_{T}}\right]\bar{\omega}^{\prime}_{l}+\left[\frac{h\left[\left(l\left(l+1\right)-2\right)F_{T}-2r^{2}\left(\rho+P\right)\right]+2rh^{1/2}F^{\prime}_{T}}{r^{2}F_{T}}\right]\bar{\omega}_{l}-\frac{2\Omega
h^{1/2}}{r}\frac{F^{\prime}_{T}}{F_{T}}=0\,,$ (79)

the behavior of this equation for our gravitational model (37) at small
$r\rightarrow 0$ takes the form

$\bar{\omega}_{l}(r)\rightarrow\alpha_{1}r^{S_{+}}+\alpha_{2}r^{S_{-}}\,,$
(80)

such that $\alpha_{1}$ and $\alpha_{2}$ are constants of integration and
$S_{\pm}=-\frac{3}{2}\pm\frac{1}{2}\sqrt{9+4\left(l(l+1)-2\right)}$. The
requirement of geometrical regularity at the origin of the star implies that
$\alpha_{2}=0$ and thus $\bar{\omega}^{\prime}_{l}(0)=0$. Using (37) and
replacing the equations (38)-(39) in (79) we finally obtain

$\bar{\omega}^{\prime\prime}_{l}=\left[\frac{-4}{r}+\frac{rh\left(\rho+P\right)+r\phi^{\prime
2}-4\xi\sqrt{h}\phi\phi^{\prime}}{2\left(\kappa+\xi\phi^{2}\right)}\right]\bar{\omega}^{\prime}_{l}-\left[\frac{h(l(l+1)-2)}{r^{2}}+\frac{4\xi\sqrt{h}\phi\phi^{\prime}-4rh(\rho+P)}{r\left(\kappa+\xi\phi^{2}\right)}\right]\bar{\omega}_{l}+\frac{4\xi\Omega\sqrt{h}\phi\phi^{\prime}}{r\left(\kappa+\xi\phi^{2}\right)}\,.$
(81)

Subsequently, one can also extract the initial behavior of $\bar{\omega}(r)$
near the center by following a similar analysis as (48)-(51) for the potential
(65)

$\bar{\omega}_{l}=\bar{\omega}_{l,c}+\frac{6\,\bar{\omega}_{l,c}\left(\rho_{c}+P_{c}\right)+2\xi\left(\Omega-\bar{\omega}_{l,c}\right)\left(-\mu^{2}+\lambda\phi_{c}^{2}\right)\phi_{c}^{2}}{15\left(\kappa+\xi\phi_{c}^{2}\right)}r^{2}\,.$
(82)

As you move farther from the star $r\rightarrow\infty$ taking into account the
sub-horizon limit and the weak field approximation, the spacetime becomes
almost flat and all the values of $\bar{\omega}_{l}$ vanish for $l\geq 2$ so
$\bar{\omega}$ depends only on $r$ and converges toward $\Omega$ at infinity
hartle1967slowly . By setting $l=1$ the solution of the radial function
$\bar{\omega}$ becomes

$\bar{\omega}=\Omega-\frac{2J}{r^{3}}\,,$ (83)

where $J$ being a constant of integration and corresponds to the angular
momentum of the star which is related to the relativistic moment of inertia
for slowly rotating systems through the relation $I=J/\Omega$. Using equation
(77) for $l=1$ one can find the integral expression for the moment of inertia

$I=\frac{1}{3}\int_{0}^{r_{s}}r^{3}\frac{F^{\prime}_{T}}{F_{T}}\left[\sqrt{h}(1-\bar{\omega})-r\bar{\omega}^{\prime}\right]dr+\frac{1}{3}\int_{0}^{r_{s}}\frac{r^{4}}{F_{T}}\sqrt{\frac{h}{f}}\left(\rho+P\right)\bar{\omega}\,dr\,.$
(84)

However in our numerical procedure we prefer to use a double shooting method
to determine the moment of inertia in terms of the mass of the star taking
into account the boundary condition at the radius of the star

$\displaystyle\bar{\omega}(r_{s})=\Omega\left(1-\frac{2I}{r_{s}^{3}}\right)\,,\qquad\frac{d\bar{\omega}(r)}{dr}\biggl{|}_{r=r_{s}}=\frac{6\Omega
I}{r_{s}^{4}}\,.$ (85)

Due to quadratic dependence on radius, the moment of inertia typically
increases in a nearly linear manner as the mass of the NS increases,
particularly for smaller masses while it decreases rapidly once the maximum
mass is approached breu2016maximum . The presence of the scalar fields
influence the behavior of the moment of inertia at larger distances and
therefore increase the rotation speed of the star in compare to GR as it shown
in Fig. (5).

Figure 5: The moment of inertia-mass relations using four EoSs and the same
parameters of Fig. (3)

We end our study by exploring the relations of the normalized moment of
inertia $\widetilde{I}=I/Mr_{s}^{2}$ and the dimensionless moment of inertia
$\overline{I}=I/M^{3}$ in terms of the stellar compactness of the star
$C=M/r_{s}$ whose forms can be written as

$\displaystyle\widetilde{I}$ $\displaystyle=$ $\displaystyle
a_{0}+a_{1}C+a_{4}C^{4}\,,$ (86) $\displaystyle\overline{I}$ $\displaystyle=$
$\displaystyle b_{1}C^{-1}+b_{2}C^{-2}+b_{3}C^{-3}+b_{4}C^{-4}\,,$ (87)

the constants $a_{i}$ for $i=0,1,4$ and $b_{j}$ for $i=1,2,3,4$ are the
fitting coefficients and they are evaluated numerically in Tables (2) and (3)
respectively. The numerical results in Fig. (6) are the plots of the best fit
of $\widetilde{I}$ and $\overline{I}$ for 03 different values of the coupling
parameter $\xi={0,0.05,0.25}$ where we notice small corrections of
$\widetilde{I}$ with the linear fit and the decrease of $\overline{I}$ as we
increase the compactness. These results show an excellent EoS universality for
lower and higher compactness while the values of the reduced $\chi^{2}_{red}$
in the last columns of Tabs. (2) and (3) are in the same orders of paper
staykov2016moment . In the lower panels we measured the relative deviations
$\left|1-I/I_{\text{fit}}\right|$ showing that the errors are below $10\%$ for
the normalized moment of inertia (86) and the dimensionless moment of inertia
(87) relations and they decrease to about $2\%$ for compactness $C=0.3$.

Figure 6: The different normalizations of the moment of inertia in terms of the stellar compactness in slow rotation approximation for the model (37). In the lower panel we plot the polynomial fits and the data $\left|1-I/I_{\text{fit}}\right|$ | $a_{0}$ | $a_{1}$ | $a_{4}$ | $\chi^{2}_{red}$
---|---|---|---|---
GR | $0.200\pm 0.002$ | $0.884\pm 0.013$ | $0.401\pm 0.218$ | $3.131\times 10^{-5}$
$\xi=0.05$ | $0.206\pm 0.002$ | $0.911\pm 0.013$ | $0.505\pm 0.219$ | $3.015\times 10^{-5}$
$\xi=0.25$ | $0.242\pm 0.002$ | $1.009\pm 0.016$ | $0.277\pm 0.308$ | $3.192\times 10^{-5}$
Table 2: The coefficients of the universal relations (86) | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $\chi^{2}_{red}$
---|---|---|---|---|---
GR | $0.843\pm 0.017$ | $0.215\pm 0.008$ | $0.0009\pm 0.001$ | $-0.0002\pm 0.00004$ | $1.403\times 10^{-2}$
$\xi=0.05$ | $0.974\pm 0.021$ | $0.171\pm 0.009$ | $0.007\pm 0.001$ | $-0.0004\pm 0.00005$ | $1.918\times 10^{-2}$
$\xi=0.25$ | $1.049\pm 0.025$ | $0.211\pm 0.011$ | $0.007\pm 0.001$ | $-0.0005\pm 0.00005$ | $2.407\times 10^{-2}$
Table 3: The coefficients of the universal relations (87)

##  VI Concluding Remarks

In the presence of a nonminimally coupled scalar field with a potential of the
form (65) we have studied the NS configuration for spherically symmetric and
static spacetime. We have modeled the stellar matter inside the star by a
perfect fluid with a set of different equations of state: SLy, BSk19, BSk20
and BSk21. The study was carried out to investigate the mass-radius relations
of the star where the numerical results were satisfactory and remained within
the range of the established constraints of the GW170817 event and the maximum
mass of the star has reached $2.37M_{\odot}$ with $\xi=0.25$ for BSk21 which
is a great achievement in the view of the recent observation of the pulsar PSR
J0952-0607.

Another intriguing aspect is that the increase of the coupling parameter $\xi$
is accompanied by an increase in the maximum mass of the star until it reaches
a specific threshold $\xi=0.25$ beyond which the maximum mass starts to
decrease. We illustrated this fact by plotting the maximum mass (as well as
the minimum radius) in terms of the coupling parameter $\xi$ in Fig. (4).
Theoretically, we have found that the spacetime outside the star is a
Schwarzschild de-Sitter type which is a natural result since we haven’t
consider a vanishing potential at asymptotical infinity. The radial profiles
of the star are shown in Fig.(1) where one can notice that the behavior of the
metric functions is almost Minkowskian outside the star as the contributions
of the mass $\mathcal{M}$ and the effective cosmological constant
$\Lambda_{\text{eff}}$ is transparent in the sub-horizon limit and weak field
approximation.

In the second part, we have extended our analysis to include slow rotation
approximation. For this reason, we have adopted the regular Weitzenböck tetrad
to generate an axial symmetric spacetime and derive the radial rotation
equation. Using the relation between the angular momentum $\bar{\omega}$ and
the moment of inertia $I$ we have plotted the relation between this latter and
the mass of the star. Clearly, the presence of the scalar fields influence the
behavior of the moment of inertia at larger distances in compare to GR and
finally we demonstrated that the universal relations (86) and (87) are indeed
independent on the equations of state.

##  Appendix A The Numerical Equations

The dimensionless torsion and boundary term are given by

$\displaystyle T(s)$ $\displaystyle=$
$\displaystyle-\frac{e^{2s}\left(h-1\right)\left(-f+fh-2f_{,s}\right)}{4\pi
fh^{2}}\,,$ (88) $\displaystyle B(s)$ $\displaystyle=$
$\displaystyle\frac{e^{-2s}\left[2f_{,s}f\left(h_{,s}+2h^{2}-4h\right)+h\left(-f_{,ss}+2f_{,s}^{2}+f_{,s}\right)+4f^{2}\left(h_{,s}+h^{2}-h\right)\right]}{8\pi
f^{2}h^{3}}\,.$ (89)

The modified TOV equations (38), (39) and (43) in terms of the new variable
$s$ become

$\displaystyle\frac{f_{,s}}{f}$ $\displaystyle=$
$\displaystyle\frac{h^{2}-1}{2}+\frac{\pi
e^{2s}h^{2}\left(4P+2\mu^{2}\phi^{2}-\lambda\phi^{4}\right)}{1+\xi\phi^{2}}+\frac{\phi_{,s}^{2}}{4\left(1+\xi\phi^{2}\right)}\,,$
(90) $\displaystyle\frac{h_{,s}}{h}$ $\displaystyle=$
$\displaystyle\frac{1-h^{2}}{2}+\frac{\pi
e^{2s}h^{2}\left(4\rho-2\mu^{2}\phi^{2}+\lambda\phi^{4}\right)}{1+\xi\phi^{2}}-\frac{2\xi\left(h-1\right)\phi_{,s}\phi}{1+\xi\phi^{2}}+\frac{\phi_{,s}^{2}}{4\left(1+\xi\phi^{2}\right)}\,,$
(91) $\displaystyle\phi_{,ss}$ $\displaystyle=$
$\displaystyle-\frac{1}{1+\xi\phi^{2}}[2\xi
h\phi(-(h-1)^{2}(1+\xi\phi^{2})-2\pi
e^{2s}h(h-1)(4P+2\mu^{2}\phi^{2}-\lambda\phi^{4})+4\pi
e^{2s}h\phi^{2}(\mu^{2}-\lambda\phi^{2}))$ (92) $\displaystyle+8\pi
e^{2s}h^{2}\phi(\mu^{2}-\lambda\phi^{2})+h^{2}\phi_{,s}(1+\xi\phi^{2}+2\pi
e^{2s}(2P-2\rho+2\mu^{2}\phi^{2}-\lambda\phi^{4}))+(h-1)\xi\phi\phi_{,s}]\,,$

and the energy density in Eq. (44) can be written as

$\displaystyle\rho_{,s}=-\frac{\left(\rho+P\right)}{4P^{\prime}(\rho)}\left[-2\left(1-h^{2}\right)+\frac{4\pi
e^{2s}h^{2}\left(4P+2\mu^{2}\phi^{2}-\lambda\phi^{4}\right)+\phi_{,s}^{2}}{1+\xi\phi^{2}}\right]\,.$
(93)

The Eq. (81) is written as

$\displaystyle\bar{\omega}_{l,ss}$ $\displaystyle=$
$\displaystyle\frac{1}{2\left(1+\xi\phi^{2}\right)}\\{\left[8\pi
e^{2s}h^{2}\left(\rho+P\right)-4\xi
h\phi_{,s}\phi-6\left(2\xi\phi_{,s}\phi+\xi\phi^{2}+1\right)+\phi_{,s}^{2}\right]\bar{\omega}_{l,s}$
(94) $\displaystyle+32h\left[\pi
e^{2s}h\left(\rho+P\right)-\xi\phi_{,s}\phi\right]\bar{\omega}_{l}+32\xi
h\phi_{,s}\phi\\}\,.$

The initial conditions (68)-(71) become

$\displaystyle f(s)$ $\displaystyle=$ $\displaystyle
f_{c}\left[1+\frac{2e^{2s_{0}}\left(6P_{c}+2\rho_{c}+2\mu^{2}\phi_{c}^{2}-\lambda\phi_{c}^{4}\right)}{3\left(1+\xi\phi_{c}^{2}\right)}\right]\,,$
(95) $\displaystyle h(s)$ $\displaystyle=$ $\displaystyle
h_{c}+\frac{2e^{2s_{0}}\left(4\rho_{c}-2\mu^{2}\phi_{c}^{2}+\lambda\phi_{c}^{4}\right)}{3\left(1+\xi\phi_{c}^{2}\right)}\,,$
(96) $\displaystyle\phi(s)$ $\displaystyle=$
$\displaystyle\phi_{c}\left[1+\frac{4}{3}\pi
e^{2s_{0}}\left(-\mu^{2}+\lambda\phi_{c}^{2}\right)\right]\,,$ (97)
$\displaystyle P(s)$ $\displaystyle=$ $\displaystyle P_{c}+\frac{\pi
e^{2s_{0}}\left(\rho_{c}+P_{c}\right)\left(6P_{c}+2\rho_{c}+2\mu^{2}\phi_{c}^{2}-\lambda\phi_{c}^{4}\right)}{3\left(1+\xi\phi_{c}^{2}\right)}\,,$
(98)

and the Eq. (72) become

$\displaystyle M=\frac{3e^{s}}{16\pi
h^{2}\left(1+\xi\phi^{2}\right)}\left[\phi_{,s}\phi-8\xi\phi\left(h-1\right)+4\pi
e^{2s}h^{2}\left(4\rho-2\mu^{2}\phi^{2}+\lambda\phi^{4}\right)\right]\,.$ (99)

Finally the Eq. (82) is written as

$\displaystyle\bar{\omega}$ $\displaystyle=$
$\displaystyle\bar{\omega}_{c}+\frac{16\pi
e^{2s}}{15\left(1+\xi\phi_{c}^{2}\right)}\left[\xi\phi_{c}^{2}\left(\mu^{2}-\lambda\phi_{c}^{2}\right)\left(\bar{\omega}_{c}-\Omega\right)+3\bar{\omega}_{c}\left(\rho_{c}+P_{c}\right)\right]\,.$
(100)

## References

  * (1) C. M. Will, “The confrontation between general relativity and experiment,” Living reviews in relativity 17 (2014) 1–117, arXiv:1403.7377 [gr-qc].
  * (2) A. G. Riess and et al, “Observational evidence from supernovae for an accelerating universe and a cosmological constant,” The astronomical journal 116 no. 3, (1998) 1009, arXiv:astro-ph/9805201.
  * (3) S. Perlmutter and et al, “Measurements of $\mathrm{\Omega}$ and $\mathrm{\Lambda}$ from 42 high-redshift supernovae,” The Astrophysical Journal 517 no. 2, (1999) 565, arXiv:astro-ph/9812133.
  * (4) D. N. Spergel and et al, “First-year Wilkinson Microwave Anisotropy Probe (WMAP)* observations: determination of cosmological parameters,” The Astrophysical Journal Supplement Series 148 no. 1, (2003) 175, arXiv:astro-ph/0302209.
  * (5) B. P. Abbott and et al, “GW150914: The Advanced LIGO detectors in the era of first discoveries,” Physical review letters 116 no. 13, (2016) 131103, arXiv:1602.03838 [gr-qc].
  * (6) B. P. Abbott and et al, “Observation of gravitational waves from a binary black hole merger,” Physical review letters 116 no. 6, (2016) 061102, arXiv:1602.03837 [gr-qc].
  * (7) B. P. Abbott and et al, “GW170817: observation of gravitational waves from a binary neutron star inspiral,” Physical review letters 119 no. 16, (2017) 161101, arXiv:1710.05832 [gr-qc].
  * (8) LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL Collaboration, B. P. Abbott et al., “Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A,” Astrophys. J. Lett. 848 no. 2, (2017) L13, arXiv:1710.05834.
  * (9) B. D. Metzger, “Welcome to the multi-messenger era! lessons from a neutron star merger and the landscape ahead,” arXiv:1710.05931 [gr-qc].
  * (10) G. J. Olmo, D. Rubiera-Garcia, and A. Wojnar, “Stellar structure models in modified theories of gravity: Lessons and challenges,” Physics Reports 876 (2020) 1–75, arXiv:1912.05202 [gr-qc].
  * (11) T. P. Sotiriou and V. Faraoni, “$f(R)$ theories of gravity,” Reviews of Modern Physics 82 no. 1, (2010) 451, arXiv:0805.1726 [gr-qc].
  * (12) A. De Felice and S. Tsujikawa, “$f(R)$ theories,” Living Reviews in Relativity 13 (2010) 1–161, arXiv:1002.4928 [gr-qc].
  * (13) A. Cooney, S. DeDeo, and D. Psaltis, “Neutron stars in $f(R)$ gravity with perturbative constraints,” Physical Review D 82 no. 6, (2010) 064033, arXiv:0910.5480 [astro-ph.HE].
  * (14) S. Arapoğlu, C. Deliduman, and K. Y. Ekşi, “Constraints on perturbative $f(R)$ gravity via neutron stars,” Journal of Cosmology and Astroparticle Physics 2011 no. 07, (2011) 020, arXiv:1003.3179 [gr-qc].
  * (15) M. Orellana, F. García, F. A. Teppa Pannia, and G. E. Romero, “Structure of neutron stars in R-squared gravity,” General Relativity and Gravitation 45 no. 4, (2013) 771–783, arXiv:1301.5189 [astro-ph.CO].
  * (16) A. V. Astashenok, S. Capozziello, and S. D. Odintsov, “Further stable neutron star models from $f(R)$ gravity,” Journal of Cosmology and Astroparticle Physics 2013 no. 12, (2013) 040, arXiv:1309.1978 [gr-qc].
  * (17) M. A. Resco, Á. de la Cruz-Dombriz, F. J. L. Estrada, and V. Z. Castrillo, “On neutron stars in $f(R)$ theories: Small radii, large masses and large energy emitted in a merger,” Physics of the dark universe 13 (2016) 147–161, arXiv:1602.03880 [gr-qc].
  * (18) A. Ganguly, R. Gannouji, R. Goswami, and S. Ray, “Neutron stars in the starobinsky model,” Physical Review D 89 no. 6, (2014) 064019, arXiv:1309.3279 [gr-qc].
  * (19) S. S. Yazadjiev, D. D. Doneva, K. D. Kokkotas, and K. V. Staykov, “Non-perturbative and self-consistent models of neutron stars in R-squared gravity,” Journal of Cosmology and Astroparticle Physics 2014 no. 06, (2014) 003, arXiv:1402.4469 [gr-qc].
  * (20) S. Capozziello, M. De Laurentis, R. Farinelli, and S. D. Odintsov, “Mass-radius relation for neutron stars in $f(R)$ gravity,” Physical Review D 93 no. 2, (2016) 023501, arXiv:1509.04163 [gr-qc].
  * (21) W.-X. Feng, C.-Q. Geng, W.-F. Kao, and L.-W. Luo, “Equation-of-state of neutron stars with junction conditions in the starobinsky model,” International Journal of Modern Physics D 27 no. 01, (2018) 1750186, arXiv:1702.05936 [gr-qc].
  * (22) A. V. Astashenok, S. D. Odintsov, and A. De la Cruz-Dombriz, “The realistic models of relativistic stars in $f(R)$$=$ $R$$+$ $\alpha$$R^{2}$ gravity,” Classical and Quantum Gravity 34 no. 20, (2017) 205008, arXiv:1704.08311 [gr-qc].
  * (23) I. Quiros, “Selected topics in scalar–tensor theories and beyond,” International Journal of Modern Physics D 28 no. 07, (2019) 1930012, arXiv:1901.08690 [gr-qc].
  * (24) E. Babichev and D. Langlois, “Relativistic stars in $f(R)$ and scalar-tensor theories,” Physical Review D 81 no. 12, (2010) 124051, arXiv:0911.1297 [gr-qc].
  * (25) H. Sotani, “Slowly rotating relativistic stars in scalar-tensor gravity,” Physical Review D 86 no. 12, (2012) 124036, arXiv:1211.6986 [astro-ph.HE].
  * (26) M. Minamitsuji and H. O. Silva, “Relativistic stars in scalar-tensor theories with disformal coupling,” Physical Review D 93 no. 12, (2016) 124041, arXiv:1604.07742 [gr-qc].
  * (27) R. Kase and S. Tsujikawa, “Neutron stars in $f(R)$ gravity and scalar-tensor theories,” Journal of Cosmology and Astroparticle Physics 2019 no. 09, (2019) 054, arXiv:1906.08954 [gr-qc].
  * (28) A. S. Arapoğlu, K. Y. Ekşi, and A. E. Yükselci, “Neutron star structure in the presence of nonminimally coupled scalar fields,” Physical Review D 99 no. 6, (2019) 064055, arXiv:1903.00391 [gr-qc].
  * (29) A. S. Arapoğlu and A. E. Yükselci, “A neutron star with a strange quark star-like mass-radius relation,” Physics of the Dark Universe 28 (2020) 100488, arXiv:2002.03025 [gr-qc].
  * (30) S. D. Odintsov and V. Oikonomou, “Neutron stars in scalar–tensor gravity with quartic order scalar potential,” Annals of Physics 440 (2022) 168839, arXiv:2104.01982 [gr-qc].
  * (31) T. Kobayashi, “Horndeski theory and beyond: a review,” Reports on Progress in Physics 82 no. 8, (2019) 086901, arXiv:1901.07183 [gr-qc].
  * (32) D. Langlois and K. Noui, “Degenerate higher derivative theories beyond horndeski: evading the ostrogradski instability,” Journal of Cosmology and Astroparticle Physics 2016 no. 02, (2016) 034, arXiv:1510.06930 [gr-qc].
  * (33) A. Cisterna, T. Delsate, L. Ducobu, and M. Rinaldi, “Slowly rotating neutron stars in the nonminimal derivative coupling sector of horndeski gravity,” Physical Review D 93 no. 8, (2016) 084046, arXiv:1602.06939 [gr-qc].
  * (34) H. Boumaza and D. Langlois, “Neutron stars in degenerate higher-order scalar-tensor theories,” Physical Review D 106 no. 8, (2022) 084053, arXiv:2207.13624 [gr-qc].
  * (35) J. B. Hartle, “Slowly rotating relativistic stars. I. equations of structure,” The Astrophysical Journal 150 (1967) 1005.
  * (36) S. S. Yazadjiev, D. D. Doneva, and D. Popchev, “Slowly rotating neutron stars in scalar-tensor theories with a massive scalar field,” The European Physical Journal C (2016) , arXiv:1602.04766 [gr-qc].
  * (37) K. V. Staykov, D. Popchev, D. D. Doneva, and S. S. Yazadjiev, “Static and slowly rotating neutron stars in scalar–tensor theory with self-interacting massive scalar field,” The European Physical Journal C 78 (2018) 1–8, arXiv:1805.07818 [gr-qc].
  * (38) T. Damour and G. Esposito-Farèse, “Nonperturbative strong-field effects in tensor-scalar theories of gravitation,” Phys. Rev. Lett. 70 (Apr, 1993) 2220–2223. https://link.aps.org/doi/10.1103/PhysRevLett.70.2220.
  * (39) M. Salgado, D. Sudarsky, and U. Nucamendi, “Spontaneous scalarization,” Physical Review D 58 no. 12, (1998) 124003, arXiv:gr-qc/9806070.
  * (40) P. Pani and E. Berti, “Slowly rotating neutron stars in scalar-tensor theories,” Physical Review D 90 no. 2, (2014) 024025, arXiv:1405.4547 [gr-qc].
  * (41) K. Yagi and M. Stepniczka, “Neutron stars in scalar-tensor theories: Analytic scalar charges and universal relations,” Physical Review D 104 no. 4, (2021) 044017, arXiv:2105.01614 [gr-qc].
  * (42) Z. Hu, Y. Gao, R. Xu, and L. Shao, “Scalarized neutron stars in massive scalar-tensor gravity: X-ray pulsars and tidal deformability,” Physical Review D 104 no. 10, (2021) 104014, arXiv:2109.13453 [gr-qc].
  * (43) G. Ventagli, G. Antoniou, A. Lehébel, and T. P. Sotiriou, “Neutron star scalarization with gauss-bonnet and ricci scalar couplings,” Physical Review D 104 no. 12, (2021) 124078, arXiv:2111.03644 [gr-qc].
  * (44) R. Kase, M. Minamitsuji, and S. Tsujikawa, “Neutron stars with a generalized proca hair and spontaneous vectorization,” Physical Review D 102 no. 2, (2020) 024067, arXiv:2001.10701v2 [gr-qc].
  * (45) D. D. Doneva and S. S. Yazadjiev, “Topological neutron stars in tensor-multi-scalar theories of gravity,” Physical Review D 101 no. 6, (2020) 064072, arXiv:1911.06908 [gr-qc].
  * (46) C. Breu and L. Rezzolla, “Maximum mass, moment of inertia and compactness of relativistic stars,” Monthly Notices of the Royal Astronomical Society 459 no. 1, (2016) 646–656, arXiv:1601.06083 [gr-qc].
  * (47) K. V. Staykov, D. D. Doneva, and S. S. Yazadjiev, “Moment-of-inertia–compactness universal relations in scalar-tensor theories and ${R}^{2}$ gravity,” Physical Review D 93 no. 8, (2016) 084010, arXiv:1602.00504 [gr-qc].
  * (48) D. Popchev, K. V. Staykov, D. D. Doneva, and S. S. Yazadjiev, “Moment of inertia–mass universal relations for neutron stars in scalar-tensor theory with self-interacting massive scalar field,” The European Physical Journal C 79 (2019) 1–11, arXiv:1812.00347 [gr-qc].
  * (49) D. D. Doneva, S. S. Yazadjiev, N. Stergioulas, and K. D. Kokkotas, “Rapidly rotating neutron stars in scalar-tensor theories of gravity,” Physical Review D 88 no. 8, (2013) 084060, arXiv:1309.0605 [gr-qc].
  * (50) S. S. Yazadjiev, D. D. Doneva, and K. D. Kokkotas, “Rapidly rotating neutron stars in R-squared gravity,” Physical Review D 91 no. 8, (2015) 084018, arXiv:1501.04591 [gr-qc].
  * (51) S. Bahamonde, K. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. Said, J. Mifsud, and E. Di Valentino, “Teleparallel gravity: from theory to cosmology,” Reports on Progress in Physics (2021) , arXiv:2106.13793 [gr-qc].
  * (52) R. Aldrovandi and J. G. Pereira, Teleparallel gravity: an introduction, vol. 173.  Springer Science & Business Media, 2012.
  * (53) M. Krššák, R. Van Den Hoogen, J. Pereira, C. Böhmer, and A. Coley, “Teleparallel theories of gravity: illuminating a fully invariant approach,” Classical and Quantum Gravity 36 no. 18, (2019) 183001, arXiv:1810.12932 [gr-qc].
  * (54) Y.-F. Cai, S. Capozziello, M. De Laurentis, and E. N. Saridakis, “$f(T)$ teleparallel gravity and cosmology,” Reports on Progress in Physics 79 no. 10, (2016) 106901, arXiv:1511.07586 [gr-qc].
  * (55) C.-Q. Geng, C.-C. Lee, E. N. Saridakis, and Y.-P. Wu, ““Teleparallel” dark energy,” Physics Letters B 704 no. 5, (2011) 384–387, arXiv:1109.1092 [gr-qc].
  * (56) D. Horvat, S. Ilijić, A. Kirin, and Z. Narančić, “Nonminimally coupled scalar field in teleparallel gravity: boson stars,” Classical and Quantum Gravity 32 no. 3, (2015) 035023, arXiv:1407.2067 [gr-qc].
  * (57) M. Hohmann, “Scalar-torsion theories of gravity. I. general formalism and conformal transformations,” Physical Review D 98 no. 6, (2018) 064002, arXiv:1801.06528 [gr-qc].
  * (58) M. Hohmann and C. Pfeifer, “Scalar-torsion theories of gravity. II. $L(T,X,Y,\phi)$ theory,” _Physical Review D_ 98 no. 6, (2018) 064003, arXiv:1801.06536 [gr-qc].
  * (59) M. Hohmann, “Scalar-torsion theories of gravity. III. analogue of scalar-tensor gravity and conformal invariants,” Physical Review D 98 no. 6, (2018) 064004, arXiv:1801.06531 [gr-qc].
  * (60) B. Li, T. P. Sotiriou, and J. D. Barrow, “$f(T)$ gravity and local lorentz invariance,” Physical Review D 83 no. 6, (2011) 064035, arXiv:1010.1041 [gr-qc].
  * (61) M. Krššák and E. N. Saridakis, “The covariant formulation of $f(T)$ gravity,” Classical and Quantum Gravity 33 no. 11, (2016) 115009, arXiv:1510.08432 [gr-qc].
  * (62) C. G. Boehmer, A. Mussa, and N. Tamanini, “Existence of relativistic stars in $f(T)$ gravity,” Classical and Quantum Gravity 28 no. 24, (2011) 245020, arXiv:1107.4455 [gr-qc].
  * (63) R.-H. Lin, X.-N. Chen, and X.-H. Zhai, “Realistic neutron star models in $f(T)$ gravity,” The European Physical Journal C 82 no. 4, (2022) 308, arXiv:2109.00191 [gr-qc].
  * (64) H. Fortes and J. Araujo, “Solving Tolman–Oppenheimer–Volkoff equations in $f(T)$ gravity: a novel approach,” Classical and Quantum Gravity 39 no. 24, (2022) 245017, arXiv:2105.04473 [gr-qc].
  * (65) H. Boumaza, “Slowly rotating neutron stars in scalar torsion theory,” The European Physical Journal C 81 no. 5, (2021) 448.
  * (66) K. Hinterbichler and J. Khoury, “Screening long-range forces through local symmetry restoration,” Physical review letters 104 no. 23, (2010) 231301, arXiv:1001.4525 [hep-th].
  * (67) R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, “PSR J0952- 0607: The fastest and heaviest known galactic neutron star,” The Astrophysical Journal Letters 934 no. 2, (2022) L18, arXiv:2207.05124 [astro-ph.HE].
  * (68) M. Hohmann, L. Järv, M. Krššák, and C. Pfeifer, “Modified teleparallel theories of gravity in symmetric spacetimes,” Physical Review D 100 no. 8, (2019) 084002, arXiv:1901.05472 [gr-qc].
  * (69) S. Bahamonde, J. G. Valcarcel, L. Järv, and C. Pfeifer, “Exploring axial symmetry in modified teleparallel gravity,” Physical Review D 103 no. 4, (2021) 044058, arXiv:2012.09193 [gr-qc].
  * (70) N. Tamanini and C. G. Boehmer, “Good and bad tetrads in $f(T)$ gravity,” Physical Review D 86 no. 4, (2012) 044009.
  * (71) M. Krššák and J. Pereira, “Spin connection and renormalization of teleparallel action,” The European Physical Journal C 75 (2015) 1–8, arXiv:1504.07683 [gr-qc].
  * (72) T. P. Sotiriou and V. Faraoni, “Black holes in scalar-tensor gravity,” Physical Review Letters 108 no. 8, (2012) 081103, arXiv:1109.6324 [gr-qc].
  * (73) E. Babichev, K. Koyama, D. Langlois, R. Saito, and J. Sakstein, “Relativistic stars in beyond horndeski theories,” Classical and Quantum Gravity 33 no. 23, (2016) 235014, arXiv:1606.06627 [gr-qc].
  * (74) P. Haensel and A. Y. Potekhin, “Analytical representations of unified equations of state of neutron-star matter,” Astronomy & Astrophysics 428 no. 1, (2004) 191–197, arXiv:astro-ph/0408324 [astro-ph].
  * (75) A. Potekhin, A. Fantina, N. Chamel, J. Pearson, and S. Goriely, “Analytical representations of unified equations of state for neutron-star matter,” Astronomy & Astrophysics 560 (2013) A48, arXiv:1310.0049 [gr-qc].
  * (76) B. F. de Aguiar, R. F. Mendes, and F. T. Falciano, “Neutron stars in the symmetron model,” Universe 8 no. 1, (2021) 6, arXiv:2112.03823 [gr-qc].