Tidal Deformability of Neutron Stars in Scalar-tensor Theories of Gravity

Gravitational waves from compact binary coalescences are valuable for testing theories of gravity in the strong field regime. By measuring neutron star tidal deformability using gravitational waves from binary neutron stars, stringent constraints were placed on the equation of state of matter at extreme densities. Tidal Love numbers in alternative theories of gravity may differ significantly from their general relativistic counterparts. Understanding exactly how the tidal Love numbers change will enable scientists to untangle physics beyond general relativity from the uncertainty in the equation of state measurement. In this work, we explicitly calculate the fully relativistic l ≥ 2 tidal Love numbers for neutron stars in scalar-tensor theories of gravitation. We use several realistic equations of state to explore how the mass, radius, and tidal deformability relations differ from those of general relativity. We find that tidal Love numbers and tidal deformabilities can differ significantly from those in general relativity in certain regimes. The electric tidal deformability can differ by ∼200%, and the magnetic tidal deformability differs by ∼300%. These deviations occur at large compactnesses (C = M/r ≳ 0.2) and vary slightly depending on the equation of state. This difference suggests that using the tidal Love numbers from general relativity could lead to significant errors in tests of general relativity using the gravitational waves from binary neutron star and neutron star black hole mergers.


INTRODUCTION AND MOTIVATION
Compact objects such as neutron stars and black holes are essential for testing general relativity (GR) in the strong field regime.Gravitational waves emitted by compact objects by LIGO-Virgo have improved our understanding of gravity in the strong field regime.The LIGO/Virgo collaboration has detected almost one-hundred compact binary coalescences to date: two binary neutron star mergers, two neutron star-black hole mergers, and more than eighty binary black hole mergers Abbott et al. (2019aAbbott et al. ( , 2021a,b),b).An independent analysis of the available data found even more events Nitz et al. (2020Nitz et al. ( , 2021a,b),b).Analysis of these events has already placed limits on possible deviations from GR Abbott et al. (2021cAbbott et al. ( , 2019b)); Mehta et al. (2022); Wang et al. (2021Wang et al. ( , 2022)).Recently, waveforms for various alternate theories of gravity have been developed and applied to parameter estimation.These waveforms allow for stringent tests of various theories of gravity and more general tests for physics beyond GR such as scalar and tensor propagation modes Chatziioannou et al. (2021); Nair et al. (2019); Mirshekari & Will (2013).
Neutron stars are also unique laboratories for studying nuclear physics at ultra-high densities.Information about neutron star matter is encoded in gravitational wave (GW)s from binary neutron star and neutron star-black hole mergers Özel & Freire (2016); Hebeler et al. (2010Hebeler et al. ( , 2013)).Neutron stars contain vital information needed to understand phases of matter encountered in Quantum Chromodynamics.The tidal deformability encodes information about the nuclear equation of state in GWs Binnington & Poisson (2009); Damour & Nagar (2009); Hinderer (2008).Studies of binary neutron star merger GW170817 have improved our knowledge of the nuclear equation of state Raaijmakers et al. (2021); Capano et al. (2020); Radice & Dai (2019); Abbott et al. (2018).Despite this, the nuclear equation of state is still unknown.Studying neutron stars in alternative theories of gravity is challenging because deviations in neutron star properties caused by non-GR effects are of the same order of magnitude as the stephanie.brown@aei.mpg.dearXiv:2210.14025v2[gr-qc] 22 Sep 2023 uncertainty in the equation of state.Understanding how the mass-radius-tidal deformability relationships deviate from GR is essential to untangling these differences.
Tidal deformability connects GWs and the nuclear equation of state.Tidal deformabilities and the associated tidal Love numbers relate an applied external tidal field to the induced internal multipole moment, measuring the magnitude of deformation under a given tidal force.Love numbers were initially defined in Newtonian gravity Love (1909); Shida (1912) and then expanded to GR by Flanagan & Hinderer (2008); Hinderer (2008).The concept was further expanded and made more concrete in several follow-up papers, including Binnington & Poisson (2009) and Damour & Nagar (2009).
This work focuses on scalar-tensor theory, one of the most natural and best studied alternate theories of gravity.The theory was initially motivated partly by Mach's principle Brans & Dicke (1961) and partly in an attempt to expand GR to five dimensions Jordan (1955).However, it is still of interest today.Scalar degrees of freedom are critical for string theory, superstring theory, and other supergravity theories Fujii & Maeda (2003).Therefore scalar-tensor theories can sometimes be used as a phenomenological proxy for more complex extensions of GR.Furthermore, scalar fields have been proposed as an alternative solution to the dark energy problem Garcí-Bellido & Quirós (1990); Boisseau et al. (2000); Clifton et al. (2012).
Scalar-tensor theories add a massless scalar field (φ) to the standard GR metric (g µν ).The metric and the scalar field are coupled into an 'effective metric' gµν = A 2 (φ)g µν .The earliest versions of this theory were presented more than half a century ago by Fierz (1956); Jordan (1955); Brans & Dicke (1961).In the simplest scalar-tensor theory, known as FJBD (Fierz, Jordan, and Brans and Dicke), the scalar field is coupled to the metric by the coupling function A(φ) = e αφ .Solar system experiments have placed stringent constraints on the value of α Shapiro (1990).These constraints also significantly limit the strong-field behavior.Damour and Esposito-Ferèse discovered the 'spontaneous scalarization' effect, which allows large deviations from GR in the strong field regime without violating the strict solar system constraints.Damour and Esposito-Ferèse defined A(φ) = e βφ 2 /2 and found that scalarization occurs for β ≲ −4. 5 Damour & Esposito-Farèse (1993).A follow-up study showed that scalarization occurs for β ≲ −4.35 Harada (1998).
In this work, we calculate the tidal Love numbers of neutron stars in scalar-tensor theories of gravity, focusing on the spontaneous scalarization case.Sec. 2 presents the equilibrium configuration for neutron stars in scalar-tensor theory.Sec. 3 discusses the first order linear time-independent perturbations upon which the tidal deformabilities depend.Sec. 4 details the method for deriving the various tidal Love numbers.Sec. 5 presents the results and demonstrates how the Love Numbers in scalar-tensor theories differ from those in GR.The paper concludes with Sec. 6, which discusses the results.

