1 Introduction

Neutron stars are ideal astrophysical laboratories to probe the nature of nuclear matter under extreme conditions of density and isospin asymmetry [50,51,52, 86, 88] as well as to test fundamental theories of strong-field gravity [72, 97, 100]. While there has been significant improvement in understanding properties of dense nuclear matter at and near nuclear saturation density, \(\rho _0 \simeq 2.5 \times 10^{14}\) g / cm\(^{3}\), our knowledge of dense matter at super-saturation densities \(\rho > \rho _0\) corresponding to the core region of neutron stars remains quite poor. In part, it is related due to the fact that current nuclear interaction models that are fitted to properties of terrestrial nuclear observables largely fail to constrain the isovector part of the nuclear interaction. This, in turn, affects model predictions for properties of nuclear matter with large isospin asymmetry that are present in the core of neutron stars [25, 44, 45]. As a result, different neutron-star matter equation of state models emerge that give rise to very different neutron-star structure properties, such as masses and radii [23, 26], moments of inertia [27, 84], and tidal deformations [22, 30]. Moreover, there is a possible degeneracy between the nuclear matter equation of state and models of gravity applied to describe the structure of neutron stars. While many studies have been devoted to break this degeneracy, it remains one of the outstanding problems to date [5, 13, 20, 38, 49, 55, 69, 79,80,81,82, 96, 98].

In this contribution, we examine effects of the nuclear equation of state uncertainties simultaneously within the general theory of relativity (GR) and the scalar-tensor (ST) theory of gravity. For over a century, GR has been tested in many astrophysical scenarios and thus far its agreement with experiments and observations have been remarkable [82]. Nevertheless, most of the general relativistic tests have been performed in the weak-field regime [72]. Neutron stars exhibit a strong curvature, and it is not yet clear whether gravitational field of such compact objects is fully described by GR. Astrophysical observations may help pin down theories of gravitation in the strong-field regime provided that the equation of state of dense nuclear matter is well known at all densities relevant to neutron stars. In this article, we show that future astrophysical and gravitational wave observations from low- and canonical-mass neutron star will constrain the equation of state of nuclear matter, whereas measurement of radii and tidal deformations of massive neutron stars will aid in putting further constraints on the ST theories of gravity.

In Sect. 2, we present the formalism of calculating the structure of neutron stars in the scalar-tensor theory of gravity. Next, in Sect. 3, we review the equation of state of neutron-star matter within the relativistic mean-field framework and discuss, in particular, the current uncertainties in the isovector sector of the nuclear interaction. In Sect. 4, we present our results for the neutron-star structure calculation in both GR and ST theory of gravity. We demonstrate that the equation of state of nuclear matter can be constrained using canonical- and low-mass neutron-star observations irrespective of the models of gravity used. We also show that once the equation of state is constrained, in conjunction with the massive neutron-star observations, one may then probe and limit the parameter space of the ST theory of gravity.

Throughout this paper, we adopt geometric units, \(c = 1 = G\), where c is the speed of light, and G is the gravitational constant, respectively.

2 Neutron-star structure in scalar-tensor theory of gravity

The scalar-tensor theories of gravitation are one of the most natural generalizations of general relativity that date back to early 1950s [8, 19, 31, 47, 48]. According to these theories, the gravitational force is also mediated by a scalar field, \(\varphi \), in addition to the second-rank metric tensor, \(g_{\mu \nu }\), present in general relativity. One can write the most general form of the action defining this theory as [14, 99]

$$\begin{aligned} S = \frac{1}{16 \pi } \int \mathrm{d}^4x \sqrt{-g_{*}}\left[ R_{*} -2 g_{*}^{\mu \nu } \varphi _{,\mu } \varphi _{, \nu } - V(\varphi )\right] + S_\mathrm{matter}\left[ \psi _\mathrm{matter}; A^2(\varphi )g_{*\mu \nu } \right] , \end{aligned}$$
(1)

where \(R^{*} \equiv g_{*}^{\mu \nu } R^{*}_{\mu \nu }\) is the curvature scalar of the so-called “Einstein metric", \(g_{*\mu \nu }\), describing the pure spin-2 excitations, whereas \(\varphi \) is a long-range scalar field describing spin-0 excitations. Here \(V(\varphi )\) is the scalar field potential, \(\psi _\mathrm{matter}\) is a collective representation of all matter fields, and \(S_\mathrm{matter}\) is the corresponding action of the matter represented by \(\psi _\mathrm{matter}\), which in turn is coupled to the “Jordan–Fierz metric” \(g_{\mu \nu }\) that is related to the “Einstein metric" through the conformal transformation:

