1 Introduction

Recently, the LIGO/Virgo scientific collaboration announced the discovery of a compact binary coalescence, GW190814, involving a 22.2–24.3\(M_\odot \) black hole (BH) and a compact object with a mass of 2.50–2.67\(M_\odot \) [1]. The secondary component is either the most massive neutron star (NS), or the lightest BH ever discovered in a double compact-object system. The existence of ultra-massive NSs has been revealed in several studies. So far, the precise determination of the mass of pulsars leads to \(M=1.97\pm 0.04\) \(M_\odot \) for PSR J1614–2230 [2], PSR J0348+0432 has \(M=2.01\pm 0.04\) \(M\odot \) [3] and the highly likely the most massive NS observed to date, PSR J0740 + 6620, with \(M=2.14\pm 0.1\) \(M_\odot \) [4].

The discovery of the gravitational wave (GW) binary GW190814 triggered intense theoretical efforts about the real nature of its secondary component, in particular, the high mass requires the compact star matter to be described by a stiff equation of state (EoS). In [5] is argued that fast rotation is capable to explain the existence of a stable \(\sim 2.6 M_\odot \) NS for moderately stiff EoS but may not be adequate for soft EoSs. In particular, several soft EoSs favored by GW170817 and with maximum spherical masses of \(\sim 2.1 M_\odot \) cannot be used to explain this source as a uniformly spinning NS. In [6] the authors infer a lower limit on the maximum mass \(M_{\mathrm{max}}\) of non-rotating NSs, using arguments based on universal relations connecting the masses and spins of uniformly rotating NSs. Furthermore, they obtain a lower bound on the dimensionless spin for the secondary companion, using the upper maximum mass constraints from the GW170817 event [7, 8], and show that the allowed range in dimensionless spins correspond to rotational frequencies much higher than the fastest millisecond pulsars known [9]. In [10] a rigorous upper bounds on maximum mass under the exclusive assumptions of causality and general relativity (GR), showing that the presence of a NS in GW190814 is not inconsistent with present observational constraints on the NS EoS. On the other hand, in [11] is showed that the stiffening of the EoS required to support ultra-massive NSs is inconsistent with either constraints obtained from the low deformability of medium-mass stars demanded by GW170817 or from energetic heavy-ion collisions. Several other methodology and speculations about GW190814 are presented in [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28].

As is well known, the two main observable of a NS, i.e., their mass and radius, both depend crucially on the choice of EoS and the gravitational theory, where gravity will determine the macroscopic hydrostatic equilibrium. The general relativity is very well tested and compatible in weak-field observations, for instance, from solar system and terrestrial experiment tests [29,30,31], but there are theoretical and observational reasons to believe that GR should be modified when gravitational fields are strong and/or on large scales. On large scales and from an observational point of view, the physical mechanism responsible for accelerating the Universe at late times is still an open question, and new degrees of freedom of the gravitational origin are alternatives to explain such an accelerated stage (see, e.g., [32, 33] for review). Theories beyond GR can serve as alternatives to explain the current tension in the Hubble constant that persists in the framework of the standard cosmological model [34,35,36]. Also, modified gravity models are motivated to drive the accelerating expansion of the Universe at early times, known as inflationary era. On astrophysical scales, compact objects such as BHs and NSs are our best natural laboratories to constrain strong gravity. In these bodies, gravity prevails over all other interactions and collapse leads to large-curvature and strong-gravity environments [37]. We refer the reader to [38] and references therein for several modified gravity scenarios motivation under the regime of strong gravitational field.

A consistent way to modify gravity is to introduce some screening mechanism in order to have viable gravity, i.e., GR, at small distances and scales and possible relevant effects on others environmental scales, for instance, on compact objects or even on cosmological level. Recall that the screening mechanisms include chameleons [39, 40], symmetrons [41], dilatons [42], Vainshtein mechanism [43], etc. At the heart of screening mechanisms lies the fact that there are 29 orders-of-magnitude separating the cosmological and terrestrial densities and 20 orders of magnitude separating their distance scales. As a result, the properties of the new degrees of freedom of the gravitational origin (a scalar field) can vary wildly in different environments [44]. Screening gravity have been intensively investigated in the literature, in the most diverse aspects in cosmology and gravitation (see, e.g., [45,46,47,48,49,50,51,52,53,54,55,56,57,58] for a short list). A compilation/list of various works on compact objects in modified theories of gravity can be found in [59].

