# Pre-main sequence evolution of low-mass stars in Eddington-inspired Born-
Infeld gravity

Merce Guerrero<EMAIL_ADDRESS>Departamento de Física Teórica and IPARCOS,
Universidad Complutense de Madrid, E-28040 Madrid, Spain  Diego Rubiera-Garcia
<EMAIL_ADDRESS>Departamento de Física Teórica and IPARCOS, Universidad
Complutense de Madrid, E-28040 Madrid, Spain  Aneta Wojnar
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Institute of
Physics, University of Tartu, W. Ostwaldi 1, 50411 Tartu, Estonia

###### Abstract

We study three aspects of the early-evolutionary phases in low-mass stars
within Eddington-inspired Born-Infeld (EiBI) gravity, a viable extension of
General Relativity. These aspects are concerned with the Hayashi tracks (i.e.
the effective temperature-luminosity relation); the minimum mass required to
belong to the main sequence; and the maximum mass allowed for a fully
convective star within the main sequence. In all cases we find a dependence of
these quantities not only on the theory’s parameter, but also on the star’s
central density, a feature previously found in Palatini $f(R)$ gravity. Using
this, we investigate the evolution of these quantities with the (sign of the)
EiBI parameter, finding a shift in the Hayashi tracks in opposite directions
in the positive/negative branches of it, and an increase (decrease) for
positive (negative) parameter in the two masses above. We use these results to
ellaborate on the chances to seek for traces of new physics in low-mass stars
within this theory.

##  I Introduction

Despite the many tests that Einstein’s General Theory of Relativity (GR) has
successfully passed Will:2014kxa , over the last decades a plethora of
modified theories of gravity has been proposed in order to address its
shortcomings Capozziello:2011et ; Nojiri:2017ncd . These include the yet
undetected dark energy/matter sources needed for the consistence of the
cosmological concordance model Bull:2015stt , or the existence of space-time
singularities at high-energy scales Senovilla:2014gza , such as the one at the
center of black holes and the Big Bang singularity. Attempts to modify GR must
come along with the necessity to comply with its well tested weak-field limit,
while deviations with respect to its predictions must be searched in those
domains where the strength of the gravitational interaction grows large
enough, for instance, via gravitational waves out of binary mergers
Ezquiaga:2020dao , gravitational lensing and shadows
EventHorizonTelescope:2020qrl , or in the structure of neutron stars
Olmo:2019flu .

From an astrophysical point of view, neutron stars are suitable objects to
test modified gravity via the opportunity in the combination of
electromagnetic radiation and gravitational waves that the newly born field of
multimessenger astronomy offers LIGOScientific:2017vwq . Among the open
problems in this field, one can mention the degeneracy of the mass-radius
relations due to the fact that the equation of state at supranuclear densities
is unknown Camenzind , the theoretical difficulties to meet the observations
of neutron stars above two solar masses Linares:2018ppq ; Antoniadis:2013pzd ,
or the so-called mass gap problem, namely, the existence of objects above the
neutron star mass limit but below the lightest black hole mass, as manifested
in the observation of gravitational waves from the merging of two objects with
2.6 and 23 $M_{\odot}$, respectively LIGOScientific:2020zkf .

We consider here the alternative path of focusing on other astrophysical
objects for which their internal structure and equation of state are better
known. Even though in such objects the gravitational interaction is less
strong as compared with neutron stars and black holes, modified gravity
effects are able to induce extra terms to the Poisson equation, see e.g.
Sakstein:2018fwz . This leads to a different stellar structure whose
associated macroscopic features can be tested. As examples, we highlight the
predictions for the time scales and effective luminosities of both main
sequence stars Chowdhury:2020wfr  and sub-stellar objects such as as brown
dwarfs Benito:2021ywe  and giant gaseous planets Wojnar:2021xbr , lithium
abundance and age determination for white dwarfs and low mass stars (LMS)
Wojnar:2020frr ; Wojnar:2020ckw ; lupa , the minimum hydrogen burning mass for
high-mass brown dwarfs to belong to the main sequence, Sakstein:2015zoa ;
Sakstein:2015aac ; Crisostomi:2019yfo ; Olmo:2019qsj ; Rosyadi:2019hdb , early
evolutionary tracks of LMS Wojnar:2020txr , and even tests with exoplanets
olek1 ; olek2 ; olek3  are at our disposal.

The main aim of the present work is to highlight the most important phases of
the early evolution of LMS, within a suitable extension of GR dubbed as
Eddington-inspired Born-Infeld gravity (EiBI) Banados:2010ix . EiBI gravity
belongs to the so-called Ricci Based Gravities (RBG), a family of viable
gravitational extensions of GR constructed in terms of scalars out of
contractions of the metric with the (symmetric part of the) Ricci tensor. RBGs
are formulated à la Palatini, that is, taking metric and affine connection as
independent entities. This fact allows them to yield second-order field
equations that do not propagate additional degrees of freedom beyond the two
polarizations of the gravitational field. This acts as a safeguard of (most
of) these theories against getting into conflict with solar system tests and
gravitational wave observations so far, while at the same time offering a
workable framework to extract new gravitational physics in the strong-field
regime.

The relative simplicity of the stellar structure equations of EiBI gravity is
a key feature that has allowed the community to scan its predictions as
compared to GR expectations within different types of stars, see e.g.
Pani:2011mg ; Delsate:2012ky ; Casanellas:2011kf ; Avelino:2012ge ;
Banerjee:2017uwz ; Wibisono:2017dkt ; Danarianto:2019mxf ; Rosyadi:2019hdb ;
Banerjee:2021auy . In particular, on its non-relativistic regime, this theory
leads to a modified Poisson equation with a single extra free parameter (for a
full account of the theory and its phenomenology we refer the reader to
BeltranJimenez:2017doy ). In our analysis of the pre-main sequence evolution
of LMS within it, we shall mainly focus on the paths followed by the
contracting star, represented by Hayashi tracks hayashi  on the Hertzsprung-
Russell (HR) diagram, and the associated limiting masses by the hydrogen
burning and the development of a fully convective core, respectively, at the
gateway of the main sequence. To simplify our analysis we shall disregard the
deuterium burning process which happens during the pre-main sequence phase and
in massive brown dwarfs, since the energy generated by this process is
significantly smaller than the one of the hydrogen ignition. Another
simplification of our analysis lies on the fact that in order to properly
incorporate lithium burning in the low-mass stars or cooling process of a
brown dwarf object, one needs to use a more realistic model of the electron
degeneracy than the one used in this work, while the choice of the opacities
is always subject to discussion.

This work is organized as follows: in Sec. II, we introduce the main
ingredients of EiBI gravity, and work out its non-relativistic equations until
arriving to the generalized Lane-Emden equation for a polytropic equation of
state, and set our simplified photospheric model for it and the convenctive
instability criterion. In Sec. III, three main aspects of the early evolution
of LMS within this theory are discussed: i) its Hayashi tracks, namely, the
effective temperature vs luminosity evolution, ii) the minimum main sequence
mass, namely, the minimum required mass for a star to stably burn sufficient
hydrogen, and iii) the maximum fully convective mass. In all these three cases
we discuss the modifications in both the positive and negative branches of the
EiBI parameter (which appears unavoidably entangled with the star’s central
density, a common feature to RBGs), with the former being the most succulent
from a physical point of view. Finally, Sec. IV contains some closing
thoughts.