$$\begin{aligned} g_{\mu \nu } \equiv A^2(\varphi )g_{*\mu \nu }. \end{aligned}$$
(2)

The field equations are easily formulated in the “Einstein metric”; however, all non-gravitational physical experiments measure quantities in the “Jordan–Fierz metric” and thus they are referred to as the “physical metric”. From here on, we refer to them as the Einstein frame and the Jordan frame, respectively.

Taking variation of the action S with respect to the metric \(g_{*\mu \nu }\) and the scalar field \(\varphi \), we find the set of field equations in the Einstein frame

$$\begin{aligned}&G_{*\mu \nu } = 8 \pi T_{*\mu \nu } + 2 \varphi _{,\mu } \varphi _{, \nu } - g_{*\mu \nu } g_{*}^{\alpha \beta } \varphi _{,\alpha } \varphi _{, \beta } - \frac{1}{2} V(\varphi ) g_{*\mu \nu } \ , \end{aligned}$$
(3)
$$\begin{aligned}&\Box _{*} \varphi - \frac{1}{4} \frac{\mathrm{d} V(\varphi )}{\mathrm{d}\varphi } = - 4 \pi \alpha (\varphi ) T_{*}. \end{aligned}$$
(4)

Here again \(T_{*\mu \nu }\) is the energy–momentum tensor in the Einstein frame, which is related to the physical energy–momentum tensor in Jordan frame \(T_{\mu \nu }\) via

$$\begin{aligned} T_{*\mu \nu } = A^2(\varphi ) T_{\mu \nu } \end{aligned}$$
(5)

and \(T_{*}\) and \(\alpha (\varphi )\) are defined as

$$\begin{aligned}&T_{*} \equiv T_{*\mu }^{\mu } \ , \end{aligned}$$
(6)
$$\begin{aligned}&\alpha (\varphi ) \equiv \frac{\mathrm{d} \ln A(\varphi )}{\mathrm{d} \varphi }. \end{aligned}$$
(7)

The energy–momentum conservation can either be expressed as \(\nabla _{\nu } T_{\mu }^{\nu } = 0\) in the Jordan frame, or

$$\begin{aligned} \nabla _{*\nu } T_{*\mu }^{\nu } = \alpha (\varphi ) T_{*} \nabla _{*\mu } \varphi , \end{aligned}$$
(8)

in the Einstein frame.

We consider the neutron star to be made of a perfect fluid. In this case, the energy density \(\mathcal {E}\), the pressure P and the 4-velocity \(u_{\mu }\) in the two frames are related via the following relations:

$$\begin{aligned}&\mathcal {E}_{*} = A^4(\varphi ) \mathcal {E}, \end{aligned}$$
(9)
$$\begin{aligned}&P_{*} = A^4(\varphi ) P, \end{aligned}$$
(10)
$$\begin{aligned}&u_{*\mu } = A^{-1}(\varphi ) u_{\mu }, \end{aligned}$$
(11)

where the subscript asterisks denote quantities in the Einstein frame, as usual. The space-time metric describing an unperturbed, non-rotating, spherically symmetric neutron star can be written as

$$\begin{aligned} \mathrm{d}s^2_{*} = -e^{2 \nu (r)} \mathrm{d}t^2 + e^{2 \lambda (r)} r^2 + r^2 (\mathrm{d}\theta ^2 + \sin ^2 \mathrm{d}\phi ^2), \end{aligned}$$
(12)

where

$$\begin{aligned} e^{-2 \lambda (r)} = 1 - \frac{2 M(r)}{r} \equiv N^2(r), \end{aligned}$$
(13)

and the function \(\nu (r)\) will be calculated later. Here we also defined \(N(r) \equiv \left( 1-2M(r)/r \right) ^{1/2}\), for convenience. Using the space-time metric above, one can arrive at the dimensionally reduced field equations [99], which in turn can be cast into equations for hydrostatic structure:

$$\begin{aligned} \frac{\mathrm{d}M(r)}{\mathrm{d}r}= & {} 4 \pi r^2 A^4(\varphi ) \mathcal {E}(r) + \frac{r^2}{2} N^2(r) \chi ^2(r) + \frac{r^2}{4} V(\varphi ), \end{aligned}$$
(14)
$$\begin{aligned} \frac{\mathrm{d}P(r)}{\mathrm{d}r}= & {} - (\mathcal {E}(r) + P(r)) \bigg (\frac{\mathrm{d} \nu (r)}{\mathrm{d}r} + \alpha (\varphi ) \chi (r)\bigg )\ , \end{aligned}$$
(15)
$$\begin{aligned} \frac{\mathrm{d}\nu (r)}{\mathrm{d}r}= & {} \frac{M(r)}{r^2 N^2(r)} + \frac{4 \pi r}{N^2(r)} A^4(\varphi ) P(r) + \frac{1}{2} r \chi ^2(r) -\frac{r}{4 N^2(r)} V(\varphi ), \end{aligned}$$
(16)
$$\begin{aligned} \frac{\mathrm{d} \chi (r)}{\mathrm{d}r}= & {} \bigg \{4 \pi A^4(\varphi ) \bigg [\alpha (\varphi ) (\mathcal {E}(r) - 3 P(r)) + r \chi (r) (\mathcal {E}(r) - P(r))\bigg ] \nonumber \\- & {} \frac{2 \chi (r)}{r} \left( 1-\frac{M(r)}{r}\right) + \frac{1}{2}r\chi (r)V(\varphi ) + \frac{1}{4}\frac{\mathrm{d}V(\varphi )}{\mathrm{d}\varphi } \bigg \}N^{-2}(r), \end{aligned}$$
(17)

where we introduced

$$\begin{aligned} \chi (r) \equiv \frac{\mathrm{d} \varphi }{\mathrm{d}r}. \end{aligned}$$
(18)

It is well known that predictions of scalar-tensor theories are physically equivalent to nonlinear modified gravity theories [58]. In particular, for the \(R^2\)-gravity, where the Ricci scalar is replaced with \(f(R) = R + a R^2\) in the Einstein–Hilbert action, the explicit form of the potential \(V(\varphi )\) can be written as

$$\begin{aligned} V(\varphi ) = \frac{1}{4a} (1-e^{-2\varphi }{\sqrt{3}})^2 \ , \quad \alpha (\varphi ) = -\frac{1}{\sqrt{3}}. \end{aligned}$$
(19)

Next, following Ref. [14], we set \(V(\varphi )=0\), and consider a coupling function of the form

$$\begin{aligned} A(\varphi ) = \exp \left( \alpha _0\varphi + \frac{1}{2}\beta _0 \varphi ^2\right) \ . \end{aligned}$$
(20)

The coupling constants \(\alpha _0\) and \(\beta _0\) are real numbers. It was shown that measurement of the surface atomic line redshifts from neutron stars could be used as a direct test of strong-field gravity theories [17]. In particular, coupling constants with \(\alpha _0=0\) and \(\beta _0=-8\) were used. Over the past decade, significant improvements have been made in constraining these coupling constants. For example, solar-system experiments, binary-pulsar and pulsar-white dwarf timing observations put some of the most stringent constraints on these constants that conservatively can be written as [14, 33]

$$\begin{aligned}&\alpha _0^2 = \partial \ln A(\varphi _0) / \partial \varphi _0 \lesssim 2 \times 10^{-5} \ , \end{aligned}$$
(21)
$$\begin{aligned}&\beta _0 = \partial ^2 \ln A(\varphi _0) / \partial \varphi _0^2 \gtrsim -5 \ . \end{aligned}$$
(22)

On the other hand, it was first shown in Ref. [14] that predictions by models with \(\beta _0>-4.35\) cannot in general be distinguished from the general relativistic results due to the so-called “spontaneous scalarization” effect. Notice that GR is automatically recovered in the limits of \(\alpha _0 = 0\) and \(\beta _0 = 0\), where the Eqs. (14–17) are reduced to the famous Tolman–Oppenheimer–Volkoff equations [65, 92].

We solve the interior and the exterior problem simultaneously using the following natural boundary conditions in the center of the star:

$$\begin{aligned} P(0) = P_\mathrm{c}, \quad \mathcal {E}(0) = \mathcal {E}(\mathrm c), \quad \chi (0) = 0. \end{aligned}$$
(23)

We also demand cosmologically flat solution at infinity to agree with the observation:

$$\begin{aligned} \lim _{r \rightarrow \infty } \nu (r) = 0, \quad \lim _{r \rightarrow \infty } \varphi (r) = 0. \end{aligned}$$
(24)

The only input required to integrate Eqs. (14–17) is the equation of state of neutron-star matter, \(\mathcal {E} = \mathcal {E}(P)\), that is trivial in the case of exterior solution. Once an EOS is supplemented, for a given central pressure \(P(0) = P_\mathrm{c}\), one can integrate them from the center of the star \(r = 0\), all the way up to \(r \rightarrow \infty \).

The stellar coordinate radius is then determined by the condition in which the pressure vanishes, i.e.\(P(r_\mathrm{s}) = 0\), where \(r_\mathrm{s}\) is the surface radius in the Einstein frame. The physical radius of a neutron star R is then found in the Jordan frame through