NEUTRON STARS IN SCALAR-TENSOR THEORY
Scalar-tensor theories are straightforward alternatives to GR.They depend on both a metric tensor (g µν ) and a massless scalar field (φ) and are typically expressed in one of two conformal frames: the Einstein frame and the Jordan frame.Historically, there has been much debate over the correct choice of frame Deruelle & Sasaki (2011), but it is now agreed that experiments measure Jordan frame quantities even though the field equations simplify in the Einstein frame Crisostomi et al. (2018); Pani & Berti (2014); Palenzuela et al. (2014); Doneva et al. (2013); Barausse et al. (2013).
In the Jordan frame, the action is where the tilde denotes Jordan frame quantities, ϕ is the Jordan frame scalar field, gµν is metric, R is the Ricci scalar, ω(ϕ) is a function of the scalar field that characterizes a specific scalar-tensor theory, and λ(ϕ) is the scalar potential.S m denotes the action of the matter, which is a function of the matter fields Ψ m and the Jordan metric gµν .Due to the ϕ R term, the gravitational constant G becomes a function of the scalar field i.e., G = G(ϕ).Throughout this work, we will continue to denote Jordan frame quantities with a tilde.The Jordan frame is the physical frame, but the field equations are typically expressed in the Einstein frame, where the metric and scalar decouple.A conformal transformation relates the two frames: (2) Using this transformation, the action can be re-written in a way that resembles the Einstein-Hilbert action: where all quantities are related to the Einstein metric g µν .φ is the Einstein frame scalar field, R is the scalar curvature, G * is the bare gravitational coupling constant which is set to 1, along with c, from here on.This paper will focus on the λ(φ) = 0 case.
The Jordan (ϕ) and Einstein (φ) frame scalar fields are related by the following equation Palenzuela et al. (2014): Much of the work presented here is applicable for any A(φ), but when necessary, the spontaneous scalarization coupling function Damour & Esposito-Farèse (1993) is used: (5) The modified field equations, derived from the Einstein frame action have the form where µν can be considered the stress-energy of the massless scalar field and has the form T µν is the stress-energy tensor in the Einstein frame, and T is the contracted stress-energy tensor T = T µ µ = g µν T µν .T µν is related to the Jordan frame stress-energy tensor ( Tµν ) in the following manner Note setting α(φ) to zero retrieves the GR field equations.We model neutron stars as static, spherically symmetric, non-rotating objects and assume that neutron star matter can be described as a perfect fluid.The stress-energy tensor for a perfect fluid is defined in the physical frame as where ũµ is the four-velocity of the fluid and ρ and p are the energy density and pressure in the Jordan frame.We assume that p and ρ are related by some barotropic equation of state so that where δ p and δ ρ are the Eulerian fluid perturbations.As the star is static, only the t component of the four-velocity is non-zero: Conservation of energy and momentum is defined in the physical or Jordan frame i.e., ∇µ T µ ν = 0. Transforming to the Einstein Frame gives The metric for a static, spherically symmetric, self-gravitating object is where ν and λ are functions of r and e −λ = 1 − 2µ(r)/r.The modified Tolman-Oppenheimer-Volkoff (TOV) or structure equations, which can be derived from the field equations and the equation for conservation of energy, have the form where µ is the mass function.ψ = ∂ r φ is used throughout this paper for improved readability.