As regards to NS observations, it seems there is a theoretical degeneration and is not clear if measurement of mass and radius constrain possible gravity effects or EoSs. As showed in [60], measurement of mass constrains gravity rather than the EoS (see also [61]). Thus, it is difficult to distinguish which of these effects is actually contributing to the actual observed mass and radius values. In this work, we analyze the NS mass-radius relation through a modification of the TOV equations induced by the presence of a possible chameleon screening (thin-shell effect), where the microscopic physics inside the NS will be modeled by realistic soft EoSs of which it is not possible to generate very massive NS with \(\sim 2.6\) \(M_\odot \) or even \(\sim 2.14\) \(M_\odot \) in GR context. Our main aim is to show that an appropriate combination of modified gravity, rotation effects and realistic soft EoSs, can have a joint effect for achieve high masses and explain GW190814 secondary component, and in return also NS like PSR J0740 + 6620.

This paper is organized as follows. In the next section, we present our theoretical framework and methodology to generate the NS mass-radius relation. In Sect. 3 we present our main results and in Sect. 4 our final comments and perspectives. (The speed of light c is set equal to unity).

2 Theoretical framework and methodology

In this section we review in a nutshell the theoretical methodology used to obtain our main results.

2.1 Screening mechanisms

The study of scalar–tensor theories has been motivated by some cosmological observation, especially in order to explain the current accelerating expansion stage of the Universe. One proposed explanation is that gravity is modified on large scales, but must be suppressed on small scales, for instance, in the solar system deviations are constrained to be subdominant by a factor of \(10^{-5}\). Screening gravity circumvent this problem by introducing non-linear modifications of the Poisson equation that dynamically suppress deviations from GR in the solar system without the need to fine-tune on the scalar mass or the coupling to matter. In short, screening mechanisms utilize non-linear dynamics to effectively decouple solar system and cosmological scales. At the heart of screening mechanisms lies the fact that there are 29 orders-of-magnitude separating the cosmological and terrestrial densities. As a result, the properties of the scalar field (new degree of freedom of gravitational origin) can vary widely in different environments.

In this work, we will consider the well-studied chameleon screening [39, 40]. In chameleon models, the mass of the scalar is an increasing function of the ambient density. The purpose of this section is review the calculation of the parameterized post-Newtonian (PPN) parameter \(\gamma \) that is relevant for chameleon theories. We follow the methodology presented in [62] and references therein.

Let us consider a specific subset of the general scalar–tensor theories, which in the Einstein frame is given by

$$\begin{aligned} S= & {} \int dx^4 \sqrt{-g}\left[ \frac{M_{\mathrm{pl}}^2}{2}R-\frac{1}{2}\nabla _\mu \phi \nabla ^\mu \phi -V(\phi )\right] \nonumber \\&+S_m[{\tilde{g}}], \end{aligned}$$
(1)

where the Jordan frame metric \({\tilde{g}}_{\mu \nu }\) is a Weyl rescaling of the Einstein frame metric \(g_{\mu \nu }\) by a conformal factor \(A(\phi )\), i.e., \({\tilde{g}}_{\mu \nu }=A^2(\phi )g_{\mu \nu }\). The coupling is given by

$$\begin{aligned} \alpha (\phi )=M_{\mathrm{pl}}\frac{d \ln A(\phi )}{d \phi }. \end{aligned}$$
(2)

Each particular scenario is set by the choice of \(A(\phi )\) and \(V(\phi )\) functions. The PPN metric for a single body, which we will refer to as body A with mass \(M_A\), reads