$$\begin{aligned} R = A\left( \varphi (r_\mathrm{s})\right) r_\mathrm{s}. \end{aligned}$$
(25)

The physical stellar mass M as measured by an observer at infinity—also known as the Arnowitt–Deser–Misner (ADM) mass—matches with the coordinate mass since at infinity the coupling function approaches unity.

3 Neutron-star matter equation of state

The structure of neutron stars is sensitive to the equation of state of cold, fully catalyzed, and neutron-rich matter. The matter inside neutron stars span many orders of magnitude in density leading to rich and exotic phases in their interiors. In the outer crust of neutron stars, the matter is organized into a Coulomb lattice of neutron-rich nuclei embedded in a degenerate electron gas [6, 75]. The nuclear composition in this region is solely determined by the masses of neutron-rich nuclei in the region of \(26 < Z \lesssim 40\) and the pressure support is primarily provided by the degenerate electrons. The equation of state for this region is, therefore, relatively well known [6, 36]. In this work, we adopt the outer crust equation of state by Haensel, Zdunik and Dobaczewski [36]. As the density increases, the nuclei in the outer crust become more and more neutron rich. At a density of about \(\rho \approx 4 \times 10^{11}\) g / cm\(^{-3}\), the nuclei in the outer crust become so neutron rich that they can no longer hold additional neutrons, and neutrons start dripping out. This region defines the boundary between the outer and the inner crust.

The inner crust extends from the neutron-drip density up to about \(\rho \approx 2/3 \rho _0\), where the uniformity in the system is restored. On the top layers of the inner crust, nucleons continue to cluster into a Coulomb crystal of neutron-rich nuclei embedded in a uniform electron gas, where, however, now the system is also in a chemical equilibrium with a superfluid neutron gas [70]. As the density continues to increase, the spherical nuclei start to deform in an effort to reduce the Coulomb repulsion. As a result, the inner crust exhibits complex and exotic structures that are collectively known as “nuclear pasta" [37, 57, 73], which emerge from a dynamical competition between the short-range nuclear attraction and the long-range Coulomb repulsion. Although significant progress has been made in simulating this exotic region [24, 41,42,43], the equation of state for this region remains highly uncertain and must be inferred from theoretical calculations. While a detailed knowledge of the equation of state for this region is important for the interpretations of cooling observations [21], its impact on the bulk properties of neutron stars is minimal [71]. For this region, therefore, we resort to the equation of state provided by Negele and Vautherin [63].

The structure of neutron stars are mostly sensitive to the equation of the state of the core, and the crust–core transition properties. In particular, most of the mass of neutron stars are contained in the liquid core that comprises a region of the star with densities as low as one-third to as high as ten times nuclear matter saturation density \(\rho _0\). For this region, we employ the equation of state generated from various refinements of relativistic mean-field model by Serot and Walecka [76, 77, 95]. For consistency, the transition density from the liquid core to the solid crust is computed using the same relativistic mean-field models, where it is done by searching for the critical density at which the uniform system becomes unstable to small amplitude density oscillations [11]. We would like to emphasize that the crust–core transition density (hence transition pressure) plays an important role in neutron star bulk properties with various models leading to crust thicknesses that differ by over 1 km and predicting significantly different crustal components of the moment of inertia that are important in interpreting observations of pulsar glitches [71].

The equation of state for the uniform liquid core is based on an interaction Lagrangian that has been accurately calibrated to a variety of ground-state properties of both finite nuclei and infinite nuclear matter. This model includes a nucleon field (\(\psi \)), a scalar–isoscalar meson field (\(\phi \)), a vector–isoscalar meson field (\(V^{\mu }\)), and an isovector meson field (\(b^{\mu }\)) [76, 77]. The free Lagrangian density for this model is given by [23, 27]

$$\begin{aligned} {\mathcal {L}}_{0}= & {} \bar{\psi }\left( i\gamma ^{\mu }\partial _{\mu }\!-\!m_b\right) \psi + \frac{1}{2}\partial _{\mu } \phi \partial ^{\mu } \phi -\frac{1}{2}m_{s}^{2}\phi ^{2} \nonumber \\&- \frac{1}{4}F^{\mu \nu }F_{\mu \nu } + \frac{1}{2}m_{v}^{2}V^{\mu }V_{\mu } - \frac{1}{4}{} \mathbf{b}^{\mu \nu }\mathbf{b}_{\mu \nu } + \frac{1}{2}m_{\rho }^{2}{} \mathbf{b}^{\mu }{} \mathbf{b}_{\mu } \;, \end{aligned}$$
(26)

where \(F_{\mu \nu }\) and \(\mathbf{b}_{\mu \nu }\) are the nuclear isoscalar and isovector field tensors, respectively,