##  II Model and basic equations

###  II.1 Action and field equations of EiBI gravity

The action of EiBI gravity is given by

$\displaystyle\mathcal{S}_{EiBI}$ $\displaystyle=$
$\displaystyle\frac{1}{\kappa^{2}\epsilon}\int
d^{4}x\left[\sqrt{-|g_{\mu\nu}+\epsilon
R_{\mu\nu}(\Gamma)|}-\lambda\sqrt{-g}\right]$ (1) $\displaystyle+$
$\displaystyle\mathcal{S}_{m}(g_{\mu\nu},\psi_{m})$

where $\kappa^{2}=8\pi G/c^{4}$ is Newton’s constant, $\epsilon$ is the
(length squared) EiBI parameter, $\lambda$ is related to the asymptotic
character of the solutions (from now on we fix $\lambda=1$ to deal with
asymptotically flat solutions), vertical vars denote a determinant, while we
reserve the notation of $g$ for the determinant of the space-time metric
$g_{\mu\nu}$, the latter being a priori independent of the affine connection
$\Gamma\equiv\Gamma_{\mu\nu}^{\lambda}$ of the (symmetric) Ricci tensor
$R_{\mu\nu}(\Gamma)$. As for the matter fields, collectively labeled by
$\psi_{m}$, they only couple to the space-time metric, but not to the
independent connection $\Gamma$, which is consistent as long as no fermionic
fields are present Afonso:2017bxr .

The field equations of action (1) are obtained as

$\frac{\sqrt{-q}}{\sqrt{-g}}q^{\mu\nu}-g^{\mu\nu}=\kappa^{2}\epsilon
T^{\mu\nu}\ ,$ (2)

where the rank-two tensor defined via

$q_{\mu\nu}\equiv g_{\mu\nu}+\epsilon\,R_{\mu\nu}(\Gamma)\ ,$ (3)

is the metric compatible with the independent connection, while
$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta\mathcal{L}_{m}}{\delta
g^{\mu\nu}}$ is the stress-energy tensor of the matter. The relation (3) can
be rewritten in the more seductive form

$q_{\mu\nu}=g_{\mu\alpha}{\Omega^{\alpha}}_{\nu}\ ,$ (4)

in such a way that the deformation matrix ${\Omega^{\mu}}_{\alpha}$ is
determined by the algebraic relation

$|\Omega|^{1/2}({\Omega^{\mu}}_{\nu})^{-1}=\delta_{\mu}^{\nu}-\kappa^{2}\epsilon\,{T^{\mu}}_{\nu}\
.$ (5)

With the help of the definitions above, the metric field equations (2) can be
recast as

${R^{\mu}}_{\nu}(q)=\frac{1}{|\Omega|^{1/2}}\left(\mathcal{L}_{G}\delta^{\mu}_{\nu}+\kappa^{2}{T^{\mu}}_{\nu}\right)\
,$ (6)

where the gravitational EiBI Lagrangian can also be expressed in terms of the
deformation matrix as
$\mathcal{L}_{G}=\tfrac{|\Omega|^{1/2}}{\epsilon\kappa^{2}}-1$. The
representation (6) of the EiBI field equations puts forward the existence of
an Einstein frame for this theory, sourced on its right-hand via additional
(non-dynamical) couplings in the matter fields, a common feature of the whole
RBG family Afonso:2018bpv . In the non-relativistic limit, all these RBG
theories share a common set of extra pieces to the Poisson equation with
theory-dependent coefficients Olmo:2021yul .

###  II.2 Non-relativistic stellar structure equations

Let us thus head to the non-relativistic limit of the field equations (6)
above. To this end, we set the following ansatz for the time-independent
metric:

$ds^{2}=-(1+2\Phi)dt^{2}+(1-2\psi)d\vec{x}d\vec{x}\ ,$ (7)

where $\Phi$ and $\psi$ are only functions of $\vec{x}$. As for the stress-
energy tensor, since the pressure is generally negligible for non-relativistic
stars, we take a relativistic pressureless fluid,
$T^{\mu\nu}=\rho\,u^{\mu}u^{\nu}$, where $\rho$ is the energy density and
$g^{\mu\nu}u_{\mu}u_{\nu}=-1$ a unit time-like vector. From Eq.(5), we can
easily find the components of the corresponding deformation matrix and apply
them to Eq.(4) to get the metric components of $q_{\mu\nu}$. Finally, if one
expands Eq.(3) up to linear order in $\Phi,\,\psi$, $\rho$ and their
derivatives, the modified Poisson equation is found to be Pani:2011mg ;
Olmo:2021yul

$\nabla^{2}\Phi=\dfrac{\kappa^{2}}{2}\rho+\dfrac{\kappa^{2}\epsilon}{4}\nabla^{2}\rho\
,$ (8)

where the second term corresponds to the EiBI correction. In a static,
spherically symmetric space-time, the above equation can be rewritten as

$\dfrac{1}{r^{2}}\dfrac{d}{dr}\left(r^{2}\dfrac{d\Phi}{dr}\right)=\dfrac{\kappa^{2}}{2}\rho+\dfrac{\kappa^{2}\epsilon}{4\,r^{2}}\dfrac{d}{dr}\left(r^{2}\dfrac{d\rho}{dr}\right)\
.$ (9)

Using the hydrostatic equation $\tfrac{d\Phi}{dr}=-\rho^{-1}\tfrac{dP}{dr}$
and integrating it over the radial coordinate $r$, the above equation
transforms into

$\displaystyle\dfrac{dP}{dr}=-\dfrac{\kappa^{2}M(r)\rho}{8\pi
r^{2}}-\dfrac{\kappa^{2}\epsilon\,\rho}{4}\dfrac{d\rho}{dr}\ ,$ (10)

where the mass function $M(r)$ is defined as

$M(r)=\int^{r}_{0}4\pi x^{2}\rho dx\ .$ (11)

The solutions of the hydrostatic equilibrium equations (10) and (11), equipped
with an equation of state [given in Sec. II.3], provide the main ingredient
from the gravitational sector for the internal and external features of a LMS
on its early evolutionary stages. As shall be discussed later, they also
contribute to the description of the boundary region between the star’s
interior and its photosphere, and have an impact on the photospheric
quantities.

###  II.3 Polytropic stars

It is now time to move to the matter description of LMS. Our simplified model
assumes a fully convective interior, from the center up to the photosphere.
Such objects are typically well described by a polytropic equation of state
which in a general case takes the form

$p=K\rho^{\frac{n+1}{n}}\ ,$ (12)

where $n$ is the polytropic index, while $K$ is the degenerate parameter which
carries the information about the microscopic features of the fluid, such as
e.g. electron degeneracy, Coulomb force, and ionization. Let us introduce the
following dimensionless variables

$\rho=\rho_{c}\theta^{n},\;\;\;P=p_{c}\theta^{n+1},\;r=r_{c}\xi\ ,$ (13)

where $\rho_{c}$ and $p_{c}$ are the star’s central density and pressure,
respectively, while $r_{c}$ is defined via the expression