$$\begin{aligned} {\tilde{g}}_{00}&=-1 +2 \frac{G^{PPN} M_A}{r} \end{aligned}$$
(3)
$$\begin{aligned} {\tilde{g}}_{ij}&=\left( 1+2{\tilde{\gamma }} \frac{G^{PPN} M_A}{r}\right) \delta _{ij}. \end{aligned}$$
(4)

We will refer to the quantity \(G^{PPN}\) as the PPN gravitational constant. \(G^{PPN}\) controls the size of effects computed using the PPN metric. It is distinct from the gravitational constant G that appears in GR and Newtonian gravity. We are interested in the regime where some body has some degree of screening. In this case, the equation inside of some screening radius is

$$\begin{aligned} \nabla ^2\phi =0 \end{aligned}$$
(5)

while outside the screening radius reads

$$\begin{aligned} \nabla ^2\phi =8\pi \alpha (\phi ) G\rho _A\quad r > r_s, \end{aligned}$$
(6)

where \(\alpha \) is the coupling function.

In order to move on, we need to define \(\alpha \). Here, we assume the chameleon field which uses a non-linear potential to make the field mass a function of the environmental density. The equation of motion reads

$$\begin{aligned} \nabla ^2\phi =-n\frac{\Lambda ^{4+n}}{\phi ^{n+1}}+\frac{\alpha \rho _A}{M_{\mathrm{pl}}}, \end{aligned}$$
(7)

where the effective potential is given by

$$\begin{aligned} V_{\mathrm{eff}}=\frac{\Lambda ^{4+n}}{\phi ^{n}}+\frac{\alpha \phi \rho _A}{M_{\mathrm{pl}}}. \end{aligned}$$
(8)

The mass-scale \(\Lambda \) can vary over many orders of magnitude, but it is often compared to the dark energy scale, since this value is relevant for the present-day cosmic acceleration. Astrophysically, the chameleon profile of a spherically-symmetric object of mass M and radius R is not sourced by the object’s mass, but rather by the mass inside a shell near the surface, a phenomenon that has been dubbed the thin-shell effect. The reason for this is the following: deep inside the object, the field minimizes its effective potential corresponding to the ambient density but, as one moves away from the center, the field must eventually evolve asymptotically to the minimum at the density of the medium in which the object is immersed (astrophysical or cosmological densities, depending on the situation to be investigated). The field can only roll once the density is low enough so that its effective mass is light enough. The radius at which this happens is typically called the screening radius \(r_{\mathrm{s}}\).

The coupling function \(\alpha \) above is a constant for chameleons [44, 62]. Ignoring possible scalar mass contribution, the solution is then

$$\begin{aligned} \phi =\phi _0-2\alpha ^2\frac{Q_A G M_A}{r}, \end{aligned}$$
(9)

where Q is the scalar charge of body A and is given by

$$\begin{aligned} Q_A=\left( 1-\frac{M_A(r_s^A)}{M_A}\right) . \end{aligned}$$
(10)

Transforming to the Jordan frame and expanding \(A(\phi )\) to first order in \(GM_A/r\), one finds Eqs. (3) and (4) with

$$\begin{aligned} G^{PPN}&=G\left[ 1+2\alpha ^2Q_A\right] \quad \mathrm{and} \end{aligned}$$
(11)
$$\begin{aligned} {\tilde{\gamma }}&= \frac{1-2\alpha ^2Q_A}{1+ 2\alpha ^2Q_A}. \end{aligned}$$
(12)

In a binary system, if we consider the orbital dynamics of a body of mass \(M_B\) orbiting the body sourcing this metric, this second body may have its own screening radius \(r_s^B\), so that the chameleon force between the two bodies A and B reads [63]

$$\begin{aligned} F=-\frac{G M_A M_B}{r^2} [1+2 \alpha ^2 Q_A Q_B ] \end{aligned}$$
(13)

with the following modified potential