$$\begin{aligned} F_{\mu \nu }= & {} \partial _{\mu }V_{\nu } - \partial _{\nu }V_{\mu } \;, \end{aligned}$$
(27)
$$\begin{aligned} \mathbf{b}_{\mu \nu }= & {} \partial _{\mu }{} \mathbf{b}_{\nu } - \partial _{\nu }{} \mathbf{b}_{\mu } \;, \end{aligned}$$
(28)

where the parameters \(m_b\), \(m_{s}\), \(m_{v}\), and \(m_{\rho }\) represent the nucleon and meson masses and may be treated as empirical constants. The interacting component of Lagrangian density can be written by the following expression [62, 76, 77]:

$$\begin{aligned} {\mathcal {L}}_\mathrm{int} = \bar{\psi }\left[ g_\mathrm{s}\phi \!-\! \left( g_\mathrm{v}V_\mu \!+\! \frac{g_{\rho }}{2}\varvec{\tau }\cdot \mathbf{b}_{\mu } \right) \gamma ^{\mu } \right] \psi - U(\phi ,V^{\mu },\mathbf{b^{\mu }}) \end{aligned}$$
(29)

that includes Yukawa couplings—with coupling parameters, \(g_\mathrm{s}\) , \(g_\mathrm{v}\), and \(g_{\rho }\)—between the nucleon and meson fields. The Lagrangian density is also supplemented by nonlinear meson interactions, \(U(\phi ,V^{\mu },\mathbf{b}^{\mu })\) that improve the phenomenological standing of the model:

$$\begin{aligned} U(\phi ,V^{\mu },\mathbf{b}^{\mu }) = \frac{\kappa }{3!} (g_\mathrm{s}\phi )^3 \!+\! \frac{\lambda }{4!}(g_\mathrm{s}\phi )^4 \!-\! \frac{\zeta }{4!} (g_\mathrm{v}^2 V_{\mu }V^\mu )^2 \!-\! \Lambda _\mathrm{v} (g_{\rho }^{2}\,\mathbf{b}_{\mu }\cdot \mathbf{b}^{\mu }) (g_\mathrm{v}^2V_{\nu }V^\nu ). \end{aligned}$$
(30)

The details on the calibration procedure can be found in Refs. [12, 45, 76, 77, 91] and references therein.

While the full complexity of the quantum system cannot be tackled exactly, the ground-state properties of the system may be computed in a mean-field approximation. In this approximation, all the meson fields are replaced by their classical expectation values and their solution can be readily obtained by solving the classical Euler–Lagrange equations of motion. The only remnant of quantum behavior is in the treatment of the nucleon field which emerges from a solution to the Dirac equation in the presence of appropriate scalar and vector potentials [76, 77]. Following standard mean-field practices, the energy density of the system is given by the following expression:

$$\begin{aligned} {\mathcal E} (\rho , \alpha )= & {} \frac{1}{\pi ^{2}}\int _{0}^{k_\mathrm{F}^{p}} k^{2}E_{k}^{*}\,\mathrm{d}k + \frac{1}{\pi ^{2}}\int _{0}^{k_\mathrm{F}^{n}} k^{2}E_{k}^{*}\,\mathrm{d}k + \frac{1}{2}m_{s}^{2}\phi _{0}^{2} + \frac{\kappa }{3!} (g_\mathrm{s}\phi _{0})^3 + \frac{\lambda }{4!}(g_\mathrm{s}\phi _{0})^4 \nonumber \\&+ \frac{1}{2}m_{v}^{2}V_{0}^{2} + \frac{\zeta }{8}(g_\mathrm{v}V_{0})^4 + \frac{1}{2}m_{\rho }^{2}b_{0}^{2} + 3\Lambda _\mathrm{v}(g_\mathrm{v}V_{0})^2 (g_{\rho }b_{0})^2, \end{aligned}$$
(31)

where \(\rho \) is the baryon density of the system, \(\alpha = ( \rho _\mathrm{n} - \rho _\mathrm{p})/\rho \) is the neutron–proton asymmetry, \(E_{k}^{*}\!=\!\sqrt{k^{2}+m_\mathrm{b}^{*2}}\), \(m_\mathrm{b}^{*}\!=\!m_\mathrm{b}-g_\mathrm{s}\phi _{0}\) is the effective nucleon mass, \(k_\mathrm{F}^{p} (k_\mathrm{F}^{n})\) is the proton (neutron) Fermi momentum. Since the mean-field approximation is thermodynamically consistent, the pressure of the system at zero temperature may be obtained either directly from the energy–momentum tensor or from the energy density and its first derivative [76, 77], that is,