$r_{c}^{2}=\dfrac{2(n+1)p_{c}}{\kappa^{2}\rho_{c}^{2}}=\frac{2(n+1)K\rho_{c}^{1/n-1}}{\kappa^{2}}\
.$ (14)

These variables allow to rewrite the hydrostatic equilibrium equation (9) in a
more suitable form under the modified Lane-Emden equation

$\dfrac{d}{d\xi}\left\\{\xi^{2}\,\dfrac{d\theta}{d\xi}\left[1+\alpha\theta^{n-1}\right]\right\\}=-\xi^{2}\theta^{n}\
.$ (15)

In this equation the EiBI corrections are encapsulated into the single
dimensionless parameter

$\alpha=\dfrac{\epsilon\,n}{2\,r_{c}^{2}}\ ,$ (16)

which depends not only on the polytropic parameters but also on the star’s
central energy density $\rho_{c}$, as given by Eq.(14). This is a general
feature of Palatini theories of gravity (at least for the RBG family), caused
by the particular way the matter fields source the new gravitational dynamics
aw_pol , and strongly departs from what happens in other theories of gravity,
including GR itself. This implies that astrophysical constraints on EiBI
gravity (and on the whole RBG family) parameter requires further information
on the star’s central density, as we shall see later.

The generalized Lane-Emden equation (15) can be used to rewrite the mass
function (11) as well as other relevant stellar features such as radius,
central density and temperature in terms of its solutions:

$\displaystyle M$ $\displaystyle=$ $\displaystyle
4\pi\rho_{c}r_{c}^{3}\omega_{n},$ (17) $\displaystyle R$ $\displaystyle=$
$\displaystyle\gamma_{n}\left(\dfrac{K}{G}\right)^{\frac{n}{3-n}}M^{\frac{1-n}{3-n}}\
,$ (18) $\displaystyle\rho_{c}$ $\displaystyle=$
$\displaystyle\delta_{n}\left(\dfrac{3M}{4\pi R^{3}}\right)\ ,$ (19)
$\displaystyle T$ $\displaystyle=$
$\displaystyle\frac{K\mu}{N_{A}k_{B}}\rho_{c}^{1/n}\theta\ ,$ (20)

where $k_{B}$ is Boltzmann’s constant, $N_{A}$ the Avogadro number and $\mu$
the mean molecular weight, while the central temperature is defined as
$T_{c}=\frac{K\mu}{N_{A}k_{B}}\rho_{c}^{1/n}$. The three remaining constants,
$\omega_{n}$, $\gamma_{n}$, and $\delta_{n}$, come from the evaluation of the
corresponding solution of the generalized Lane-Emden equation (15) at the
star’s surface, $\xi_{R}$, via the expressions

$\displaystyle\omega_{n}$ $\displaystyle=$
$\displaystyle\left[-\xi^{2}\,\dfrac{d\theta}{d\xi}\left(1+\alpha\theta^{n-1}\right)\right]=\left[-\xi^{2}\,\dfrac{d\theta}{d\xi}\right]_{\xi=\xi_{R}},$
(21) $\displaystyle\gamma_{n}$ $\displaystyle=$
$\displaystyle(4\pi)^{\frac{1}{n-3}}(n+1)^{\frac{n}{3-n}}\,\omega_{n}^{\frac{n-1}{3-n}}\xi_{R},$
(22) $\displaystyle\delta_{n}$ $\displaystyle=$
$\displaystyle-\dfrac{\xi_{R}}{3\left[\dfrac{d\theta}{d\xi}\left(1+\alpha\theta^{n-1}\right)\right]}=-\frac{\xi_{R}}{3\dfrac{d\theta}{d\xi}\Big{|}_{\xi=\xi_{R}}},$
(23)

where the last equalities come from applying the boundary condition of the
surface, i.e., $\theta(\xi_{R})=0$.

In general, analytic solutions to the generalized Lane-Emden equation (15) are
not possible, so one has to resort to a numerical resolution procedure. For
the sake of this work, let us take the value of the polytropic index $n=3/2$,
which is the one suitable to describe the convective interior of a LMS. In
such a case, one can approximate the central behaviour of the solution of the
Lane-Emden equation (15) by

$\theta(\xi\approx
0)=1-\dfrac{\xi^{2}}{6(1+\alpha)}\sim\text{exp}\left(-\dfrac{\xi^{2}}{6(1+\alpha)}\right),$
(24)

where the initial conditions $\theta(0)=1$ and $\theta^{\prime}(0)=0$ have
been applied. This result will prove its usefulness later.

###  II.4 Simple photospheric model

As already mentioned, the model described above holds up to the photosphere,
which is the outer, luminous layer that delimits the star. It is formally
defined as the radius for which the so-called optical depth equals the value
2/3, that is (see e.g. hansen )

$\tau(r)=\int_{r_{ph}}^{\infty}\kappa_{op}\,\rho\,dr=\frac{2}{3}\ ,$ (25)

where $\kappa_{op}$ is dubbed as the opacity, a phenomenological quantity
playing a key role in the characterization of the star. Later on, we shall use
various opacity models, depending on the physical features of the material
filling the star. The photosphere is so close to the surface of the star that
its radius, $r_{\text{ph}}$, can be well approximated by the star’s radius,
$R$. The photospheric temperature, often called “effective temperature”, is
the one appearing in the Stefan-Boltzmann equation

$L=4\pi\sigma R^{2}T^{4}_{\text{eff}}\ ,$ (26)

where $\sigma$ is the Stefan-Boltzmann constant and $L$ the luminosity.
Therefore, we shall assume the star to radiate its energy as a black body with
a temperature $T_{\text{eff}}$.

The hydrostatic equilibrium equation (10) can be conveniently rewritten in the
following way

$p^{\prime}=-\rho\,\left(g+\dfrac{\kappa^{2}\epsilon\rho^{\prime}}{4}\right)\
,$ (27)

where primes indicate radial derivatives, while $g$ is the surface gravity
defined as

$g\equiv\dfrac{\kappa^{2}M(r)}{8\pi r^{2}}\sim\dfrac{\kappa^{2}M}{8\pi
R^{2}}=\text{const}\ .$ (28)

Taking up to two derivatives in the equation above and using the definition
(11), one can combine the resulting expressions to find

$\rho^{\prime}=-\dfrac{M(r)}{2\pi r^{4}}\ ,$ (29)

so that Eq.(27) evaluated in the photospheric region becomes

$p^{\prime}_{ph}=-\rho g\,\left(1-\frac{\epsilon}{R^{2}}\right)\ ,$ (30)

where in this equation we have set units $\kappa^{2}=8\pi G$ (and assumed
$c=1$ from now on). This equation can be integrated with the help of (25),
providing the photospheric pressure as

$p_{ph}=\dfrac{2g\left(1-\dfrac{\epsilon}{R^{2}}\right)}{3\kappa_{op}}\ .$
(31)

It is clear now that the most relevant element of the photosphere’s modelling
is its opacity. Depending on the physical conditions, mainly contained within
the pressure and temperature regimes of the considered stages of the stellar
evolution, we shall use different analytical expressions, which approximately
reflect how opaque matter is to the electromagnetic radiation.

###  II.5 Convective instability - modified Schwarzschild criterion