$$\begin{aligned} V(r)=\frac{GM}{r}\left[ 1+2\alpha ^2\left( 1-\frac{M(r_{\mathrm{s}})}{M}\right) \right] , \end{aligned}$$
(14)

where we have once again ignored the mass of the scalar. Thus, the physical quantity that is measurable in these theories is

$$\begin{aligned} G^{PPN}=G [1+2 \alpha ^2 Q_A Q_B], \end{aligned}$$
(15)

with corresponding PPN parameter \(\gamma \) given by

$$\begin{aligned} \gamma = 2[1+2\alpha ^2Q_AQ_B]^{-1}-1. \end{aligned}$$
(16)

Notice that other approaches have been developed in the literature, see for instance [44, 54, 55]. In what follows, let us quantify how this approximation presented above can modify spherical compact stars.

2.2 Modified TOV

In GR, the hydrostatic equilibrium of a star is described by the Tolman–Oppenheimer Volkoff (TOV) equations (see, e.g., [64]), namely

$$\begin{aligned} \frac{dP}{dr}= & {} - \frac{G M(r) \rho (r)}{r^2}\left( 1 + \frac{P(r)}{\rho (r)} \right) \nonumber \\&\times \left( 1 + \frac{4 \pi r^3 P(r)}{M(r)} \right) \left( 1 - \frac{2GM(r)}{r} \right) ^{-1} \end{aligned}$$
(17)

and

$$\begin{aligned} \frac{dM}{dr} = 4 \pi r^2 \rho (r), \end{aligned}$$
(18)

where \(\rho (r)\) is the energy density, P(r) is the pressure and M(r) is the mass within the radial coordinate r.

The mass and the radius of a star are the two more obvious observables. Evidently, these quantities depends on the theory of gravitation and microscopic physics, i.e, the EoS of the matter of the star. In view of the formalism described in the previous section, we can notice that the chameleon screening can be quantified in the hydrostatic equilibrium rescaling G by \(G^{PPN}\) in Eqs. (17) and (18). \(G^{PPN}\) is also directly connected with the PPN parameter \(\gamma \), since \(\alpha ^2 Q_A Q_B\) can be estimated.

In this work, we will restrict to investigating NSs. Thus, we need to assume a EoS that describes the microscopic physics of the matter of the NS to complete the system of equations above. Let us assume a realistic modeling of this physics, following the methodology described in [65], where the nuclear matter inside NS are built by joining together different polytropic phases on a sequence of different density intervals (piecewise EoS), given by \(\rho > \rho _0\), if for a set of dividing densities \(\rho _0< \rho _1 < \rho _2\cdots \), the pressure and energy density are everywhere continuous and satisfy

$$\begin{aligned} P(\rho ) = K_i \rho ^{\Gamma _i},\quad d\left( \frac{\epsilon }{\rho } \right) =-P d \left( \frac{1}{\rho }\right) , \quad \rho _{i-1}< \rho < \rho _i,\nonumber \\ \end{aligned}$$
(19)

where \(\epsilon \) is the energy density and fixed by the first law of thermodynamics. The numerical modeling of EoSs is well summarized in [65] (see also [66]). In this work, we will assume the low-density EoSs given by SLy, MS2, GNH3, ALF2, BBB2, H4, PS, WFF3, AP2 and FPS models. In total, we have 10 soft EoSs in our sample. In addition to several key properties of each ones, it should be stressed that these EoSs are consistent with constraints from GW170817 (see [67] and references therein). The best fit parametric configuration that describe these EoSs modeling are summarized in Table III in [65].