STATIONARY PERTURBATIONS
In this section, we compute the linear, time-independent scalar and spacetime perturbations following the method initially laid out by Thorne & Campolattaro (1967).The complete system of time-dependent perturbations in scalar-tensor theory was calculated in Sotani & Kokkotas (2005), and the perturbation equations in this section have been cross checked with the extant results.
We use the Regge-Wheeler gauge Regge & Wheeler (1957), which separates the metric perturbation h µν into its even and odd parity components h µν = h + µν + h − µν .Sotani & Kokkotas (2005) demonstrated that the metric in both frames can be written in the Regge-Wheeler gauge using the proper redefinition of the metric components between frames.
For this analysis, as we are interested in time-independent perturbations, all perturbations (H 0 , H 2 , K, h 0 , and h 1 ) are functions of r only.Furthermore, the tr term (H 1 ) that is typically present in the Regge-Wheeler gauge vanishes.
The Einstein metric can be written in the following way: where and where Y ℓm (θ, ϕ) is the spherical harmonic function for l, m, and sym indicates that the metric is symmetric.
The explicit form of the conformal transformation between the Jordan and Einstein frame perturbation ( hµν → h µν ) is needed to determine the Jordan frame tidal deformability.The conformal transformation is obtained by perturbing Eq. (2) and substituting in the Regge-Wheeler metric Sotani & Kokkotas (2005).This gives where δA is the variation of the conformal factor; it is a function of the scalar field perturbation δφ = δφ(r).The relationship between δA and δφ depends on the functional form of the conformal factor.In the case of spontaneous scalarization δA = βA(φ)φδφ.The explicit relationships between the individual metric perturbations are We dropped the ℓm subscripts from H 0 , H 2 , K, h 0 , and h 1 for readability and will continue to do so throughout this work.
The complete set of perturbation equations needed to calculate the tidal deformability are laid out in Appendix B.
In GR, the full system of time-independent perturbed equations can be reduced to one differential equation for each parity: one for the even parity tensor perturbation (H) and one for the odd parity tensor perturbation (h).Scalar-tensor theories have an additional equation for the scalar field, which is of even parity.The metric and the scalar field are decoupled in the Einstein frame; therefore, the equations for H and δφ decouple.The Jordan frame perturbation H depends on both H and δφ.

NEUTRON STAR TIDAL DEFORMABILITY
We derive and compute the scalar-tensor tidal Love numbers and tidal deformabilities using the method developed by Hinderer (2008) and extended in Binnington & Poisson (2009); Damour & Nagar (2009).
Tidal deformabilities (e.g., λ) relate an applied external tidal field ( E ij...k ) to the induced multipole moment (Q ij...k ).To linear order in E ij...k , the tidal deformability is a proportionality constant between the two Hinderer (2008) i.e., Both E ij...k and Q ij...k can be decomposed into tensor harmonics where Y ℓm ij...k (θ, ϕ) are the even parity tensor spherical harmonics defined by Thorne (1980).This means that the tensor relation in Eq. ( 20) can be expressed as a scalar relation To calculate λ, it is sufficient to calculate one non-vanishing E m Hinderer (2008).
A scalar tidal deformability λ (φ) is defined analogously i.e. where ..k are the scalar tidal and multipole terms.The external tidal field and the induced multipole moment affect space-time in and around the neutron star.Outside the star, the large r behavior of the metric can be written in terms of E ij and Q ij Hinderer (2008); Thorne (1998).For example, the metric expansion for a spherically symmetric star of mass µ in a quadrupolar tidal field E ij for large r is where n i is the unit radial vector.