Another crucial information in the description of stellar interiors is how the
energy is transported through different regions of a star. Since our LMS is
modelled by a fully convective sphere enveloped by a radiative photosphere,
one needs a formal criterion encapsulating the physical conditions responsible
for any of those energy transports. This is given by the so-called
Schwarzschild criterion, turning out to be dependent on the underlying theory
of gravity, as shown in Wojnar:2020ckw . Therefore, the heat is transported
via radiative processes when the temperature gradient is smaller than the
adiabatic one

$\nabla_{rad}<\nabla_{ad}\ .$ (32)

For the rest of our setup we shall model the photosphere’s matter as an ideal,
monoatomic gas, for which it can be shown that the adiabatic gradient has a
constant value, $\nabla_{ad}=0.4$ hansen . On the other hand, the radiative
gradient is defined as

$\nabla_{rad}:=\left(\dfrac{d\ln T}{d\ln P}\right)_{rad}.$ (33)

To find its form and dependence on EiBI gravity parameter, we need to analyze
the radiative heat transport equation, which is given by

$\dfrac{\partial T}{\partial
m}=-\dfrac{3}{64\pi^{2}a}\dfrac{\kappa_{rc}l}{r^{4}T^{3}}\ ,$ (34)

where $\kappa_{rc}$ is the radiative or/and conductive opacity, $l$ the local
luminosity, while $a=7.57\times 10^{-15}\frac{erg}{cm^{3}K^{4}}$ represents
the radiation density constant. Combining the above expression with the
modified hydrostatic equilibrium equation (10) differentiated with respect to
the mass, one gets

$\dfrac{\partial T}{\partial P}=\dfrac{3\kappa_{rc}\,l}{16\pi
r^{2}aT^{3}}\left(\dfrac{Gm}{r^{2}}+\dfrac{\kappa^{2}\epsilon\rho^{\prime}}{4}\right)^{-1}\
.$ (35)

Using this result into the definition (33) provides the radiative temperature
gradient for EiBI gravity as

$\nabla_{rad}=\dfrac{3\kappa_{rc}\,l\,p}{16\pi
r^{2}aT^{3}}\left(1-\dfrac{\epsilon}{r^{2}}\right)^{-1}\ ,$ (36)

and therefore, depending on the value of the EiBI parameter $\epsilon$, this
modification has a (des-)stabilizing effect, altering a radiative region
development.

##  III Early evolution of low-mass stars

With the formalism developed in the previous section, we are ready to study
the early evolution of a LMS within EiBI gravity. When the baby proto-star
approaches the main sequence, it is a luminous but otherwise cold (sub-)
stellar object. Similarly as for the later phases, such an object can be
accommodated on the HR diagram; thus, it can be found on the right-hand side
part of the evolutionary diagram, above the main sequence band. The
evolutionary path that it follows is called a Hayashi track hayashi ,
described by a relation between the effective temperature, luminosity, mass,
and metallicity, where the last one is responsible for the shape of the curve.
However, since we are dealing with a toy-model description to understand the
new features brought by the gravitational corrections of the EiBI gravity,
that aspect will not be apparent in our subsequent analysis.

A stellar object will leave its Hayashi track when any of the following
processes happens:

  * •

Radiative core development: Since the luminosity decreases as the baby star
follows the Hayashi track down but the effective temperature remains almost
constant, this means that, from the Stefan-Boltzman law, the star is
contracting. Therefore, it may happen that the star’s interior becomes
radiative, as a consequence of increasing its interior temperature. In such a
situation, it will reach a minimum and follow an almost horizontal line before
getting to the main sequence, moving to higher effective temperatures. This
stage of the early evolution is called a Henyey track henyey ; hen2 ; hen3 ,
and it will not be studied here; however, in Sec. III.3 we shall discuss in
detail the onset of the radiative core development as a boundary condition of
the fully convective star on the main sequence.

  * •

Hydrogen ignition: When the central temperature and pressure increase, the
conditions present in the stellar core can become sufficient to ignite
hydrogen and stop further gravitational contraction. If the process is stable,
in other words, when the energy radiated away through the photosphere is
balanced by energy produced by the hydrogen burning in the core, the star has
evolved to the next stage of the stellar evolution, that is, the main sequence
phase, which we analyze in Sec. III.2.

  * •

Contraction stops at the onset of electronic degeneracy: This process will
happen when none of the above ones takes place - that is, the interior of such
an object is too cold to start hydrogen burning. Apart from the light elements
burnt in the initial phase, those objects do not possess any source of energy
production in their cores and, therefore, they will cool down with time when
electron degeneracy pressure balances the gravitational contraction
Burrows:1992fg . Such objects are called brown dwarfs, and will be discussed
somewhere else in detail.

###  III.1 Hayashi tracks

In what follows, we will now focus on a simple description of the Hayashi
tracks in EiBI gravity. As mentioned before, the objects following this stage
of the evolution are fully convective, and we shall also assume that their
interiors are made of a fully ionized monatomic gas with temperature $T$ and
mean molecular weight $\mu$. In such a situation, the equation of state can
still be formally recast as polytropic (12):

$p=\tilde{K}T^{1+n}\ ,$ (37)

where the polytropic constant $\tilde{K}$ is related to the degenerate one by

$\tilde{K}=\left(\dfrac{N_{A}k_{B}}{\mu}\right)^{-(n+1)}K^{-n}\ .$ (38)

Note that in the relation (37) we have used the ideal gas law given by

$\rho=\frac{\mu p}{N_{A}k_{B}T}\ .$ (39)

It is worth stressing that $K$ depends on the theory of gravity, since it can
be expressed with respect to the solution of the modified Lane-Emden equation
(15) via

$K=\left[\frac{4\pi}{\xi_{R}^{n+1}(-\theta^{\prime}_{n}(\xi_{R}))^{n-1}}\right]^{\frac{1}{n}}\frac{G}{n+1}M^{1-\frac{1}{n}}R^{\frac{3}{n}-1}\
.$ (40)

Note also that the equation of state (37) is valid up to the photosphere,
since above the interior-photosphere boundary the energy transport is ruled by
radiative processes instead. In such a region, we shall use a simplified
relation for the absorption law, given by the Kramer formula:

$\kappa_{\text{abs}}=\kappa_{0}p^{i}T^{j}\ .$ (41)

For cold stars, whose effective temperatures lie inside the range
$3000\lesssim T\lesssim 6000$K, the surface layer is dominated by H- opacity
hansen . Considering the hydrogen mass fraction as $X\approx 0.7$, the opacity
is given by

$\kappa_{H^{-}}=\kappa_{0}\rho^{\frac{1}{2}}\,T^{9}\,\,\text{cm}^{2}\text{g}^{-1}\
,$ (42)

with $\kappa_{0}\approx 2.5\times 10^{-31}\left(\frac{Z}{0.02}\right)$, where
the metal mass fraction $Z$ (or metallicity) is an important element in the
stellar modelling. Its value is typically taken within the range
$0.001\lesssim Z\lesssim 0.03$ cantiello ; as an example, the solar
metallicity is $Z=0.02$. Since that region can be also modelled as an ideal
gas, the opacity (42) can be expressed as