$$\begin{aligned} P (\rho , \alpha ) = \rho \frac{\partial \mathcal {E}(\rho , \alpha )}{\partial \rho } -\mathcal {E}(\rho , \alpha ). \end{aligned}$$
(32)

It is often useful to expand the energy per nucleon of the system in even powers of \(\alpha \):

$$\begin{aligned} E/A(\rho , \alpha ) - m_b = E_\mathrm{SNM}(\rho ) + \alpha ^2 S(\rho ) + \mathcal {O}(\alpha ^4), \end{aligned}$$
(33)

where \(E_\mathrm{SNM}(\rho )\) is the energy per nucleon of symmetric nuclear matter, whereas \(S(\rho )\) is referred to as the nuclear symmetry energy, a quantity that represents the increase in the energy of the system as it departs from the symmetric limit of equal number of neutrons and protons [40, 93].

Further, it is customary to characterize the behavior of both symmetric nuclear matter and the symmetry energy in terms of a few bulk parameters near nuclear saturation density \(\rho _0\):

$$\begin{aligned}&E_\mathrm{SNM}(\rho ) = \varepsilon _0 + \frac{1}{2} K x^2 + \cdots , \end{aligned}$$
(34)
$$\begin{aligned}&S(\rho ) = J + L x + \frac{1}{2} K_\mathrm{sym} x^2 + \cdots , \end{aligned}$$
(35)

where \(x \equiv (\rho -\rho _0)/3\rho _0\) is a dimensionless parameter, \(\varepsilon _0\) and K represent the energy per nucleon and the incompressibility coefficient of symmetric nuclear matter, respectively, whereas J, L, and \(K_\mathrm{sym}\) are the magnitude, slope and curvature of the symmetry energy at saturation. The bulk parameters of symmetric nuclear matter are relatively well constrained [22, 28, 28]. On the other hand, the density dependence of the nuclear symmetry energy remains unconstrained due to lack of sensitive isovector nuclear probes [29, 54]. In this contribution, we will primarily concentrate on the impact of variations of density slope of the symmetry energy L that is closely related to the pressure of pure neutron matter at saturation density.

We assume the neutron-star matter to consist of neutrons, protons, electrons, and muons in chemical equilibrium. We do not consider any “exotic” degrees of freedom, such as hyperons, meson condensates, or quarks. The electrons and muons are assumed to behave as relativistic free Fermi gases (with \(m_{e}\!\equiv \!0\)). The muons appear in the system only after the electronic Fermi momentum becomes equal to the muon rest mass. The total energy density and pressure of the star are obtained by adding up the nucleonic and leptonic contributions.

4 Results

In this work, we use the FSUGold2 parametrization that was introduced by Ref. [12] to specifically apply for both finite nuclei and neutron stars. Of particular importance to this study is the role of omega-meson self-interactions, as described by the parameter \(\zeta \) in the interaction Lagrangian density, that is used to tune the equation of state at high density to reproduce the maximum mass of a neutron star. It was first shown by Müller and Serot that using different values of \(\zeta \) one can reproduce the same observed nuclear matter properties at nuclear saturation, yet produce maximum neutron-star masses—using general relativity only—that differ by almost one solar mass [62]. For example, models with \(\zeta \!=\!0\) predict the maximum neutron-star masses of about \(2.8 M_{\odot }\). Note that this tuning primarily affects the equation of state of symmetric nuclear matter at high density, which is relevant to the core of neutron stars. On the other hand, by including the nonlinear coupling constant \(\Lambda _\mathrm{v}\), Horowitz and Piekarewicz showed that one can modify the density dependence of the symmetry energy [45]. It was shown that tuning \(\Lambda _\mathrm{v}\) provides a simple and efficient method of controlling the density dependence of symmetry energy without compromising the success of the model in reproducing well-determined ground-state observables. The original FSUGold2 model has a relatively stiff symmetry energy with the density slope of \(L = 112.8\) MeV. Following the same method as outlined in Ref. [27], we obtain a family of “FSUGold2" parametrizations with \(L = 47\), 60, 80, 100 MeV. We emphasize that in doing so, predictions for properties of finite nuclei, such as the binding energies and charge radii of closed shell nuclei, remain intact.

The maximum mass of a neutron star predicted by the original FSUGold2 with \(\zeta = 0.0256\) and using the general relativistic TOV equations is \(M_\mathrm{max} = 2.07 M_{\odot }\), which is consistent with the observations of highly precise measurements of two massive neutron stars made at the Green Bank Telescope [4, 18]. Indeed, the maximum possible mass of a neutron star may not be very far from this value as was recently shown by Refs. [60, 61, 74], \(M_\mathrm{max} \approx 2.17 M_\mathrm{Sun}\). This in turn suggests that the \(\zeta \) parameter that controls the stiffness of the equation of state of symmetric nuclear matter is already well constrained, and future observations of maximum mass may put even tighter constraints.