Electric Type Love Numbers
In GR, the electric or even parity Love numbers are calculated from the g tt = g (0) tt + h tt component of the metric and are based on a single, second order linear differential equation for H = H 0 = H 2 .However, in scalar-tensor theory, there are two types of even parity perturbations: scalar and tensor.The even parity metric tidal Love numbers k ℓ define how the body responds to a change in the metric.The scalar tidal Love numbers κ ℓ define how the body responds to a change in the applied scalar field.As the scalar field and the metric are not coupled in the Einstein frame, a change in the matter field does not induce a scalar perturbation, and vice-versa1 .The perturbation equations for the tidal Love number calculation must be derived carefully to first order in either the scalar perturbation or the metric perturbation but not both.This approach differs from previous approaches in Yazadjiev et al. (2018) and Pani & Berti (2014).
There are two master equations: one second order linear differential equation for the tensor perturbation H = H 0 = H 2 , which comes from the perturbation of the field equation (Eq.( 6a)), and one for the scalar perturbation δφ which comes from the scalar wave equation (Eq.( 6b)).
While the differential equation for φ can be derived directly from the scalar wave equation (see Eq. ( B35)), the differential equation for H is derived from the a system of equations Eqs.(B29) to (B34) and is obtained by the following steps (which have been widely used in GR Hinderer (2008); Damour & Nagar (2009); Binnington & Poisson (2009)): This gives where a prime (') denotes the derivative with respect to r and λ refers to the metric function and not the tidal deformability.In the case of spontaneous scalarization, Eq. ( 27) becomes External to the star, Eqs. ( 26) and ( 28), reduce to Eq. ( 29a) depends on the ψ and so is coupled to the scalar wave equation (Eq.( 6b)).As long as ψ > 0, there is no analytical solution to Eq. (29a).Only approximate solutions exist at the surface of the star (p = 0) because ψ ̸ = 0. Since φ asymptotically approaches a constant value φ ∞ , the derivative ψ vanishes at large r.In this regime, Eq. (29a) has an exact solution.The solutions to Eqs. (29a) and ( 29b) are where P m ℓ and Q m ℓ are the associated Legendre functions of the first and second kind.In order to determine c 1 , c 2 , d 1 , and d 2 , we match the asymptotic behavior of the two solutions i.e.
to the expansion of the g tt and the scalar field component of the metric (Eq.( 25)) respectively.This gives us c 1 , c 2 , d 1 , and d 2 in terms of the tensor tidal deformability λ and the scalar tidal deformability λ (φ) respectively.For example, in the ℓ = 2 case, we have By requiring continuity of the logarithmic derivatives and thus of H, δφ, and their derivatives at the surface of the star, it is possible to determine λ and λ (φ) in terms of µ, r, and either y or w respectively.This is done by substituting Eq. ( 32) and Eq.(30a) or Eq. ( 33) and Eq.(30b) into Eq.( 34) and solving for λ or λ (φ) .The tidal Love numbers are connected to the tidal deformabilities by the following equations Lastly, we can define the dimensionless tidal deformabilities: While this approach is sufficient to define λ, difficulty arises in numerically calculating λ and k ℓ because Eq. ( 30a) is only a solution to Eq. (29a) in the large r limit.It is not a solution near the surface of the star where numerical matching is typically done.In this region, an exact solution does not exist.While there is not an exact solution, an approximate series solution to Eq. ( 29a) can be constructed order by order in powers of r/µ.The leading order behavior of H is However, the leading order solution alone is not accurate enough for our purposes.The solution to Eq. ( 29a) is a linear superposition of the growing and diminishing solutions with two coefficients a ℓ and a −(ℓ+1) which are determined by the boundary conditions.
To create a more accurate solution, we construct two series solutions by adding higher-order terms.There is one growing and one diminishing solution that correspond to the two terms in Eq. (37).From there, higher-order terms are added to construct a solution with the form For numerical purposes, the series is truncated at order n = 13.This ensures the series has converged within 0.5%.
Note that H has only two degrees of freedom (a ℓ ,a −(ℓ+1) ).All other constants, α + and α − , are functions of these two.The constants are determined by substituting one series solution, either growing or decaying, into Eq.(29a) and solving for the coefficients order by order.

Magnetic Type Love numbers
The odd parity or magnetic Love numbers j l and their associated tidal deformabilities σ l are functions of the odd parity metric perturbation h − µν (h 0 , h 1 ).The odd parity metric perturbations h 0 and h 1 (Eq.( 17)) are coupled only to the explicit fluid velocity perturbation U (r)Y ℓm : The odd parity metric perturbations can be constrained by three equations, which come from the tϕ, rϕ, and θϕ components of the perturbation equations (see Eqs. (B36) to (B38)).There are multiple approaches to the magnetic Love number in the literature, but two are worthy of note Pani et al. (2018).The earliest two publications on magnetic tidal deformabilities Binnington & Poisson (2009) and Damour & Nagar (2009) have approaches that are fundamentally different and whose results do not agree.The first approach developed by Binnington and Poisson Binnington & Poisson (2009) assumes a strictly static fluid i.e., h 0t = h 1t = U = 0.The second approach from Damour and Nagar Damour & Nagar (2009) assumes an irrotational fluid.Instead of initially setting h 0t = h 1t = U = 0, this approach calculates the full Regge-Wheeler equation and then takes the static limit (ω → 0).Note that these approaches seem equivalent at a surface level but do not lead to the same answer because the irrotational approach picks up a non-vanishing term from the fluid velocity perturbation.This section lays out both approaches and clarifies the subtle differences between them.

Static Approach
In this section, we apply the static method derived in Binnington & Poisson (2009) to scalar-tensor theory; we assume that the perturbations are strictly static i.e., h 0t = h 1t = U = 0 and consider the odd parity perturbation equations Eqs.(B36) to (B38).Under this assumption, Eq. (B37) becomes and Eq.(B38) becomes independent of h 0 and constrains only h 1 .The final remaining equation, Eq. (B36), yields a second order differential equation for h 0 .
In the region exterior to the star, µ(r) = µ and p = ρ = 0 and Eq. ( 41) takes on a simpler form.
This equation differs from the GR equation by the factor of −2e −λ ψ 2 Binnington & Poisson (2009).Eq. ( 42) is now coupled to the scalar wave equation (Eq.( 6)) and no longer has an exact solution.This result differs from the f (R) results Yazadjiev et al. (2018).The coupling functions A(φ) differ between the two theories.Additionally, in f (R) theories, φ approaches zero as r → ∞, whereas in the theories considered here φ approaches a constant, non-zero value.Rather than matching solutions at the surface, Yazadjiev et al. (2018) matches the numerical solution to an analytical solution at some r match , beyond which the ψ term can be neglected.r match is defined by the Compton wavelength of the scalar field.
Since Eq. ( 42) is true for all r > r s , where r s is the surface of the star, including the large r regime where r ≫ r s and ψ → 0, there is an exact solution in the large r limit.This is sufficient to define the tidal Love numbers.
However, as was the case with the electric tidal deformability, there is no analytical solution at the surface of the star.Furthermore, the static approach is considered less physically relevant than the irrotational approach.
While an approximate solution to the static case could be constructed using the method laid out in Sec.4.1 or a method similar to that used in the f (R) case by Yazadjiev et al. (2018), this work focuses on the irrotational case because it is more realistic and has an analytical solution Pani et al. (2018); Landry & Poisson (2015); Shapiro (1996).

Irrotational Approach
In the irrotational approach, which was initially presented in Damour & Nagar (2009), it is assumed that the perturbations have a standard e −iωt time dependence i.e. h i (r, t) = h i (r)e −iωt .
Under this assumption, Eq. (B38) can be rewritten as where Ψ is defined such that Assuming that h 0 (r, t) = h 0 (r)e −iωt , Eq. ( 43) can be used to define h 0 : It is evident from this equation that h 0 is not well defined in the ω → 0 limit Pani et al. (2018).Substituting Eq. ( 45) and Eq. ( 44) into Eq.(B36) gives the following master equation: This agrees with Equation 40 in Sotani & Kokkotas (2005).
Since it is assumed that the neutron star is static, we are interested in the ω → 0 limit.The master equation becomes Outside of the star, this equation simplifies further: Interestingly, this equation, unlike the static master equation (Eq.( 41)), does not depend explicitly in φ or ψ.Therefore external to the star, the solution to Eq. ( 47) is known and identical to the GR solution.All non-GR effects arise from matching the internal and external solutions at the surface of the star.
We briefly demonstrate the difference in the static and irrotational solutions in scalar-tenor theory.The method presented in Pani et al. ( 2018) is applied to the scalar-tensor problem.
Using the axial component of the stress-energy tensor conservation equation (Eq.( 12)) and assuming ω ̸ = 0, one finds that Substituting Eq. ( 49) into Eq.( B36) and then taking the ω → 0 limit, the following differential equation for h 0 is obtained There is a sign change in the 8πA 4 (φ)(p + ρ) term between Eq. ( 41) and Eq. ( 50).The difference occurs because U (r) = 0 for the static approach and U (r) ̸ = 0 in the irrotational case.So while there is irrotational fluid motion in one case, the other has a completely static fluid.As a result of this difference, the irrotational tidal Love numbers are negative while the static Love numbers are positive.
Returning to the main goal of this work, calculating j ℓ , σ ℓ , and Σ ℓ we use the similarity between the irrotational master equation and its GR counterpart to define the ℓ = 2 solution.Eq. ( 48) has an exact solution of the form where R = r/µ and F is a hypergeometric function.For ℓ = 2, F is expressible in terms of simple functions.b q and b p are determined by the boundary conditions at the surface of the star.Since both ψ and ψ ′ are required to be continuous at the surface of the star, the logarithmic derivative y odd = rψ ′ /ψ must also be continuous at the surface of the star.j ℓ , σ ℓ , and Σ ℓ are therefore defined to be and where R s = r s /µ and C is the compactness.j ℓ and σ ℓ are, at a glance, identical to their GR counterparts, but the scalar-tensor and GR values differ because all non-GR effects are contained in the value of y odd s calculated by integrating Eq. ( 47) along with the modified TOV equations (Eq.( 14)) inside the star.

Electric Love Numbers
This section presents the electric tidal Love numbers and the associated tidal deformabilites and compares them to the GR results.There are two degrees of freedom needed to define a specific case of spontaneous scalarization: β and φ ∞ .β is constrained by binary pulsar experiments to β < −5 at the 1σ level Freire et al. (2012).In this work we use several values of β to demonstrate the results: β = −4.5, −5, −5.5, −6.Generally, the figures compare only the β = −4.5 and β = −6.This gives two sets of results, one conservative and one optimistic.The value of the scalar field at infinity, φ ∞ , is tightly constrained by the Cassini experiment Bertotti et al. (2003).That experiment directly constrains the Brans-Dicke parameter ω BD to be > 4 × 10 4 .The value of the scalar field at infinity is related to the Brans-Dicke parameter by the equation This constrains φ ∞ to < 2.7 × 10 −3 and < 2.0 × 10 −3 for β = −4.5 and β = −6 respectively.We use φ ∞ = 10 −3 for all results presented.Changing φ ∞ to 2.0 × 10 −3 increases the deviation from GR. Conversely, changing φ ∞ to 10 −4 decreases the deviation from GR.These differences grow with increasing compactness but are less than 1% for the values considered.Using Eq. (30a) and Eq. ( 35), it is possible to define the tidal Love number in the large r limit.The scalar tidal Love number can be similarly calculated.
The ℓ = 2 tidal Love numbers are defined as follows where y = rH ′ /H and w = rδφ ′ /δφ.y is traditionally evaluated at the star's surface for numerical applications.However, Eq. ( 56a) is not valid when r = r s .Close to the star ψ ̸ = 0, and the solution to Eq. (29a) can only be approximated.After constructing a series solution that is accurate to better than 0.5% for even the largest values of ψ considered (see Sec. 4.1), we compared values from the exact solution (Eq. ( 56a)) evaluated at the surface to the values of k 2 calculated from the approximate solution.The tidal deformabilities agreed to better than 3.7% for all equations of state and values of β explored.The percent difference between the approximate and the exact values is strongly dependent on the compactness and increases with increasing compactness.For the vast majority of the parameter space explored, the difference between scalar-tensor theory and GR is larger than the difference between approximate and exact solutions.Exceptions occur for β = −6 where the scalar-tensor theory and GR curves intersect.This can be seen in Fig. 1.
Fig. 1 shows how the electric tidal Love numbers and tidal deformabilities differ in scalar-tensor theory and GR.Three different equations of state are considered: FPS, SLy, and MS1.These equations of state cover a wide range of stiffness and support a maximum mass of > 1.8M ⊙ .FPS and SLy are both within constraints from analyses of GW170817 Capano et al. (2020); Abbott et al. (2017).However, as NICER results favor stiffer equations of state, we include MS1 Raaijmakers et al. (2021Raaijmakers et al. ( , 2020)); Bogdanov et al. (2019a,b).
Fig. 1 plots the physical or Jordan frame values, which are related to their Einstein frame counterparts by Eqs. ( C56) and (C57).In Fig. 1a, Fig. 1b, and Fig. 1c, the observables λ2 , Λ2 , and k2 are plotted against the neutron star's compactness ( C), in this case defined as the Jordan frame TOV mass ( M ) over the Jordan frame radius rs .In Fig. 1d, Fig. 1e, and Fig. 1f, the percent difference between scalar-tensor theory and GR is shown, also as a function of compactness.Note that Fig. 1e and Fig. 1f are essentially identical.This is due to the definition of the dimensionless tidal deformability (Eq.( 36)).As the tidal Love number and the dimensionless tidal deformability are related by a factor of 3 2 C5 and C is the x-axis variable, the factors of 3 2 C5 cancel out.This can easily be shown by substituting the definition of the tidal Love number into the equation for the percent difference and forcing C GR = C.This same phenomenon appears in Fig. 3.
It is clear that the spontaneous scalarization effect can lead to significant deviations from the GR tidal deformabilities.It is also clear that the deviations are strongly dependent on the objects compactness and the coupling constant.For the case where β = −6, the tidal Love number and dimensionless tidal deformability differ at most by ∼ 25% and the tidal deformability differs by up to ∼ 200%.The peak occurs around C ≈ 0.25 for the tidal deformability and ≈ 1.9 for the Love number, with the exact value varying by equation of state.In the more conservative case where β = −4.5, this reduces to ∼ 15% for the tidal Love number and ∼ 20% for the tidal deformability, and the peaks occur around C ≈ 0.3 and ≈ 2.3 respectively.
The tidal deformability curve for scalar-tensor theories has a different shape than those in GR: a second peak appears.This peak is small for the weak coupling case, but for more negative coupling constants, the second peak is clear.This second peak is caused by the spontaneous scalarization effect, which causes large deviations from GR in conditions with strong gravitational fields Damour & Esposito-Farèse (1993, 1996).
As the difference between scalar-tensor theory and GR is much greater than the difference between the two methods of calculating k ℓ , we consider Eq. (56a) evaluated at the surface of the star to be sufficiently accurate for GW parameter estimation with current detectors.
We show the Jordan frame ℓ = 2 scalar tidal Love numbers and tidal deformabilities in Fig. 2. Scalar tidal Love numbers will effect scalar GW emission Bernard (2020).We find that scalar tidal deformabilities are much smaller than the electric tidal deformabilities, around two orders of magnitudes smaller even for strongly scalarized cases.Additionally, the scalar tidal deformabilites and tidal Love numbers depend strongly on the coupling constant, with strong scalarization leading to negative scalar tidal love numbers.

Magnetic Love Numbers
This section presents the magnetic tidal Love number and the associated tidal deformabilites in scalar-tensor theory and compare them to the GR results.
The exact equations for the magnetic tidal Love numbers j ℓ and tidal deformabilities σ ℓ can be determined by substituting Eq. (51) into Eq.( 52) and Eq. ( 53).
The explicit equation for the ℓ = 2 or quadrupolar tidal Love number is where C = µ/r s is the Einstein frame compactness and y = y odd (r s ) = r s Ψ ′ /Ψ is the logarithmic derivative a the surface.Fig. 3 shows the Jordan frame l = 2 love numbers, tidal deformabilities, and the difference between the GR and scalar-tensor tidal effects.The Jordan frame values are related to their Einstein frame counterparts by Eq. (C48a) and Eq.(C48b).The conformal transformations are derived in Appendix C.
It is clear that tidal Love numbers and tidal deformabilities differ between GR and scalar-tensor theory.For the optimistic case where β = −6, the tidal Love number has a maximum deviation of ∼ 17% and the tidal deformability has a maximum deviation of ∼ 300%.This maximum deviation occurs at C ≈ 0.2 for the Love number and ≈ 0.24 for the tidal deformability.In the more conservative case where β = −4.5, the peak occurs at C ≈ 0.23 for all tidal properties and the deviation changes to ∼ 2.5% and ∼ 15% for j 2 and σ 2 respectively.
In GR empirical relationships between the ℓ = 2 dimensionless magnetic and electric tidal deformabilities have been found Forteza et al. (2018).The dimensionless magnetic tidal deformability Σ 2 and the dimensionless electric tidal deformability Λ 2 have a quasi equation of state independent relationship: We find that the scalar-tensor tidal deformabilities can be fit to a similar relationship, with the coefficients depending on the value of β.Regardless of equation of state, R 2 > 0.99 for all cases.There does not appear to be a similar relationship between the scalar deformability and the electric tidal deformability.The scalar tidal deformability depends strongly on β and Λ (φ) /Λ can take on different shapes that are equation of state dependent.

DISCUSSION
This work presents the electric, magnetic, and scalar tidal Love numbers and tidal deformabilities.We find that the electric and magnetic tidal effects may differ significantly from their general relativistic counterparts (∼ 200 and ∼ 300 for electric and mag- netic respectively).These large deviations occur at larger compactnesses (≳ 2) and are caused by the spontaneous scalarization effect.The exact deviation and the compactness where this maximum deviation occurs are equation of state dependent.
This paper approaches tidal effects through the lens of GW parameter estimation.The mass-radius-tidal deformability relationships explored in this paper can be applied directly to GW parameter estimation of GWs from binary neutron star and neutron star-black hole systems.The ℓ = 2 dimensionless electric tidal deformability is the leading order tidal effect for GWs.Given that this number can vary by ∼ 25 between scalar-tensor theory and GR, it may be necessary to take modified tidal effects into account when doing tests of GR using GWs from systems with neutron stars.
We present an analytical expression for the magnetic tidal Love numbers in scalar-tensor theory for the first time.The results establish that the magnetic Love numbers are only implicitly dependent on the scalar field and have an analytical solution.This is in agreement with Sotani & Kokkotas (2005), which shows that the time-dependent perturbation equation is only implicitly dependent on the scalar field.However, this was discussed only in the context of perturbations and not of tidal Love numbers.
The magnetic Love numbers in this paper can be compared to their f (R) counterparts because f (R) theory and scalar-tensor theory are mathematically similar.Differences arise when calculating the tidal Love numbers in part due to the behavior of the scalar field at infinity.In f (R) theory, the scalar field and its derivative go to zero at infinity, and the tidal deformability can be evaluated at some distance away from the neutron star where both the scalar field and its derivative are sufficiently small.This is different from the scalar-tensor theories considered in this paper where the scalar field asymptotically approaches a constant.Additionally, the coupling function differs between theories with A(φ) ∝ e αφ in f (R).Despite this, the perturbation equations inside the star should agree when the correct substitutions for A(φ) and α(φ) have been made because they are mathematically similar.However, our perturbation equation differs from Eq.( 23) in Yazadjiev et al. (2018) by a negative sign.Comparing the tidal Love numbers themselves, shown in Fig. 3, with the results from Yazadjiev et al. (2018), it is clear that the qualitative features are consistent, with the deviation from GR increasing with compactness.However, the difference between GR and scalar-tensor theory are smaller than those between GR and f (R), at least for physically allowed values of β and φ ∞ .
This paper also includes the even parity tidal Love numbers and tidal deformabilities.The ℓ = 2 electric tidal Love numbers in scalar-tensor theory were initially presented in Pani & Berti (2014) in the context of the so-called "I-Love-Q" relations.The methods in this paper differ significantly from those in Pani & Berti (2014).
To begin, Pani & Berti (2014) use the most general stationary axisymmetric metric that includes first order rotation terms rather than the stationary Schwarzshield metric used in this paper.In addition to this, there is a fundamental difference between the definitions of the tidal Love numbers and tidal deformabilities between this work and theirs.This based on the way that the even parity perturbation equations are treated.There are both metric and scalar perturbations in the even parity case, and the relationship between them is not trivial.In the Einstein Frame, the metric tensor and the scalar field are not coupled.As a change in the metric should not affect the scalar-field and vice-versa, it is important to construct two independent first-order perturbation equations.One for the metric perturbation and one for the scalar.This differs from the approach in Pani & Berti (2014), where the two even parity equations are coupled.It is unsurprising, then, that Eqs. ( 26) and ( 27) are different from the equations presented in Pani & Berti (2014).The resulting tidal Love numbers must also differ.The definition for the scalar tidal Love number in this paper also differs from that in Pani & Berti (2014).Pani & Berti (2014) does not include a source term in their definition of the scalar tidal Love number, and we do.This is because they are considering the perturbation in the scalar field produced by a change in the metric rather than by a change in the scalar field.This paper also includes the ℓ = 3, 4 even parity Love numbers and tidal deformabilities in Sec.D, which have not been presented before.
The results demonstrate that tidal Love numbers and tidal deformabilities can differ significantly between scalar-tensor theory and GR.This is consistent with other results in the literature, which show that tidal Love numbers in f (R) theory and scalar-Gauss-Bonnet gravity Yazadjiev et al. (2018); Saffer & Yagi (2021) also differ significantly from their general relativistic counterparts.As GWs emitted by neutron stars depend on the tidal deformability, it is essential to take the changes in the mass, radius, and tidal deformability into account when studying GWs from neutron stars in theories beyond GR.The allowed deviations from GR in the GWs are smaller than or similar to the uncertainty in the tidal deformability measurements from binary neutron star and neutron star black hole mergers.By taking the modified tidal deformability into account, the small deviations from GR in the waveform can be more accurately determined.
The equation for the scalar perturbation δφ is derived by perturbing scalar wave equation Eq. ( B35)).