$\kappa_{H^{-}}=\kappa_{g}p^{\frac{1}{2}}\,T^{8.5}\,\,\text{cm}^{2}\text{g}^{-1}\
,$ (43)

where we have redefined
$\kappa_{g}=\kappa_{0}\left(\tfrac{\mu}{N_{A}k_{B}}\right)^{\frac{1}{2}}\approx
1.371\times 10^{-33}Z\mu^{\frac{1}{2}}$. Particularizing the relation (25) to
${H^{-}}$ opacity, one finds that the photospheric pressure is

$p_{ph}=\dfrac{2g\left(1-\dfrac{\epsilon}{R^{2}}\right)}{3\kappa_{H^{-}}}\ .$
(44)

Applying the solution of the modified Lane-Emden (15) for $n=3/2$ to the above
expression, the Stefan-Boltzmann law (26) and the opacity expression (43), the
photospheric pressure above takes the form

$p_{ph}=8.11279\times
10^{14}\left[\frac{M\,\beta}{\sqrt{\mu}LT^{4.5}Z}\right]^{2/3},$ (45)

where we have redefined the brackets appearing in Eq.(44) as

$\beta=1-\frac{2\alpha}{3^{5/3}\delta^{2/3}\omega^{2/3}}.$ (46)

For simplicity, we have removed the sub-indices $n$ appearing in the relations
(21) and from now on we will understand them as their values for $n=3/2$.

The above photospheric pressure must be matched to the pressure of the ideal
gas given by the relation (37) evaluated at the photosphere. The latter
yields, after using the Stefan-Boltzmann law (26), the effective temperature
under the form

$T_{eff}=9.1960\times
10^{-6}\left(\dfrac{\mu^{5}L^{3/2}Mp_{ph}^{2}}{-\theta^{\prime}\xi_{r}^{5}}\right)^{1/11}\
,$ (47)

which after using the derived photospheric pressure (45) provides the result

$T_{ph}=2482.10\mu^{\frac{13}{51}}\left(\frac{L}{L_{\odot}}\right)^{\frac{1}{102}}\left(\frac{M}{M_{\odot}}\right)^{\frac{7}{51}}\frac{\beta^{\frac{4}{51}}}{Z^{\frac{4}{51}}(\sqrt{-\theta^{\prime}}\xi_{R}^{5})^{\frac{1}{17}}}\
,$ (48)

where we have re-scaled both mass and luminosity to their solar values,
$\\{M_{\odot},L_{\odot}\\}$. Let us notice that the numerical value in the
above expression is too low; it should be almost a twofold larger. The reason
of this reduced value lies in the simplifications we have made, mainly related
to the photospheric modelling. Notwithstanding, this analytical formula allows
us to track down the modifications introduced by EiBI gravity to the early
stage of the stellar evolution. Therefore, for a given star with mass $M$,
uniform mean molecular weight $\mu$, and metallicity $Z$, the above expression
gives the corresponding Hayashi track. These almost vertical lines are
evolutionary tracks of infant stars with masses supposedly below $\sim
0.5M_{\odot}$, though such a limiting mass also depends on the theory of
gravity. Its shape and position on the HR diagram depend not only on the
metallicity, but also on the theory of gravity, which in the present case is
encapsulated in the parameter $\beta$ appearing in Eq.(46) but also through
the solutions of the extended Lane-Emden equation (15). This gravity model-
dependence could actually allow us to constraint the theory’s parameter when a
more realistic atmosphere analysis is performed, since none such star is
placed in the so-called Hayashi forbidden zone, being found in the low
temperatures. As presented in Fig. 1, EiBI gravity shifts the tracks in
opposite ways depending on the sign of the theory’s parameter $\epsilon$,
either in the direction of the forbidden zone (for $\epsilon>0$), which lies
in the region of lower temperatures, or against it (for $\epsilon<0$).

Figure 1: The piece of the HR diagram representing shifted Hayashi tracks by
the modifications introduced by EiBI gravity model in logarithmic scale. The
curves are given by the equation (48) taking $M/M_{\odot}=1/2$, for some
chosen values of the parameter $\alpha$ defined in Eq.(16) as compared to the
GR/Newtonian curve, $\alpha=0$.

###  III.2 Minimum main sequence mass

Using the ingredients introduced in the previous section, we are now capable
to compute the minimum main sequence mass (MMSM). This is the minimal mass
required by a star to ignite sufficiently stable thermonuclear reactions in
its interior to compensate photospheric energy losses. Even though the central
temperature can be sufficient to start the p-p chain, it is not necessarily
enough to complete it. The thermonuclear rates depend mainly on the
temperature and density, in such a way that the energy generation rate can be
well approximated by power laws of the form (see Burrows:1992fg  for details)

$\dot{\epsilon}_{pp}=\dot{\epsilon}_{c}\left(\frac{T}{T_{c}}\right)^{s}\left(\frac{\rho}{\rho_{c}}\right)^{u-1}\
,$ (49)

where the two exponents can be phenomenologically fitted as $s\approx 6.31$
and $u\approx 2.28$ at the transition mass of the core, while the function

$\dot{\epsilon}_{c}=\dot{\epsilon}_{0}\,T_{c}^{s}\,\rho_{c}^{u-1}\,\text{ergs
g}^{-1}\text{s}^{-1}\ ,$ (50)

with $\dot{\epsilon}_{0}\approx 3.4\times 10^{-9}$ in suitable units. The
corresponding luminosity of the hydrogen burning is found as

$L_{pp}=\int\dot{\epsilon}_{pp}\,dM=4\pi\dot{\epsilon}_{c}r_{c}^{3}\rho_{c}\int_{0}^{\xi_{R}}\theta^{n\left(u+\frac{2}{3}s\right)}\,\xi^{2}d\xi\
,$ (51)

where we have used the fact that $(T/T_{c})=(\rho/\rho_{c})^{2/3}$ along the
adiabatic core. The last integral can be easily computed by using the
approximation (24) and the definition (17), which yields the result

$L_{pp}=\dfrac{6\sqrt{3\pi(\alpha+1)^{3}}}{\omega_{3/2}(2s+3u)^{3/2}}\dot{\epsilon}_{c}M\
.$ (52)

In this approximation, we have taken into account the fact that most of the
hydrogen burning will be produced in a region near to the core of the star.
Now, considering that the fraction of hydrogen in a high-mass brown dwarf is
of $75\%$, and that the number of barions per electron can be approximated to
$\mu_{e}=1.143$, besides setting the following degenerate polytropic constant

$K=\dfrac{(3\pi^{2})^{2/3}\hbar}{5m_{e}m_{H}^{5/3}\mu_{e}^{5/3}}\left(1+\frac{\alpha_{d}}{\eta}\right),$
(53)

where $\hbar$ is the reduced Plank constant, $m_{e}$ is the electron mass and
$m_{H}$ is the proton mass. Then, the luminosity (52) can be recast as

$L_{pp}=1.54\times
10^{7}L_{\odot}\dfrac{\delta^{5.49}\,(1+\alpha)^{3/2}}{\gamma^{16.46}\omega}M_{-1}^{11.97}\,\dfrac{\eta^{10.15}}{\left(\alpha_{d}+\eta\right)^{16.46}}\
,$ (54)