On the other hand, if the NS is spinning, its maximum mass can be higher than \(M_{\mathrm{TOV}}\) due to the additional rotational support against gravitational collapse to a BH. In particular, [68] computes the maximum mass allowed by uniform rotation, \(M_{\mathrm{max}}\), purely in terms of the maximum mass of the non-rotating configuration. The importance of this universal relation, connecting the dimensionless spin on the stability line and the maximum dimensionless spin at the mass shedding limit, is that it allows us to calculate \(M_{\mathrm{max}}\) sustainable through fast uniform rotation, finding that for any EoS it is about \(20\%\) larger than the maximum mass supported by the corresponding non-rotating configuration (see [6, 68] for details). Thus, we complement our modeling considering \(M_{\mathrm{max}} = \xi M_{\mathrm{TOV}}\), with \(\xi = 1.203 \pm 0.022\) [6]. Therefore, our full NS model will be constituted from some soft EoSs, rotation effects and chameleon screening. Note that this constraints on \(\xi \) is obtained by assuming GR. In Sect. 3.1, we relax this condition and obtain \(\xi \) within the theoretical framework here investigated. We will find \(\xi =1.14^{+0.080}_{-0.082}\) in the modified gravity we investigated. Thus, these values are fully compatible with each other and this choice will not biased the main results to be presented below.

For the convenience of analysis, let us consider that there is a simple relationship between the expectation mass inside the screening radius in the form \(M(r_{\mathrm{s}})_i = \beta _i M_i\), where i runs over the bodies A and B, in case of a binary system. In principle, we could consider more complicated dependency on the screening mass, for instance, like \(M(r_{\mathrm{s}})_i = \beta _i r_i^n M_i\), where r is the NS radius and n some power law dependency. This type of correction will take PPN corrections predict more complicated relation/dependence on NS mass values. On the other hand, we expect that the screening radius be close to the surface of the star, so the approximation \(M(r_{\mathrm{s}})_i = \beta _i M_i\) serves as a basis to quantify these screening effects within a simple approach, and to investigate possible new effects on the M–R relation. In what follows, we summarize our main results.

Table 1 Summary of the maximum mass values that can be obtained using 10 soft EoSs in different theoretical situations, namely, general relativity, general relativity with rapid uniform rotation effects, chameleon screening and chameleon screening with rapid uniform rotation effects
Full size table

3 Results

Table 1 summarizes the maximum mass of a NS taking into account 4 possible theoretical configurations and 10 input realistic soft EoSs, from which we can quantify how each ingredient is contributing to \(M_{\mathrm{max}}\) estimates. Taking into account the screening mechanisms, we consider a NS with screening radius \(r_s\) and \(M_A\), orbiting a generic body B with \(M_B\), so given the relation previously mentioned, \(M(r_{\mathrm{s}})_i = \beta _i M_i\), we have \(G^{PPN}/G = 1 + {\bar{\alpha }}\), with \({\bar{\alpha }} = 2\alpha ^2(1 - \beta _B - \beta _A + \beta _B \beta _A)\), which is a global constant in our case. In all estimates in Table 1, from the screening and screening + spinning cases, we assume a small and conservative total correction with \({\bar{\alpha }}=-0.05\). As we will see below (subsection A), through a statistical fit, effective corrections on \(G^{PPN}\) are of this order of magnitude. i.e., \({\bar{\alpha }}\) \(\sim -0.05\) at \(\sim 1\sigma \). Corrections \({\bar{\alpha }} > 0\) tends to decrease the expectation value of maximum mass. Figure 1 shows this using SLy EoS + spinning effects, quantifying how much \(\gamma \) and \(G^{\mathrm{PPN}}\) corrections can influence \(M_{\mathrm{max}}\) prediction. Important to note that the expected corrections for \(\gamma \) and \(G^{\mathrm{PPN}}\) within screening gravity are not the same as measured on earth or solar system. Without screening mechanisms, we would have to tune \(\gamma \) value to satisfy terrestrial and solar system experiment bounds, but with screening mechanisms this bound can be automatically satisfied for this screened region with \(M(r_s) \simeq M(r)\) on these scales. Here, we are relaxing this condition in order to have new \(\gamma \) and \(G^{\mathrm{PPN}}\) expected values on a NS under thin-shell effect. Of course, any new \(\gamma \) and \(G^{\mathrm{PPN}}\) expected values should not deviate significantly from solar-system tests. In all our calculations, we are assuming a maximum deviation of 5% when applied to NS structure equations.