On the other hand, the original FSUGold2 model predicts the corresponding general relativistic radius of a canonical 1.4 solar mass neutron star to be \(R =14.11\) km. Unfortunately, direct determination of the neutron-star radii at present is not quite satisfactory. Early attempts by Özel and collaborators to determine simultaneously the mass and radius of three X-ray bursters resulted in stellar radii to be between 8 and 10 km [66]. Later, Steiner et al. supplemented Özel’s study with additional neutron stars and concluded that the most probable radius of a 1.4 \(M_{\odot }\) lies in the range of 10.4–12.9 km [86, 87]. Nevertheless, this more conservative estimate has been put into question by Suleimanov et al., who used a more complete model of the neutron-star atmosphere and obtained a radius greater than 14 km for a single source studied [89]. Recognizing this situation and the many challenges posed by the study of X-ray bursters, Guillot et al. concentrated on the determination of stellar radii by studying quiescent low-mass X-ray binaries (qLMXB) in globular clusters. By explicitly stating all their assumptions, in particular the common radius assumption, they were able to determine a rather small neutron-star radius of \(9.4 \pm 1.2\) km [34, 35]. And recently, an improved precision over previous measurements was obtained by incorporating distance uncertainties to the globular cluster M13 that suggests the neutron-star radius to be \(12.3^{+1.9}_{-1.7}\) km and/or \(15.3^{+2.4}_{-2.2}\) km depending on the composition of the atmosphere being as H or He, respectively [78].

Since the density slope of the symmetry energy L is related to the pressure of pure neutron matter at saturation density, and the pressure of the neutron-rich matter is strongly correlated with the neutron-star radius [9, 28], we build a family of the FSUGold2 parametrizations using different values of L. While these parametrizations provide an accurate description of ground-state properties of finite nuclei, they also predict stellar radii that differ by over 1 km. The uncertainties in the density dependence of the symmetry energy exhibited by the RMF models here broadly bracket the uncertainties stemming from various nuclear many-body models. This includes from non-relativistic mean-field Skyrme interactions to microscopic calculations of the equation of state that are based on the nucleon–nucleon interactions and consider other degrees of freedom such as pions, \(\Delta \)-resonances, and hyperons [10, 25, 56, 83]. Note, however, that none of our parametrizations can produce neutron-star radii that are smaller than 12 km, which remains one of the biggest challenges today [26, 46].

Fortunately, the prospects for precision neutron-star mass–radius measurements have never been better, especially with the upcoming Neutron Star Interior Composition Explorer (NICER) X-ray timing mission that can put some of the tightest constraints in the near future. On the other hand, the recent gravitational wave observation from a binary neutron-star merger [1] has already put some indirect constraint on the neutron-star radii. In particular, using the upper limit on the tidal deformability measurement determined by LIGO and Virgo Collaboration, Refs. [3, 30, 90] have significantly constrained the allowed EoS models, which in turn predicted a radius of a canonical neutron star to be smaller than about 14 km. The tidal deformability is an intrinsic neutron-star property that describes the tendency of a neutron star to develop a mass quadrupole as a response to the tidal field induced by its companion [16, 32]. The dimensionless tidal polarizability \(\Lambda \) is defined as follows:

$$\begin{aligned} \Lambda = \frac{2k_{2}}{3}\left( \frac{R}{M}\right) ^{5} \;, \end{aligned}$$
(36)

where \(k_{2}\) is the second Love number [7, 15]. The tidal deformability is highly sensitive to the stellar radius, \(\Lambda \!\sim \!\!R^{5}\) and, therefore, its measurement can be used as a proxy to constrain the neutron-star radius that has been notoriously difficult to measure in the past [34, 35, 39, 53, 64, 66, 67, 85, 86, 89].

Table 1 Predictions for the masses and radii of neutron stars in both general relativity and scalar-tensor theory of gravity using a family of FSUGold2 interaction whose isovector coupling constants are tuned to give different values of the density slope of the symmetry energy L
Full size table
Fig. 1
figure 1

The mass-versus-radius relation calculated using the family of FSUGold2 parametrizations considered in this work labeled with the corresponding L-values. Here M stands for either general relativistic TOV-mass (solid lines) or the ADM-mass (dashed lines) predicted by the scalar-tensor theory. For the coupling constants of scalar-tensor theory an upper observational bounds of \(\alpha _0 = \sqrt{2.0} \times 10^{-5}\) and \(\beta _0 = -5.0\) are used

Full size image