where we have defined here, by convenience, $M_{-1}=M/(0.1M_{\odot})$. In
order to find the MMSM we need to equal this hydrogen burning luminosity to
the one of the photosphere. For the purpose of computing the latter, we take
Eq.(44) and assume again that the components of the stellar atmosphere behave
as ideal gas, that is

$\dfrac{\rho_{ph}k_{B}T_{ph}}{\mu
m_{H}}=\frac{2g\,\left(1-\dfrac{\epsilon}{R^{2}}\right)}{3\kappa_{op}}\ .$
(55)

This equation will allow us to get a relation between $\rho_{ph}$ and $M$, but
before going that way, let us first rewrite the surface gravity $g$ defined in
Eq.(28) as

$g=\frac{G^{3}M^{5/3}}{\gamma^{2}K^{2}}\ ,$ (56)

and consider that the photospheric temperature can be found from the matching
of the specific entropies of the gas/metallic phases there, which yields
Burrows:1992fg

$T_{ph}=\frac{1.8\times 10^{6}\rho_{ph}^{0.42}}{\eta^{1.545}}\ .$ (57)

Replacing the above two equations into (55), we find

$\rho_{ph}=2.957\times 10^{-5}\dfrac{\eta^{1.09}G^{2.11}M^{1.17}\left(\mu
m_{H}\beta\right)^{0.70}}{\left(\gamma
K\right)^{1.41}(k_{B}\kappa_{op})^{0.70}}\ .$ (58)

Inserting this result back to Eq.(57), the photospheric temperature becomes

$T_{ph}=2.254\times 10^{4}\frac{G^{0.89}M^{0.49}\left(\mu
m_{H}\beta\right)^{0.30}}{\eta^{1.09}\left(\gamma
K\right)^{0.59}(k_{B}\kappa_{R})^{0.30}}\ .$ (59)

Therefore, the photospheric luminosity, given by $L_{ph}=4\pi R^{2}\sigma
T_{ph}^{4}$, can be expressed in terms of the star mass, $M$, as follows

$L_{ph}=0.534L_{\odot}\dfrac{M_{-1}^{1.31}\beta^{1.18}}{\eta^{3.99}\gamma^{0.37}(\alpha_{d}+\eta)^{0.37}\kappa_{-2}^{1.18}}\
,$ (60)

where we have defined the quantity
$\kappa_{-2}=\kappa_{R}/(10^{-2}\text{cm}^{2}\text{g}^{-1})$. Finally,
equalling the hydrogen burning luminosity (54) with the photospheric one (60),
we find the following expression

$M_{-1}^{\text{MMSM}}=0.227\dfrac{\gamma^{1.51}\omega^{0.09}(\alpha_{d}+\eta)^{1.51}\beta^{0.11}}{(\alpha+1)^{0.14}\delta^{0.51}\eta^{1.33}\kappa_{-2}^{0.11}},$
(61)

where the EiBI dependences enter in this expression both via the coefficient
$\alpha$ in Eqs.(16) and Eq.(46) and via the parameters
$\\{\omega,\gamma,\delta\\}$ obtained from the resolution of the extended
Lane-Emden equation (15). This is the MMSM for EiBI gravity under the
assumptions and simplifications above. In order to compute it for different
values of the EiBI parameter, the main obstacle here is the fact that $\alpha$
depends on the central density of the star, as can be seen from Eqs.(14) and
(16). Therefore, for the sake of our calculations we shall take the maximal
value for the central density of $\rho_{c}\sim 10^{3}\text{g/cm}^{3}$
Burrows:1992fg , which allow us to compute the MMSM and thus set bounds on the
size of $\alpha$ as coming from observational constraints. In Table 1 we
actually compute the set of $\\{\gamma,\omega,\delta\\}$ values for several
choices of the parameter $\alpha$ in order to find the corresponding MMSM. For
$\alpha=0$ (GR case) we get $M^{\text{MMSM}}\approx 0.084M_{\odot}$, which is
somewhat halfway between other analytical calculations Burrows:1992fg  and the
results of numerical simulations Kumar . For non-vanishing values of $\alpha$,
this table shows that for positive (negative) $\alpha$ the MMSM is larger
(smaller). Thus, the positive branch of $\alpha$ is the most interesting one
for our purposes, since it allow us to constrain its size via comparison with
the observations of the less massive main-sequence stars ever observed, which
corresponds to the $0.0930\pm 0.0008M_{\odot}$ of the M-dwarf star G1 866C
Segransan . This way, in Fig. 2 we numerically depict the evolution of the
MMSM with $\alpha>0$. Our results within our simplified model points that near
values of the parameter $\alpha\gtrsim 0.1$ the model is likely to run into
conflict with observations, therefore setting a bound to the combination of
the EiBI parameter and the star’s central density, the latter to be estimated
by other means. On the other hand, in the negative branch we run into a
problem related to the fact that when the parameter reaches
$\alpha\lesssim-0.003$ there are non-physical solutions, since below that
value the sign would change in the bracket of the extended Lane-Emden equation
(15). Let us however note that, due to the same reasons stated above, this
feature does depend on the star’s central density, hence for a less dense or a
denser core we would deal with different singular values of the parameter and,
therefore, we do not extract any conclusion on the limit of validity of this
branch within the formalism presented here.

$\alpha$ | $\gamma_{3/2}$ | $\omega_{3/2}$ | $\delta_{3/2}$ | $\text{MMSM}/M_{\odot}$
---|---|---|---|---
0.100 | 2.49 | 3.01 | 5.74 | 0.0933
0.010 | 2.37 | 2.74 | 5.96 | 0.0852
0.001 | 2.36 | 2.72 | 5.99 | 0.0845
0 (GR) | 2.36 | 2.71 | 5.99 | 0.0844
-0.001 | 2.36 | 2.71 | 5.99 | 0.0843
-0.003 | 2.35 | 2.69 | 6.03 | 0.0837
Table 1: The MMSM (in units of solar masses) computed with Eq.(61) for several
values of the parameter $\alpha$ defined in (16), including the intermediate
values of the parameters $\\{\gamma,\omega,\delta\\}$ obtained from the
resolution of the extended Lane-Emden equation (15). Figure 2: The evolution
of the MMSM (in units of solar masses) with the parameter $\alpha$. The range
in this plot goes from $\alpha\in(-0.003,0.10)$, with the lower bound given by
a well-defined solution of the extended Lane-Emden equation (15).

###  III.3 Fully convective stars on the main sequence and radiative core
development

We will now focus on the final parts of the Hayashi tracks. Recalling that
during this evolutionary phase the proto-star is fully convective, it might
happen that the inner temperature increases enough to satisfy the conditions
for radiative core development and, therefore, the object can have much more
complex structure than the one we consider. Because of that, as discussed in
Sec.III, the star can either enter the Henyey evolutionary phase, represented
by the almost horizontal to the main sequence lines, or it can stop
contracting on the onset of the radiative core development as the nuclear
processes start balancing the gravitational attraction. In the last situation,
which is our concern now, the star begins its main sequence evolutionary
phase, being however still fully convective. Using the previous results found
in this paper, we can obtain the Maximum Fully Convective Mass (MFCM) of a
star on the main sequence.