Figure 2 shows the NS mass–radius relation from all the EOSs analyzed here, and summarized in Table 1 for spinning + chameleon screening combination, in direct comparison with the expected masses of the PSR J0740 + 6620, \(2.14 \pm 0.10\) \(M\odot \), a NS in a binary system with a white dwarf companion, and the mass of the lighter component of the GW190814 event, \(2.6 \pm 0.10\) \(M\odot \), if we interpret this object as being a NS. It is worth mentioning that in the results shown in Fig. 2 and Table 1, we assume an effective \({\bar{\alpha }}=-0.05\) (corresponding to an effective \(G^{\mathrm{PPN}}/G = 0.95\)).

Fig. 1
figure 1

Illustrative example of a neutron star maximum mass as a function of PPN parameter \(\gamma \) (blue) and \(G^{PPN}\) (orange) assuming EoS SLy

Full size image
Fig. 2
figure 2

Left panel: Neutron star mass–radius relation for different soft EoSs (SLy, MS2, ALF4, GNH3, BBB2) within Spinning + PPN corrections framework, with highlighted regions representing the neutron star mass expectation from the binary system J0740 + 6620 and GW190814. Right panel: Same as in left panel, but from other soft EoSs (H4, PS, WFF3, AP2 and FPS)

Full size image

The consideration above has an observational support. For instance, astrophysical tests on chameleon theories using distance indicators in the nearby universe are bounds to an effective values of \(\Delta G/G \in [0.11{-}0.45]\) [69]. The observation of GW170817 event imposed strong constraints on modified gravity/dark energy scenarios, inferring that the speed of GWs propagation is equal to the speed of light, so limiting the theoretical parametric space for various theories. Chameleon theories in principle has already this physical characteristic and this constraints does not affect in general any aspect in these theories. The authors in [70] calculate that the effects of a plausible variation in \(G_N\) on the period decay of a local binary system are many orders of magnitude suppressed with respect to the effects of a change in \(G_N\). As a result, \({\dot{G}}_N /G_N\) is constrained to \(\sim 10^{-3}\), while \(G_N\) is unconstrained by observations of the binary orbital decay rate. Here, the notation \(G_N\) is the gravitational coupling between two matter sources. Note that our approach is linked to the physical variable \(G_N\) and not \({\dot{G}}_N\). In general, all our considerations about \(G^{\mathrm{PPN}}/G\), i.e., all range of values assumed to the prior are compatible with the bounds derived from other observations.

Fig. 3
figure 3

Statistical reconstruction in the plan \(M^{MG}_{\mathrm{max}} - M^{GR}_{\mathrm{max}}\) using Gaussian processes for machine learning between the maximum mass estimates summarized in Table 1. The solid blue curve is the Gaussian process mean. Both \(M_{\mathrm{max}}\) are in units of \(M_{\odot }\)

Full size image

It is important to note that all EoSs assumed here share a relatively low maximum mass and thus they are called soft [65]. Soft EoSs consistent with GW170817 (see Table 1) are unable to provide enough mass to explain the secondary in GW190814, even for a NS endowed with maximum uniform rotation [5]. We noticed that assuming non-rotating mass configuration, it is troublesome to predicting high mass in GR, and all these EoSs have difficulties in reaching the mass expected by J0740 + 6620.

Thus, in this simplest case, all of these EoSs are excluded. But, adding rotation effects, all these EoSs can explain J0740 + 6620 easily. The authors in [71, 72] were the first to show that spinning up a NS uniformly can increase its mass by up to \(\sim 20\%\). Thus, soft EoSs consistent with rapid uniform rotation is enough to explain that compact object. The analysis frame change the situation if we look at GW190814 secondary component. When considering GR + spinning, all these soft EOSs can not achieve very high masses, and all them should be ruled out. This means that the secondary compact object in GW190814 can not be explained by an EoS like SLy (EoS favored by GW170817). Also, as argued in [5], such EoSs like SLy must now be rejected because their mass-shedding limit is below the lower limit mass of this object in GW190814. But, we find that modifying the hydrostatic equilibrium equations of NSs, incorporating a possible thin-shell effect (in our study from chameleon screening), we can reach larger masses. We note that chameleon screening can achieve high masses and explain GW190814 with soft EoSs like SLy, BB2, ALF4 and H4 (see Fig. 2), which in GR context is not possible. The BBB2, ALF4 and H4 EoSs predictions live in the lower limit of the error bar of GW190814, while the SLy EoS easily archive higher values and can explain GW190814 best fit and more massive objects. Therefore, in presence of new degrees of freedom of gravitational origin like a chameleon screening, it is possible to maintain these soft EoSs to describe the internal structure of these compact objects.