In Table 1, we provide with our calculations for neutron-star masses and radii predicted by both general relativity and the scalar-tensor theory of gravity. Notice that the variations of the density slope of the symmetry energy as depicted in various L-values predict neutron-star radii that are different by over 1 km. Nevertheless, both general relativity and scalar-tensor theory predict very similar stellar radii for most neutron stars, unless their mass is close to two solar mass. This suggests that measurements of canonical neutron-star radii, in particular, would be unable to constrain the parameters of the scalar-tensor theory but would be extremely useful in constraining the equation of state of nuclear matter. Moreover, given that the density slope of the symmetry energy is relatively insensitive to the maximum stellar mass which mostly constrains the equation of state of symmetric nuclear matter at high densities, one can use radii measurements to place significant constraints primarily on the density dependence of the nuclear symmetry energy. For massive neutron stars, the radii predictions in two models of gravity start to deviate. This phenomenon of the so-called “spontaneous scalarization” was first observed by Damour and Esposito-Farèse [16].

For completeness, in Fig. 1, we also show the full mass-versus-radii relation predicted in both GR and ST. Since recent measurements of tidal deformability suggested that neutron stars with radii \(R \gtrsim 14\) km may have already been ruled out [2, 3, 30, 59, 61, 90, 94], in this figure, we do not display predictions from the original FSUGold2 parametrization. It is safe to say that given that multi-messenger era of gravitational wave astronomy and X-ray observations of neutron stars is in its infancy, future observations of neutron stars in mass-range of \(M \lesssim 1.4 M_{\odot }\) will undoubtedly put tighter constraints on the nuclear equation of state (below \(\rho \approx 2.5 \rho _0\) corresponding to the central density of a neutron star) but not on the coupling constants of the scalar-tensor theory.

It is particularly interesting to compare the radii of two-solar mass neutron stars predicted by both theory of gravitation. As an example, we use the soft equation of state with \(L=47\) MeV and find that stellar radii differ by 5.6% (see Table 1). Since massive neutron stars are rarely found in nature—lots of theoretical work and sufficient observational data may be required to distinguish general relativity from the scalar-tensor theory using radii observations from massive stars alone—we note that this may not be the case if one considers tidal deformability measurements which scales as \(R^5\). In this contribution, we did not directly calculate tidal deformabilities in the scalar-tensor theory that in addition to radii also depends on the tidal Love number \(k_2\) [68]. It has been shown, however, that within the general relativistic framework \(k_2\) is less sensitive to the radius of a neutron star, and the \(R^5\) scaling behavior in tidal deformability is quite robust [30]. It is worth to point out that the ratio \((R^\mathrm{ST}_{2.0}/R^{GR}_{2.0})^5 = 1.312\) suggests that the corresponding tidal deformabilities in two models may differ by over 30%. We will explore this in more details in future work. This is intriguing because while mass–radii measurements of low-mass neutron stars will constrain the equation of state, tidal deformability measurements of massive stars would enable to test models of gravity in strong regime.

5 Conclusions

We have examined bulk properties of neutron stars in general relativity and the scalar-tensor theory of gravity. In particular, we discussed the current uncertain role of the density dependence of the nuclear symmetry energy in the determination of the equation of state of neutron-star matter. Our analysis shows that future observations of neutron-star radii with masses \(M \lesssim 1.4 M_{\odot }\) will enable to primarily constrain the equation of state of dense nuclear matter but not the coupling constants of the scalar-tensor theory of gravity. In particular, we have found that in conjunction with the maximum stellar mass measurements, this will lead to placing tight constraints on the density dependence of the nuclear symmetry energy.

Recently, it had been discussed that future observations of massive neutron stars may constrain the maximum sound velocity as well as the coupling parameter in the scalar-tensor theory of gravity [82]. While we confirm this result, we would also like to emphasize that once the equation of state is constrained using canonical- and low-mass neutron-star observations, and after a consensus has been reached on the upper limit of the neutron-star mass, only then one can place useful constraints on the coupling parameters of the scalar-tensor gravity. We would like to note that our work does not take into account the possibility that the core of neutron stars may have exotic degrees of freedom, such as hyperons, meson condensates, or quarks. In particular, the central density in neutron stars with \(M \gtrsim 1.4 M_{\odot }\)—for which spontaneous scalarization occurs—reaches nuclear densities of \(\rho _\mathrm{c} \gtrsim 2.5 \rho _0\), where our knowledge of strong interaction is limited. Certainly there could still be a degeneracy between the scalar-tensor model of gravity and models of strong interaction at high densities. Much collaborative theoretical and observational efforts of both nuclear physics and gravitational physics community are, therefore, required on this front.