B.3. Equations for Odd Parity
The static and irrotational methods used in this paper differ in their treatment of time derivatives.Even though the tidal Love numbers themselves are time-independent, we present the time-dependent equations in this section.
Combining the δG µν with the matter stress-energy tensor and scalar stress-energy tensor terms results in the following three equations: • Equation B36 is δG tϕ = 8πδT tϕ + T The tidal Love numbers in this paper were derived in the Einstein frame; however, as experiments measure Jordan frame quantities, it is necessary to obtain the Jordan frame quantities using a conformal transformation.We assume here that the Jordan frame metric gµν is related to the Einstein frame metric g µν by a conformal factor A(φ): where A(φ) = e 1 2 βφ 2 . By construction, the Einstein frame metric is asymptotically flat.This implies that where η µν = diag(−1, 1, 1, 1) is the Minkowski metric.As the Jordan frame metric is also asymptotically flat or Minkowskian, the r and t components must be related to their Einstein frame counterparts in the following way: r = A(φ)r and t = A(φ)t.Furthermore, the effective gravitational constant G is no longer a constant in the Jordan frame and is not necessarily equal to the bare gravitational constant G which appears in the Einstein frame equations.The relationship between the two is known Palenzuela et al. (2014): We need the conformal transformations for the perturbations between the two frames to transform the tidal Love numbers and tidal deformabilities from the Einstein frame to the Jordan frame.These are presented in Sec. 3. it is possible to define the Jordan frame tidal deformabilty λJ : where E denotes Einstein frame tensor quantities and φ denotes Einstein frame scalar quantities.From this equation, it is clear that the Jordan frame tidal deformability is related linearly to the even parity scalar and tensor tidal deformabilities.The exact relationship is λJ = A 2ℓ+1 (φ ∞ )(λ E + λ (φ) ) . .