Note that the same EoS in screening + spinning context can explain both NSs in J0740 + 6620 and GW190814 at the same time, changing only the radius prediction. For instance, from EoS SLy analysis, one obtains [9.16–11.27] km and [10.86–11.18] km for GW190814 and J0740 + 6620, respectively.

On the other hand, it is interesting to derive a general relation between the maximum NS mass from GR and the modified gravity theory from some EoS sample. We will obtain this relationship based on our 10 soft EoSs sample [65], where the main results are summarized in Table 1.

To get this relation, let us use a Gaussian Processes (GP). The GP consists of generic supervised learning method designed to solve regression and probabilistic classification problems, where we can interpolate the observations and compute empirical confidence intervals and a prediction in some region of interest [73]. Let us consider that our 10 estimates in Table 1 (the case with rotation effects) are data points to be trained between a certain mass range. Note that 10 data points, from 10 different soft EoSs, give us reasonable information to get a good prediction within the mass range where soft EoS there due the robustness of this methodology. The GP method is the state-of-the-art to obtain statistical information and model prediction from some previously known information or data.

Figure 3 show the data reconstruction. The solid blue curve is the GP mean prediction and the shaded areas is the statistical confidence level (CL) at 68%. Each red point represents the mass information (data point) in the plane \(M^{MG}_{\mathrm{max}} - M^{GR}_{\mathrm{max}}\), in the case with rotation, which has been trained to obtain the reconstruction and prediction between these mass estimates. During the reconstruction, we use the most popular covariance functions in the literature, the standard Gaussian Squared-Exponential. To obtain these results, we make use of some numerical routines available in [74]. Interesting to note that from a soft EoS sample, in general, we can obtain high mass, like a 3 solar mass, within modified gravity approach even at 68% CL.

In Fig. 3, we also show a cubic fit that best fit the plan \(M^{MG}_{\mathrm{max}} - M^{GR}_{\mathrm{max}}\) in the mass range \(M \in [2.10, 2.46]\) \(M_{\odot }\) in GR prediction. We find

$$\begin{aligned} M^{MG}_{\mathrm{max}}= & {} -320.38 + 417.26 M_{\mathrm{max}}\nonumber \\&- 179.96 M^{2}_{\mathrm{max}} + 25.90 M^{3}_{\mathrm{max}}, \end{aligned}$$
(20)

where \(M_{\mathrm{max}}\) is the maximum mass in GR.

Thus, within the theoretical framework investigated here, Fig. 3, summary the most general relations and prediction analysis between \(M^{MG}_{\mathrm{max}}\) and \(M^{GR}_{\mathrm{max}}\) within the soft EoS approach.

Fig. 4
figure 4

Two-dimensional marginalized distributions in the parametric space in the plan \(G^{\mathrm{PPN}}/G - \xi \) at 68% CL and 95% CL

Full size image

3.1 Rotation effects in modified gravity

The spin of the secondary object is a quantity that has not been constrained by the observations of GW190814. Furthermore, we know that the condition for rapid uniform rotation contributes significantly to increase the predicted maximum mass. In the results previously presented, we assume a rotation value that was obtained in GR. Once that the gravity theory has been modified, it is necessary to re-calculate the rotational effect in the new gravity framework. Let us obtain a new estimate on \(\xi \), in the situation we can still explain GW190814 secondary component, if this compact object was a NS, within the screening + rotation framework. For this, we carry out a simple \(\chi ^2\) fit given by