In order to determine when the radiative processes take over in the core, one
needs to analyze the Schwarzschild criterion, which tells us that the
radiative processes start when $\nabla_{ad}=\nabla_{rad}$. In our simplified
modelling above we assumed that the star is made of an ideal, monoatomic gas,
providing that $\nabla_{ad}=0.4$, while the radiative temperature gradient was
already derived in Eq.(36). Applying the homology law, together with Eqs.(37)
and (38), we can express the latter as

$\nabla_{rad}=5.21177\times
10^{69}\frac{L\xi^{5}(-\theta^{\prime})\kappa_{o}}{\mu^{5}M^{2}R^{3}T^{3.5}\beta}\
,$ (62)

where $L$ is the local luminosity (here evaluated at the core). Substituting
the central temperature from (20) and subsequently the Stefan-Boltzmann law,
the above expression yields

$\nabla_{rad}=8.99\times
10^{-13}\left(\frac{L}{L_{\odot}}\right)^{1.25}\frac{\xi^{10.83}\left(-\theta^{\prime}\right)^{2.17}\kappa_{0}}{\delta^{2.33}\mu^{8.5}M_{-1}^{5.5}T}\
.$ (63)

Equaling this result with $\nabla_{ad}=0.4$ one finds the maximum luminosity
for a fully convective star on the onset of the radiative core development

$L=2.0827\times
10^{9}L_{\odot}\dfrac{\beta^{0.8}\delta^{1.87}\mu^{6.8}T^{0.8}}{\xi^{8.67}\left(-\theta^{\prime}\right)^{1.73}\kappa_{0}^{0.8}}M_{-1}^{4.4}\
.$ (64)

Now, by equaling this luminosity to the one of the hydrogen burning given by
Eq.(60) yields the MFCM

$M_{-1}=1.91\dfrac{\beta^{0.11}\gamma^{2.17}\mu^{0.90}T^{0.11}\omega^{0.13}(\alpha_{d}+\eta)^{2.17}}{\left(\alpha+1\right)^{0.20}\delta^{0.48}\eta^{1.34}\xi^{1.14}\left(-\theta^{\prime}\right)^{0.23}\kappa_{0}^{0.11}}\
,$ (65)

where similar comments as on the sources of EiBI corrections of the MMSM above
apply here. To make quantitative estimates of this mass, let us first assume
the usual values for LMS as $\alpha_{d}=4.82$, $\eta=9.4$, $\mu=0.618$ and
$T_{\text{eff}}=4000$K. In addition, we have to set the opacity, keeping the
Kramers’ form written in Eq.(41) with $i=1$ and $j=-4.5$; there are the total
bound-free and the free-free opacities (see e.g. hansen  for details)

$\displaystyle\kappa_{0}^{bf}$ $\displaystyle\approx$ $\displaystyle 4\times
10^{25}\mu\dfrac{Z(1+X)}{N_{A}k_{B}}\,\text{cm}^{2}\text{g}^{-1},$ (66)
$\displaystyle\kappa_{0}^{ff}$ $\displaystyle\approx$ $\displaystyle 4\times
10^{22}\mu\dfrac{(X+Y)(1+X)}{N_{A}k_{B}}\,\text{cm}^{2}\text{g}^{-1}.$ (67)

$\alpha$ | $M_{bf}/M_{\odot}$ | $M_{ff}/M_{\odot}$
---|---|---
0.100 | 0.108 | 0.225
0.010 | 0.0994 | 0.207
0.001 | 0.0985 | 0.205
0.000 | 0.0984 | 0.205
-0.001 | 0.0983 | 0.204
-0.003 | 0.0976 | 0.203
Table 2: Numerical values for the MFCM (in solar mass units), using the total
bound-free and the free-free opacities, defined in Eq. (66) and (67),
respectively, for different values of the composite EiBI parameter $\alpha$
appearing in Eq.(16). Figure 3: The dependence of the (normalized) MFCM, for
both opacity models (66) and (67), on the parameter $\alpha\in(-0.003,0.10)$.

Once everything is settled, in Table 2 we calculate the MFCM for several
values of the parameter $\alpha$ for both opacity models. Similarly as with
the MMSM above, the MFCM increases (decreases) with positive (negative)
gravitational parameter (note that the parameters $\\{\omega,\gamma,\delta\\}$
are those appearing in Table 1). In addition, in Fig. 3 we depict the
evolution with $\alpha$ of the two MFCM masses, corresponding to each opacity.
In both table and plot it is clearly seen that the choice of the opacity model
significantly affects (roughly a factor two) the value of the MFCM. As for the
negative branch, we find the same feature as with the MMSM, namely, the fact
that for $\alpha\lesssim-0.003$ the extended Lane-Emden equation (15) fails to
provide a non-singular solution and, as such, those values are disregarded in
our analysis of the MFCM.

##  IV Discussion and conclusion

In this work we have discussed several aspects of the early evolutionary
phases of low-mass stars within an extension of GR dubbed as Eddington-
inspired Born-Infeld gravity. Such an extension is governed by a single
parameter which manifests, in the stellar structure equations of non-
relativistic stars, as an extra piece to the Poisson (Lane-Emden) equation
when a polytropic equation of state is considered. Supplemented with a
simplified photospheric model, and equipped with a criterion for conventive
instability, we investigated three features of such an early evolution of LMS.

The first feature deals with the effective temperature-luminosity relations in
the evolutionary path of a proto-star, the so-called Hayashi tracks. We have
shown that positive (negative) values of the EiBI parameter shift the
corresponding Hayashi track in the sense of larger (smaller) effective
temperature for a fixed luminosity. The second feature is the minimum required
mass for a star to stably burn enough hydrogen to compensate photospheric
losses, allowing it to belong to the main sequence. In this case, positive
(negative) values of the EiBI parameter yield larger (smaller) minimum main-
sequence masses, the former allowing to place constraints on the theory’s
parameter via comparison with the lowest-mass main-sequence stars every
observed. This poses a difficulty for this theory, since such constraints act
upon a combination of the EiBI parameter and the star’s central density. This
dependence of the stellar features not only on global quantities (such as the
total mass) but also on local ones is a common feature of the RBG family,
therefore forcing us to live with it. The third feature deals with the
development of a radiative core at the end of the Hayashi track, entering the
main-sequence phase while still being fully convective. We found the maximum
value of the mass for this to happen, again observing an increase (decrease)
of this mass with positive (negative) EiBI parameter. Note, however, that for
all these three features the main astrophysical ingredient determining their
absolute values is the opacity, whose modelling is always a delicate issue.
Its influence is obvious in the last feature (the MFCM), where two different
models of opacities (bound-free and free-free ones) result in up to a factor
two in the absolute value of this quantity.

The results found in this work highlight the viability of using metric-affine
gravities of the RBG type to study modifications to the stellar model
predictions of GR, in particular, within the non-relativistic regime. This is
so because in such a regime, RBG modifications to the usual Poisson equation
typically occur via a single additional parameter Olmo:2021yul , allowing to
study the phenomenology of several types of stars, particularly low-mass
stars, without ruining the consistence of the theory with weak-field limit
observations. As mentioned above, the main bottleneck in order to place
observational constraints upon any such theories is the determination of the
central density, which up to now we have been only able to fix by taking its
assumed values within GR, though more reliable theoretical procedures to deal
with this issue are being investigated. We are also working in other aspects
of non-relativistic objects and low-mass stars in different RBGs, and we hope
to report on this soon.