Figure 1 .
Figure 1.Panel (a) shows the ℓ = 2 Jordan frame tidal deformability in cgs units, (b) shows the ℓ = 2 Jordan frame dimensionless tidal deformability, (c) shows the ℓ = 2 tidal Love number, (d) shows the percent difference between the cgs tidal deformability in scalar-tensor theory and GR, (e) shows the percent difference between the dimensionless tidal deformability in scalar-tensor theory and GR, and (f) shows the percent difference between the tidal Love numbers in scalar-tensor theory and GR.All are shown as a function of the Jordan frame compactness.The value of the scalar field at infinity φ∞ for all cases presented here is 10 −3 .Three realistic nuclear equations of state (SLy, FPS, and MS1) are shown in black, blue, and red respectively.The results for β = −6 and β = −4.5 are shown with dash-dot and dotted line styles.The purple lines in (d), (e), and (f) indicate the percent difference between the analytical and approximate approaches to calculating the tidal Love number and tidal deformability.

Figure 2 .
Figure 2. Panel (a) shows the ℓ = 2 Jordan frame scalar tidal deformability in cgs units, (b) shows the ℓ = 2 Jordan frame dimensionless scalar tidal deformability, (c) shows the scalar tidal Love number.All are shown as a function of the Jordan frame compactness.The value of the scalar field at infinity φ∞ for all cases presented here is 10 −3 .Three realistic nuclear equations of state (SLy , FPS, and MS1) are shown in black, red, and blue respectively.The results for β = −6 and β = −4.5 are shown with dash-dot and dotted line styles.

Figure 3 .
Figure 3. Panel (a) shows the ℓ = 2 Jordan frame tidal deformability in cgs units, (b) shows the ℓ = 2 Jordan frame dimensionless tidal deformability, (c) shows the tidal Love number, (d) shows the percent difference between the cgs tidal deformability in scalar-tensor theory and GR, (e) shows the percent difference between the dimensionless tidal deformability in scalar-tensor theory and GR, and (f) shows the percent difference between the tidal Love numbers in scalar-tensor theory and GR.All are shown as a function of the Jordan frame compactness.The value of the scalar field at infinity φ∞ for all cases presented here is 10 −3 .Three realistic nuclear equations of state (SLy, FPS, and MS1) are shown in black, red, and blue respectively.The results for β = −6 and β = −4.5 are shown with dash-dot and dotted line styles.

Figure 4 .
Figure 4. Quasiuniversal relations between the Jordan Frame dimensionless magnetic quadrupolar tidal deformability (Σ2) and the Jordan Frame dimensionless electric quadrupolar tidal deformability Λ2.Only the irrotational magnetic tidal deformability is shown in this figure.Three equations of state (SLy, FPS, and MS1) are shown in different line styles, but they are indistinguishable.The color of the lines corresponds to different values of β.