$$\begin{aligned} \chi ^2 = \frac{(M_O - M_{\mathrm{max}}(G^{\mathrm{PPN}}, \xi ))^2}{\sigma ^2_O}, \end{aligned}$$
(21)

where \(M_{\mathrm{max}}(G^{\mathrm{PPN}}, \xi )\) is the maximum mass expected assuming \(G^{\mathrm{PPN}}\), \(\xi \) as free parameter, where we assume the prior \(\xi \in \) [1, 1.2]. The range in \(\xi \) represent non-rotating with \(\xi = 1\) and maximum rotation with \(\xi = 1.2\). This upper prior on \(\xi \) can be interpreted that 20% is probably no longer applicable and this correction is an upper limit. The quantities \(M_O\), \(\sigma ^2_O\), are the mass and the associate error bar at 1\(\sigma \) confidence level (CL) from the GW190814 secondary component observation.

Note that in this approach we are obtaining a new \(\xi \) constraints parametrically through the relation \(M^{MG}_{\mathrm{max}} = \xi M^{MG}_{\mathrm{TOV}}\). Figure 4 shows the parametric space in the plan \(G^{\mathrm{PPN}}/G - \xi \) at 68% CL and 95% CL, after assuming \(\xi \) as a free parameter. We find \(G^{\mathrm{PPN}}/G=0.951^{+0.039}_{-0.041}\) (with the corresponding \(\gamma =0.90 \pm 0.093\)) and \(\xi =1.14^{+0.080}_{-0.082}\). During the analysis we assume SLy as input EoS. Note that this value is statistically compatible with the value inferred in GR. Thus, we conclude that the values and analysis derived in the previous section are not biased from this perspective.

On the other hand, other independent measurements, from other observational perspectives, have been carried out recently reporting, \(\gamma =0.87^{+0.19}_{-0.17}\) [75] and \(|\gamma -1| < 0.2 (\Lambda /100)\) kpc [76], with \(\Lambda =10 {-}200\) kpc. Our constraints are in total agreement with these values. All these constraints suggests no deviation from GR. But, despite the statistical bounds be compatible with GR, note that the gravitational relaxing induced by \(G^{\mathrm{PPN}}\) is enough to generate NS with \(M \sim 2.6\) \(M_\odot \) or more, using a soft EoS. In this case, the NS can rotate below of the maximum limit, which in turn will generate an effective \(G^{\mathrm{PPN}}\) slightly smaller, but statistically equivalent even at 68% CL when considering the maximum value. From Eq. (11), if we assume the coupling parameter to be \(\alpha \simeq 1\) [44, 62], we find \(M(r_s)/M \simeq 0.98\), showing that the screening mass is near to the surface, as expected.

4 Final remarks

We consider that the gravitational interaction can deviate minimally from that predicted from GR, through a PPN correction induced by a thin-shell effect via a chameleon screening. We analyze its impacts on the NS mass–radius relation, motivated to check if it is possible to generate high NS masses (when compared to GR prediction). Our main conclusion is that from a combination of modified gravity, rotation effects and realistic soft EoSs, it is possible to achieve high masses like \(M \sim 2.6\) \(M_\odot \) or more, and explain GW190814 secondary component. In return, it is also possible to explain the most NS massive confirmed to date, i.e., PSR J0740 + 6620. Then one sees here the following interesting aspect: one can not ruled out some soft EoSs without first taking into account effects coming from alternative theories of gravity. Although our theoretical approach be simple, this consequence is clear (see, e.g., Table 1 and Fig. 3). Therefore, gravity can play an important role, without the need to invoke very stiff EoSs or unusual aspects in nuclear matter.

It might be interesting to see how more sophisticated gravity theories behave in this perspective, and if it is possible to obtain new and accurate observational bounds on the free parameters that characterize such scenarios using possible massive NS observations.