## Acknowledgements

MG is funded by the predoctoral contract 2018-T1/TIC-10431 and acknowledges
further support by the European Regional Development Fund under the Dora Plus
scholarship grants. DRG is funded by the Atracción de Talento Investigador
programme of the Comunidad de Madrid (Spain) No. 2018-T1/TIC-10431, and
acknowledges further support from the Ministerio de Ciencia, Innovación y
Universidades (Spain) project No.
PID2019-108485GB-I00/AEI/10.13039/501100011033, and the FCT projects No.
PTDC/FIS-PAR/31938/2017 and PTDC/FIS-OUT/29048/2017. AW is supported by the EU
through the European Regional Development Fund CoE program TK133 “The Dark
Side of the Universe”. This work is also supported by the Spanish projects
Nos. FIS2017-84440-C2-1-P and PID2020-117301GA-I00, the project
PROMETEO/2020/079 (Generalitat Valenciana), and the Edital 006/2018 PRONEX
(FAPESQPB/CNPQ, Brazil, Grant 0015/2019). MG thanks the Laboratory of
Theoretical Physics at the University of Tartu, and AW the Departament of
Theoretical Physics at the Complutense University of Madrid, for their kind
hospitality while doing this work.

## References

  * (1) C. M. Will, Living Rev. Rel. 17 (2014) 4.
  * (2) S. Capozziello and M. De Laurentis, Phys. Rept. 509 (2011) 167.
  * (3) S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692 (2017) 1.
  * (4) P. Bull, et al. Phys. Dark Univ. 12 (2016) 56.
  * (5) J. M. M. Senovilla and D. Garfinkle, Class. Quant. Grav. 32 (2015) 124008.
  * (6) J. M. Ezquiaga and M. Zumalacárregui, Phys. Rev. D 102 (2020) 124048.
  * (7) D. Psaltis, et al. [Event Horizon Telescope], Phys. Rev. Lett. 125 (2020) 141104.
  * (8) G. J. Olmo, D. Rubiera-Garcia and A. Wojnar, Phys. Rept. 876 (2020) 1.
  * (9) B. P. Abbott, et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 119 (2017) 161101.
  * (10) M. Camenzind, “Compact Objects in Astrophysics” (Springer, 2007).
  * (11) M. Linares, T. Shahbaz and J. Casares, Astrophys. J. 859 (2018) 54.
  * (12) J. Antoniadis, et al. Science 340 (2013) 6131.
  * (13) R. Abbott, et al. [LIGO Scientific and Virgo], Astrophys. J. Lett. 896 (2020) L44.
  * (14) J. Sakstein, Int. J. Mod. Phys. D 27 (2018) 1848008.
  * (15) S. Chowdhury and T. Sarkar, JCAP 05 (2021) 040.
  * (16) M. Benito and A. Wojnar, Phys. Rev. D 103 (2021) 064032.
  * (17) A. Wojnar, Phys. Rev. D 104 (2021) 104058.
  * (18) A. Wojnar, Phys. Rev. D 103 (2021) 044037.
  * (19) A. Wojnar, Int. J. Geom. Meth. Mod. Phys. 18 (2021) 2140006.
  * (20) L. Sarmah, S. Kalita and A. Wojnar, [arXiv:2111.08029 [gr-qc]].
  * (21) J. Sakstein, Phys. Rev. Lett. 115 (2015) 201101.
  * (22) J. Sakstein, Phys. Rev. D 92 (2015) 124045.
  * (23) M. Crisostomi, M. Lewandowski and F. Vernizzi, Phys. Rev. D 100 (2019) 024025.
  * (24) G. J. Olmo, D. Rubiera-Garcia and A. Wojnar, Phys. Rev. D 100 (2019) 044020.
  * (25) A. S. Rosyadi, A. Sulaksono, H. A. Kassim and N. Yusof, Eur. Phys. J. C 79 (2019) 1030.
  * (26) A. Wojnar, Phys. Rev. D 102 (2020) 124045.
  * (27) A. Kozak and A. Wojnar, Phys. Rev. D 104 (2021) 084097.
  * (28) A. Kozak and A. Wojnar, [arXiv:2110.15139 [gr-qc]].
  * (29) A. Kozak and A. Wojnar, [arXiv:2111.11199 [gr-qc]].
  * (30) M. Banados and P. G. Ferreira, Phys. Rev. Lett. 105 (2010) 011101.
  * (31) P. Pani, V. Cardoso and T. Delsate, Phys. Rev. Lett. 107 (2011) 031101.
  * (32) J. Casanellas, P. Pani, I. Lopes and V. Cardoso, Astrophys. J. 745 (2012) 15.
  * (33) T. Delsate and J. Steinhoff, Phys. Rev. Lett. 109 (2012) 021101.
  * (34) P. P. Avelino, Phys. Rev. D 85 (2012) 104053.
  * (35) S. Banerjee, S. Shankar and T. P. Singh, JCAP 10 (2017) 004.
  * (36) M. D. Danarianto and A. Sulaksono, Phys. Rev. D 100 (2019) 064042.
  * (37) C. Wibisono and A. Sulaksono, Int. J. Mod. Phys. D 27 (2018) 1850051.
  * (38) P. Banerjee, D. Garain, S. Paul, R. Shaikh and T. Sarkar, [arXiv:2105.09172 [gr-qc]].
  * (39) J. Beltran Jimenez, L. Heisenberg, G. J. Olmo and D. Rubiera-Garcia, Phys. Rept. 727 (2018) 1.
  * (40) C. Hayashi, Publication of the Astronomical Society of Japan, Vol.13, 450-452 (1961).
  * (41) V. I. Afonso, C. Bejarano, J. Beltran Jimenez, G. J. Olmo and E. Orazi, Class. Quant. Grav. 34 (2017) 235003.
  * (42) V. I. Afonso, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 97 (2018) 021503.
  * (43) G. J. Olmo, D. Rubiera-Garcia and A. Wojnar, Phys. Rev. D 104 (2021) 024045.
  * (44) A. Wojnar, Eur. Phys. J. C 79 (2019) 51.
  * (45) C. J. Hansen, S. D. Kawaler, V. Trimble, Stellar Interiors. Physical Principles, Structure, and Evolution, 2nd edition (Springer 2004).
  * (46) L. Henyey, R. Lelevier, and R.D. Levee, Publications of the Astronomical Society of the Pacific 67.396 (1955): 154\.
  * (47) L. Henyey, J.E. Forbes, N.L. Gould, The Astrophysical Journal 139 (1964) 306.
  * (48) L. Henyey, M.S. Vardya, P. Bodenheimer, The Astrophysical Journal 142 (1965) 841.
  * (49) A. Burrows and J. Liebert, Rev. Mod. Phys. 65 (1993) 301.
  * (50) M. Cantiello, J. Braithwaite, The Astrophysical Journal 883 (2019) 106.
  * (51) S. Kumar, Astrophys. J. 137 (1963) 1121.
  * (52) D. Segransan, et al. Astron. Astrophys. 364 (2000) 665.