Multipole moments and universal relations for scalarized neutron stars

In recent years there has been a surge of interest in what has come to be known as the `universal relations' between various global properties of neutron stars. These universal relations are equation of state independent relations between quantities such as the moment of inertia $I$, the tidal deformability or Love number $\lambda$, and the quadrupole $Q$ (I-Love-Q relations), or the relativistic multipole moments (3-hair relations). While I-Love-Q relations have been studied extensively in both general relativity and various alternatives, 3-hair relations have been studied only in general relativity. Recent progress on the definition of the multipole moments of a compact object in the case of scalar-tensor theories allows for the study of 3-hair relations in modified theories of gravity. Specifically, the aim of this work is to study them for scalarized stars in scalar-tensor theories with a massless scalar field that admit spontaneous scalarization. We find that the 3-hair relations between the mass and angular momentum moments that hold in general relativity hold for scalarized stars as well. The scalar moments also exhibit a universal behaviour, which is equation of state independent within one specific theory, but differs between different theories. Combining astrophysical observations one can in principle measure the different properties of scalarized neutron star and tell different theories apart.

I. Introduction 1 II. Stars in scalar-tensor theory 2 III. Mass, angular momentum and scalar field moments in scalar-tensor theory 4 IV. EoS independent behaviour of scalarized stars. 5 V. Relating moments to observables and comparison to GR. 9 A. Observables and moments 9 B. Measuring the scalar charge and β. 10 1. Setting up the problem and constraints. 10 2. Using universal relations to overcome degeneracies. 10 VI. Conclusions 11 Acknowledgments 12 A. Calculation of the multipole moments 12

I. INTRODUCTION
Neutron stars have been studied extensively in scalar-tensor theories (STT) of gravity. In contrast to black holes that do not possess scalar hair in these theories [1][2][3], for neutron stars the matter acts as a source of the scalar field and supports nontrivial scalar configurations. The simplest class of STT, Brans-Dicke theory, always leads to the development of a nontrivial scalar field for all matter configurations and the differences with general relativity (GR) can be considerable. This is actually a disadvantage, as Brans-Dicke theory deviates from Einstein's theory of gravity in the weak field regime where GR is tested with very high precision. Hence, weak-field test can already place tight constraints that leave little room for strong-field deviations. This argument can be circumvented in a specific class of theories that are perturbatively equivalent to GR in the weak field regime but exhibit significant deviations in the strong field regime [4]. In particular, neutron stars in such theories can exhibit spontaneous scalarization. 1 The essence of spontaneous scalarization is that, once a star exceeds a compactness threshold, it is energetically more favorable to develop a nontrivial scalar configuration [4,9]. Scalarized neutron stars have been examined further [10,11] including slow [12][13][14][15] and rapid [16,17] rotation.
As in GR, the properties of neutron stars in STT depend on the specific equation of state (EOS) that one selects to de-scribe their interior. The various uncertainties in the microphysics result in a proliferation of EOSs. The properties of the star and the relevant observables depend both on the choice of the EOS and on the theory of gravity and these dependencies appear to be degenerate. That is, using a different theory of gravity, such as STT, is hard to distinguish from changing the EOS. A way around this problem has already been proposed and continues to be developed in the form of universal, or EOS independent, relations between various quantities that characterise the structure of neutron stars (for some recent reviews see [18,19]). Such universal relations have been studied in STT as well, though not as widely as in GR. A particular class of the universal relations that has been studied in GR but not in STT are the 3-hair relations between the multipole moments of a neutron star (for the results in GR see [20][21][22]). The goal of this work is to extend the GR results on 3-hair relations to the case of scalarized stars and to explore the possibility of having a universal description of the moments of the additional degree of freedom, i.e., the scalar field. We will also explore how the scalar moment are related to observables and whether we can use such observables to put constraints on the parameters of the STT, as well as the degree of scalarization of a neutron star.
Since the class of STT discussed above cannot be constrained by the weak field experiments, one has to use observations involving strong field effects, such as the gravitational wave emission of neutron stars located in close binary systems leading to shrinking of their orbits. These observations pose strong constraints on the theory [12] and the latest results lead to tight bounds on the free coupling parameters [23,24]. Thus, the nonrotating and slowly rotating scalarized neutron stars do not differ significantly from the GR case and it would be very difficult to probe the presence of a scalar field. Only the rapidly rotating case can lead to larger deviations from the non-scalarized solutions [16]. This motivates an investigation of the universal relations for rapidly rotating stars, which is the focus of this paper.
We should note that there is another way to evade the strong constraints coming from the binary pulsar observations, namely the inclusion of nonzero scalar field mass. This would lead to a finite range for the scalar field of the order of its Compton wavelength and can reconcile the theory with observations for a much larger range of parameters. Nonrotating neutron stars in massive STT were examined for the first time in [25,26] and the results were extended in [27] and [28] for slow and rapid rotation respectively. The studies indeed showed that the neutron stars can differ dramatically from the pure general relativistic case. Defining the multipole moments in these theories is more complicated though and we will leave it for future studies.
The rest of this work is organised in the following way, Section II gives a brief description of the formalism for construction neutrons stars in STT, while Section III gives a brief description of the calculation of the moments. Section IV presents the results on the various universal relations between the multipole moments, and Section V discusses how the various relations could be used to extract information about the moments and the particular STT from observations. Finally, we end with our conclusions.

II. STARS IN SCALAR-TENSOR THEORY
The general form of the Einstein frame action in STTs with a massless scalar field is [29][30][31] where G * is the bare gravitational constant, R is the Ricci scalar curvature with respect to the Einstein frame metric g µν , the matter fields are collectively denoted by Ψ m and their action is S m . In the Einstein frame the scalar field ϕ is directly coupled to the matter via the function A(ϕ). This function plays the role of a conformal factor that relates the Einstein frame metric g µν to the Jordan frame metric g µν = A 2 (ϕ)g µν . By definition, the matter fields couple minimally to the Jordan frame metric and this guarantees that the weak equivalence principle is satisfied. We have chosen to work in the Einstein frame, as in this frame the field equations have the same structure as in GR and this simplifies calculations. Moreover, the multipole moments presented below have been previously defined and calculated in the Einstein frame [32]. We stress that any physical quantities in the Jordan frame can be expressed in terms of these moments [33].
In what follows we use geometrical units c = G * = 1 and the dimensional quantities are given in km. We will focus on stellar configurations that are stationary and axisymmetric and we will describe the matter in the Einstein frame as a perfect fluid with pressure p and energy density ε. The spacetime metric can then be written in the following general form All metric functions γ, σ, ω and α, as well ϕ, p and ε, will depend only on the coordinates r and θ. For numerical calculations it is more convenient to use the angular coordinate µ = cos θ instead of θ. Using our ansätze, the field equations that one obtains from varying the action with respect to the metric take the form ∆(σe γ/2 ) = e γ/2 8π(ε + p)e 2α 1 + υ 2 where the differential operator ∆ is defined as The last field equation (6) for the metric function α is of first order compared to the second order field equations for the rest of the metric functions. The field equation for the scalar field is The above system of equations has to be supplemented with equations that describe the dynamics of the fluid, namely the equation for hydrostatic equilibrium and the equation of state (EOS) for nuclear matter. The latter is a relation between pressure and energy density and we impose it in the Jordan frame, since the matter couples minimally to the Jordan frame metric. This minimal coupling also implies that the fluid will satisfy the usual conservation laws in terms of the Jordan frame variable,p andε. The equations above have been given in the Einstein frame and p and ε are related top andε as follows One can use these relations to express the EOS and the equation for hydrostatic equilibrium in terms of p and ε in order to work exclusively with Einstein frame variables. We find it more convenient to work directly withp andε. Hence, in the numerical implementation we use eqs. (9) to express p and ε in termsp andε in all the equations above and we express the equation for hydrostatic equilibrium in the form where we have introduced the coupling function k(ϕ) = d ln(A(ϕ))/dϕ. The Einstein frame four velocity u µ is defined as where the proper velocity is v = (Ω − ω)r sin θe −σ and Ω is the fluid angular velocity Ω = u φ /u t .
What is left to be fixed then is the particular form of the Einstein frame coupling function. We will work with the standard choice k(ϕ) = βϕ where β is a constant. One of the most important properties of this class of scalar-tensor theories is that it is perturbatively equivalent to GR in the weak field regime, while in the strong field regime nonlinear effects lead to nonuniqueness of solutions and spontaneous scalarization [4]. In the calculations below we will allow also for nonzero background value of the scalar field ϕ ∞ in some cases.
We solve the field equations using a modification of the RNS code (see [34] for the original GR version of the RNS code while the STT extension can be found in [16]). This code is based on the KEH method [35] with certain modifications introduced in [36]. A key property of this method is that the field equations are presented in an integral form. This turns out to be very useful for the calculation of the multipole moments, as explained in Appendix A.

III. MASS, ANGULAR MOMENTUM AND SCALAR FIELD MOMENTS IN SCALAR-TENSOR THEORY
Here we give a brief description of the framework and the general results for the moments in the Einstein frame [32] for a STT with a massless scalar field. More details on the particular calculation of the moments employed in the RNS code can be found in Appendix A, while a general review of the calculation in GR can be found in [19].
When discussing the multipole moments it is more convenient to use the following form of the metric that is written again in quasi-isotropic coordinates similar to the metric used by the RNS code (2), but with the new functions B = e γ and ν = (γ + σ)/2, i.e., ds 2 = −e 2ν dt 2 + r 2 sin 2 θB 2 e −2ν (dϕ − ωdt) 2 +e 2α (dr 2 + r 2 dθ 2 ). (12) The field equations for this metric are directly related to the ones given in the previous section, i.e., eqs. (3)- (6). Note that the Einstein frame field equations (3)-(5) are identical to their GR counterparts 2 (given in [36]), while eq. (6) and the equation for hydrostationary equilibrium (10) have some additional contributions involving derivatives of the scalar field ∂ i ϕ. Therefore, as discussed in more detail in [16,32], in the case of a massless scalar field the multipole moments can be calculated in the same way as in GR with scalar field corrections entering through the Ricci tensor and the equation for α, i.e., eq. (6). Similarly, the vacuum field equations for the metric functions ν, B and ω, which are used to define the moments, are the same as in GR and can be found in [37]. One can easily show that the asymptotic expansion of the metric functions and the scalar field admits the following ansatz in terms of the Legendre Polynomials P l (µ), their derivatives dP l (µ) dµ , and the Gegenbauer polynomials T where the coefficients in these expansions are of the form The calculated multipole moments of the spacetime are combinations of the expansion coefficient in (17)-(20) (see discussion for the GR case in [19]) and below we give explicitly the first few multipole moments using the formalism developed in [32]. We should note that even though the calculation of the metric coefficients (17)- (19) is the same as in GR, the coefficients in the scalar field expansion enter explicitly in the multipole moments given below.
Mass (monopole): Scalar monopole: Angular momentum (dipole): Mass quadrupole: Scalar quadrupole: Spin octupole: Mass hexadecapole: Scalar hexadecapole: Spin 2 5 -pole: These are all the non-zero multipole moments up to S 5 for a stationary and axisymmetric spacetime with equatorial symmetry and in the presence of a scalar field with the same symmetries.
As emphasized earlier already, these moments are the Einstein frame moments. Defining the moments in the Einstein frame is straightforward, while attempting to do so in the Jordan frame appears to be significantly harder. The Jordan frame is related to the Einstein frame through a conformal transformation that depends on the scalar field. In terms of the multipole moments, a conformal transformation of the metric would generally result in a mixing of the moments, with the new moments being combinations of the old ones. This is clearly not an essential redefinition of the multipole moments. In our specific case we would additionally have the mixing of mass and angular momentum moments with scalar field moments, due to the conformal factor being a function of ϕ. Any physical quantity that one would wish to express in term of some Jordan frame moments (assuming that they can be rigorously defined), can be always reexpressed in terms of the Einstein frame moments, using the relations between Einstein and Jordan frame variables. A further advantage of the Einstein frame moments is the following. In the context of STT the functional form of the conformal factor A 2 (ϕ) is specific to a theory or a class of theories and can be parameterised in terms of appropriate parameters or coupling coefficients of the theory. In the selected formulation, these coupling coefficients of a specific theory appear as the coefficients that mix the Einstein frame moments, instead of being hidden in some Jordan frame moments. This gives a more transparent handle on a specific STT (see for example [33]). Therefore, while the choice made here does not lose in generality, it can be further argued to be multiply advantageous.
As a last note, we mention that below we will use the reduced moments defined as, where n ≥ 0, i is the imaginary unit, and j ≡ J/M 2 . A similar normalisation will be used for the scalar moments, but this will be further explained in the following section.

IV. EOS INDEPENDENT BEHAVIOUR OF SCALARIZED STARS.
To explore the existence of universal relations between the various moments of scalarized stars, similar to the 3-hair relations between the moments in GR, we have constructed sequences of scalarized models using various EOSs. For these stars we have calculated the mass and angular momentum moments up to M 4 , as well as the scalar moments up to W 4 , following the procedure outlined in the previous section and the expressions given there. While the mass and angular momentum moments can be directly compared to their GR counterparts, the scalar moments don't have a GR counterpart and are in this sense novel features.
We use several equations of state in order to cover a wide range of stiffness. These are the APR4 [38], SLy4 [39], A [40], FPS [41] and the zero temperature limit of the Shen EOS [42,43]. APR4 and Sly4 are modern realistic EOS that are in agreement with the observations. EOS A and FPS are too soft and already excluded by observations, as they do not reach two solar masses [24,44]. The Shen EOS does reach the two solar mass barrier, but it is stiffer and leads to somewhat larger radii, so it is disfavoured by observations [45][46][47]. We have included softer and stiffer EOS even though they are ruled out or disfavoured, as our main goal is to demonstrate the universality of the relations given below. Hence, it is instructive to use a broader set of EOS in order to verify that this universality is not simply a residual effect from considering EOSs with very similar properties.
The scalarized models have been constructed assuming values of β in the range between −4 and −10 covering a big part of the parameter space. We should note that the current observational limit is β > −4.5 [24,44] for theories with a massless scalar field. Nevertheless, we have again decided to include larger values of |β|, to demonstrate that the universality persists for significantly scalarized stars and it is not an artefact of very weak scalarization. It is worht mentioning that considering values of β lower than −4.5 is justified if one includes a mass for the scalar field. In that case, the scalar field is confined within its Compton wavelength and therefore, for large enough scalar field masses, the emission of scalar gravitational radiation is suppressed and binary pulsar observations cannot set as tight constraints on the parameter β [25][26][27] as in the massless case. One should note however that defining the multipole moments in the case of massive scalar field is much more involved because of the finite range of the scalar field and its exponential decay at infinity. This remains an open problem which we plan on addressing in future work. One more issue we should address at this point is that of the asymptotic value of the scalar field, which in the class of models that we are investigating, is constrained to be almost zero by observations. Nevertheless, we have also calculated models with a non-zero asymptotic value of the scalar field ϕ 0 in order to have a more complete investigation of scalarized stars. For these latter models we have used a somewhat larger value of ϕ 0 (i.e., ϕ 0 = 0.03), similar to previous studies [48], in order to have a better assessment of how that would affect the behaviour of the universal relations.
For the particular choice of the coupling function the field equations are invariant under the transformation ϕ → −ϕ. Thus the neutron star solutions with opposite signs of the scalar field are otherwise indistinguishable (e.g. the metric functions describing the two solutions are the same). Therefore, in the presented results we have chosen arbitrarily one particular sign of ϕ and normalised the scalar field multipole moments accordingly. In any case, solutions with the opposite sign for the scalar field also exist and would simply result to scalar moments with an opposite sign.
We now proceed with the presentation of our results. Our first results concern the mass and angular momentum moments of scalarized stars and their behaviour with respect to their GR counterparts. These results are shown in Fig. 1, where we have plottedS 3 againstM 2 andM 4 againstM 2 . The figures include models within the full range of the β parameter that we have used both zero and non-zero asymptotic values of the scalar field.
As shown in Fig. 1, theS 3 −M 2 andM 4 −M 2 relations of GR [20,22] hold for scalarized stars as well. It would be useful at this point to contemplate on this very interest- ing result. Considering the Einstein frame multipole moments as GR moments with some additional corrections due to the scalar field, one might be tempted to conclude that this results is expected. Indeed, in certain cases it has been argued [48] that the main effect of the scalar field is to stiffen the Einstein frame EOS with respect to the prescribed Jordan frame EOS [c.f. eqs. (9)]. However, the presence of the scalar field is not in general trivially equivalent to an EOS change, since the gradient of the scalar field itself also acts as a source in the field equations. Furthermore recent studies on I − Q relations for scalarized stars have shown, in contrast to what we find here for the 3-hair relations, that for large values of |β| the scalarized I − Q relations can somewhat deviate from the corresponding GR I − Q relations [14,49]. Therefore, what we find here for the 3-hair relations and for values of β as much as β = −10 is quite intriguing. Overall it seems that the 3hair universal relations are quite less sensitive to the choice of β than the I − Q universal relations, being in a sense more universal with respect to different theories of gravity.
The results presented here though, do not eliminate the possibility that stars with an extremely high degree of scalarization in the context of STT or in the context of exotic object in GR (e.g. a mixed boson-neutron star) could deviate from these relations. Such objects are beyond the scope of this investigation.
Having seen how the mass and angular momentum moments behave we turn to the scalar field and its moments. The first quantity of interest is the scalar charge or scalar monopole W 0 . The degree of scalarization of a neutron star will depend on the value that we choose for the parameter β, with more negative values leading to more scalarized stars and therefore larger values of the scalar monopole as well as the higher order scalar moments. In Fig. 2  We should recall at this point some general properties of the models that will help the reader interpret the plots. As we mention above, the scalar moments are given in terms of j andM 2 . Increasing value of j corresponds to increasing rotation rate of the star and the higher the degree of scalarization the higher the maximum spin that the models can have. Neutron stars in GR can have a spin up to j max = 0.7 independent of the EOS (see [50,51] for more details) but scalarized stars can have larger spins. Larger values ofM 2 correspond to less compact objects of lower mass (values larger than 10 usually correspond to masses around or less than 1M ), while the more compact objects with masses close to the maximum mass have the smallest value ofM 2 , which tends to 1. The plots show that for largeM 2 the models are not scalarized, while, asM 2 decreases, at some point stars are spontaneously scalarized and the scalar monopole becomes non-zero. Eventually at small enough values ofM 2 , the models become unscalarized again. The degree of scalarization of a neutron star is not independent of the spin. More rapidly rotating neutron stars tend to be more scalarized. To counter this effect to some degree we have chosen to normalise the scalar monopole asW 0 ≡ −W 0 /(j a M ), where the spin weight is a = 0.3. The spin normalisation was introduced initially in the hope of eliminating the spin dependence, but this has not been possible for any value of a. In spite of this, we have decided to keep this normalisation for all the moments in order to minimise their variation due to the spin. This point will be further discussed when it arises again. Fig. 2 shows that, within the same theory, i.e., for the same value of β, all models fall on the same surface independent of the EOS, which means that the scalarized monopole has a universal behaviour. For different theories (different βs), or for different asymptotic values of the scalar field, the surfaces are different. Unfortunately, the surfaces shown in Fig. 2 are not easy to fit with some simple function. Spontaneous scalarization is a phase transition that occurs at a threshold and finding . The scalar charge demonstrates a universal behaviour, i.e., all the EoSs form the same surface. The relevant surface though, changes depending on the value of β of the theory and the asymptotic value of the scalar field φ0. The quantity that is plotted is the reduced scalar chargē W0 ≡ −W0/(j a M ), where a = 0.3. The reason that this scaling that includes the spin parameter was chosen is because the degree of scalarization has also a spin dependance, so the idea was to try to flatten out the effect. The same scaling with respect to j seems to work for all the different theory cases. The two upper plots correspond to exactly the same data (the middle plot is the surface formed by the points of the top plot). some empirical relation that would express this threshold in terms of the moments is not straightforward. One last thing to note is that in the case where the asymptotic value of the scalar field is not zero, the models are scalarized even for small compactnesses (or largeM 2 ) as we can see in the bottom plot of Fig. 2.
We now turn our attention to the next scalar moment, the scalar quadrupole. The behaviour of the reduced scalar quadrupoleW 2 ≡ W 2 /(j a M 3 ) is similar to what we saw for the scalar monopole and is presented in Fig. 3. As for the monopole, we have assigned a spin weight to the normalisation of W 2 which is a = 5/3. One could assume that the scalar quadrupole would be driven by the mass quadrupole of the star and therefore the spin dependance would be ∼ j 2 , but as it turns out, the behaviour is more complicated than that. For this reason we have chosen to normalise the multipole in this way in order to reduce the variation due to the spin, as we did for the scalar monopole,. Similarly to the monopole, different choices of β and ϕ 0 correspond to different surfaces in the parameter space, while all EOSs for the same theory fall on the same surface.
The last scalar moment that we have calculated from the numerical models is the scalar hexadecapole a = 3.6, are given in Fig. 4. Again we observe a behaviour similar to the previous two cases. The bottom line of this analysis is that the scalar moments in the Einstein frame demonstrate an EOS independent behaviour for the same parameters β and ϕ 0 following the same surfaces in the respective parameter spaces, while as we change β the moments fall on clearly separate surfaces. Therefore, while EOS uncertainties can be circumvented just as in GR, as the scalar moments demonstrate universal behaviour, one can identify different theories (different βs) as they correspond to different surfaces forW 0 , W 2 , andW 4 .

V. RELATING MOMENTS TO OBSERVABLES AND COMPARISON TO GR.
In the previous section we saw how the moments in the Einstein frame exhibit universal behaviour with respect to the different EOSs of nuclear matter. We also saw that in the case of mass and rotation moments the universal relations found in GR also capture the behaviour of the moments of scalarized stars independently of the specific theory chosen. The latter property doesn't hold for the scalar moments. While they are EOS independent within a specific theory, they do depend on the choice of a particular theory (reflected on the choice for the value of β). But as we have already mentioned, the Einstein frame moments are not directly observable and if we want to connect our results to astrophysical observations we will have to calculate physical quantities in the Jordan frame. The transformation to the Jordan frame is a conformal transformation of the form,g µν = A 2 (ϕ)g µν , and we have defined earlier in Section II the coupling function k(ϕ) = d ln(A(ϕ))/dϕ.
It is common in the literature to use the Damour & Esposito-Farése notation for the asymptotic expansion of this quantity, i.e., k(ϕ)| ∞ = α, (dk/dϕ) ∞ = β, d 2 k/dϕ 2 ∞ = γ, and so on. In our case and for what follows, due to the form of the coupling function that we have been using, α, γ (in this notation) and all the higher derivatives will be zero. Furthermore, since current constraints point towards a zero asymptotic value of the scalar field we will assume that at infinity we have ϕ 0 = 0. We stress that we have instead used α, γ to denote metric functions above.
Returning to the question of observables, the natural choice is to consider observables that are related to the geodesics of the spacetime. Such observables and their connection to the moments both in GR and in STTs are discussed in what follows.

A. Observables and moments
It was shown by Ryan [52] that there are quantities associated to the geodesics of a GR spacetime that can be expressed in terms of expansions where the coefficients depend on the multipole moments. The expansions of these same quantities have also been calculated for STTs with a massless scalar field [33]. These quantities are : (i) the change of energy per unit mass of a test particle (Ẽ) per logarithmic orbital frequency interval for equatorial circular geodesics, denoted by ∆Ẽ; (ii) the ratio of the periastron precession frequency (Ω r ) of a slightly eccentric equatorial orbit over the orbital frequency (Ω) of the corresponding circular orbit, Ω r /Ω; (iii) the ratio of the nodal precession frequency (Ω z ) of a slightly offequatorial orbit over the orbital frequency of the corresponding circular equatorial orbit, Ω z /Ω. The expansion parameter is U ≡ (M Ω) 1/3 , which corresponds to the orbital velocity of the test particle. The quantityM ≡ M − αW 0 corresponds to the Keplerian mass that one would measure from the motion of a companion star, if the system was part of a binary.
Here, we briefly present these expansions in GR and in scalar-tensor theory, as they were derived in [33,52], up to the same corresponding order in U and taking into account the constraints we have from the ansatz that we have used. The energy change per logarithmic orbital frequency change in GR up to O U 8 is given by the expression, while the corresponding expression in scalar-tensor theory is, after setting α = γ = 0 as discussed above, where we haveM = M − αW 0 = M . Similarly, the ratio Ω r /Ω in GR is, while the corresponding expression in scalar-tensor theory is Finally, the ratio Ω z /Ω in GR is, while the corresponding expression in scalar-tensor theory is ∆Ẽ is a quantity that is more immediately relevant to gravitational waves and extreme mass ratio inspirals, while the other two quantities can be also relevant to systems such as X-ray binaries, where one observes quasi-periodic oscillations (QPOs) of the X-ray spectrum of the accretion disc around the compact object. If one were to assume, for example, the relativistic precession model for QPOs, by Stella and Vietri [53,54], then one could associate specific QPO frequencies to Ω r , Ω z , and Ω. 4 The relevant observations could then be fitted to recover the coefficients of the expansions.
B. Measuring the scalar charge and β.
Inspecting the expansions in GR and the corresponding expansions in STT reveals that it is possible to distinguish between the two theories, either by comparing the coefficients of the same order between the two theories or by comparing different order coefficients against each other. As we saw in the previous subsection one could expand Ω r /Ω and Ω z /Ω in terms of powers of Ω as, (Ω r /Ω) = C a Ω a/3 , and (Ω z /Ω) = F a Ω a/3 , where the coefficients C a , for example, will be C a =M a/3 f a (M , β, W 0 , S 1 , M 2 , W 2 , . . .). These coefficients could be used to measure the various parameters. The frequencies that are most commonly observed in low mass X-ray binaries (LMXBs) are the two larger ones, i.e., Ω and Ω r . These are observed as pairs of kHz QPOs, while occasionally one also observes a third low frequency QPO, which is assumed to be Ω z . Since the most common occurrence is the former one, we will start assuming that only Ω r and Ω are known. We will then explore how far one can go by using either additional information from Ω z or the universal behaviour we have described previously.

Setting up the problem and constraints.
In GR one could independently measure the mass from the lowest order term in Ω r /Ω, since we have that C GR 2 = 3M 2/3 . In scalar-tensor theories however that term has additional contributions due to the scalar field and is of the form C STT 2 = 3 − βW 2 0 2M 2 M 2/3 . One could go around this problem if an independent measurement of the massM were available. For example, since this sort of QPO producing X-ray sources are LMXBs, the mass could be estimater from the Keplerian motion of the companion and the compact object (M is the Keplerian mass after all). In that case, the estimation of C STT 2 would provide a measurement of the combination βW 2 0 , but more importantly would immediately tell us that we have a deviation from GR. In GR the higher order coefficients would enable us to measure the higher order moments, while along the way we would find coefficients that would serve as consistency checks, such as C GR 5 which is a consistency check on the measurement of S 1 from C GR 3 . In scalar-tensor theories things are a little more complicated. The coefficients C STT  4 There are other models as well, such as some of the models derived from discoseismology, where oscillations of the disc can be associated to the geodesic frequencies [55][56][57][58][59]. knew the massM independently, but they could also be used to determine the mass since if we combine them we can arrive to the expression, which relates the mass to these coefficients. Therefore, even for a system where the mass is unknown, one can estimate it as long as one can accurately estimate the coefficients up to C STT

5
. This then allows to estimate S 1 as well from C STT 3 . Up to this point we haveM , βW 2 0 , and S 1 . Turning to the coefficient C STT 4 that contains M 2 we notice that we cannot estimate it independently. We can only estimate it in combination with W 0 , i.e., 9 2 The problem lies with our inability so far to separate β and W 0 . Aiming to break the degeneracy between M 2 and W 0 by using higher order terms seems a difficult task with uncertain conclusion. For instance, while the next order term, C STT 6 , includes all the relevant terms, it also includes the scalar quadrupole W 2 that first appears in the expansion at that order.
The situation for measuring the multipole moments and the parameters of the particular STT improves dramatically when we have information for both Ω r /Ω and Ω z /Ω from a specific system. In that case, we can use the same analysis presented for Ω r /Ω to estimateM , βW 2 0 , S 1 , and the combina- , not present in GR, can serve as an independent verification of the deviation from GR, as well as a consistency check up to that order. Therefore additional information from more frequencies allows for the breaking of degeneracies and performing more tests on deviations from GR.
2. Using universal relations to overcome degeneracies.
In the above discussion we showed that in order to break the M 2 − W 0 degeneracy one would have to consider both Ω r /Ω and Ω z /Ω, but to reach to that conclusion we did not take into account the results of Section IV and the universal behaviour of W 0 /M . In fact, if we were to consider only Ω r /Ω, the universal behaviour of W 0 /M and the fact that it can be expressed as some function of j and M 2 , could be used to break the M 2 − W 0 degeneracy even without considering Ω z /Ω. In what follows we will describe the algorithm that can be used to do this.
Lets assume that from observations of Ω r and Ω we have estimated the first coefficients of the Ω r /Ω expansion, i.e., C STT will provide the massM of the neutron star, as we describe above. The additional information that we have is that, and W0 using the knowledge of the massM (from eq. 38), the angular momentum S1 (from eq. 40), the constraint for βW 2 0 from eq. (39), the constraint from eq. (41), and the universal relation betweenW0, the spin parameter j = S1/M 2 , and the reduced quadrupoleM2, as it is given in Figure 2.
The plot on the left shows eq. (39) for different values of β (horizontal lines), the universal relationW0 = f (j,M2) for the value of the spin j that has been estimated (in this case j = 0.35) and for different values of β, and eq. (41) which corresponds to the red dashed line. If the results are consistent, then the universal relation cross section curve and the eq. (39) curve that correspond to the same β, should intersect with the eq. (41) curve at the same point. This then indicates the value of β and the corresponding values ofM2 andW0. The plot on the right is a magnification of the region where the three curves intersect. The kinks on the curves are due to interpolation errors. where,C Eq.(40) straightforwardly gives the spin parameter of the compact object, while Eq. (39) can be interpreted as a bond between W 0 and β and Eq. (41) relates W 0 to M 2 . Hence, one needs one more bond between these quantities in order to be able to determine them uniquely. The universal relation provide it as follows. First one expresses the equation in terms of the variables of Section IV by using the relations W 0 /M 2 = j 0.6W 2 0 and (M 2 /M 3 ) = −j 2M 2 . One can then effectively consider all quantities as having being uniquely determined, except ofW 0 andM 2 , that instead just satisfy a bond. The additional bond is provided by the universal relation betweenW 0 and (j,M 2 ) shown in Figure 2. For the estimated value of j, the result is a cross section of the surfaces shown in Figure 2, which amounts to having for different values of β different curves relatingW 0 toM 2 .
We can now plot all these constraints on aW 0 −M 2 plot, an example of which is shown in Figure 5. The plot shows theW 0 =const. lines that result from the constraint (39) for the different values of β = −4, −4.5, −5, and −6. It also shows the universalW 0 = f (M 2 ) curves for the corresponding values of β (we note that for the spin in this example, the β = −4 models are unscalarized). The last curve shown is the constraint resulting from Eq. (41) which is independent of β and is therefore a single curve (red dashed curve). In order to have a consistent solution of all of the constraints, the curves of Eq. (39) and Eq.W 0 = f (M 2 ) that correspond to the same β must intersect the curve for Eq. (41) at the same point (or approximately the same point), just as the example in Figure 5 shows. From the intersection point one can identify the value of β (= −5 in this example) as well as the values ofM 2 (∼ 2) andW 0 (= 0.4).
We have therefore presented an algorithm that makes use of the universal relations for the scalar monopole, in order to measure the parameters of a given scalarized neutron star and the parameters of the corresponding STT from a set of astrophysical observations (pairs of QPO frequencies in this case) that otherwise would not be possible due to degeneracies.

VI. CONCLUSIONS
Universal or EOS independent relations between global properties of neutron stars have proven to be a versatile tool for inferring the properties of neutron stars. These relations have been extensively studied in GR and particular flavours of them, such as the I-Love-Q relations, have been studied in a variety of modifications to GR. The 3-hair universal relations found in GR have been a more difficult problem to tackle, mainly due to the intricacies of defining multipole moments in modifications to GR. STT of gravity with a massless scalar field is a class of theories where a definition of moments is already available [32].
In this work we have computed the multipole moments for scalarized stars in these theories and have shown that they continue to exhibit universal properties. Specifically, we have found that the mass and angular momentum moments follow the same universal 3-hair relations as their GR counterparts [20][21][22], independent of the value of the β parameter and of the asymptotic value ϕ 0 of the scalar field. Furthermore we have found that the scalar field moments for every given combination of β and ϕ 0 exhibit universal behaviour in terms of the spin parameter j and the reduced quadrupoleM 2 . That is, when each moment is plotted in terms of j,M 2 it falls on the same surface independent of the EOS. In addition, different values of β and ϕ 0 result in different surfaces in the three dimensional space formed by each scalar moment and the two parameters j,M 2 . This appears to be related with the known fact that the degree of scalarizations depends on both the asymptotic value of the scalar field and the value of β.
Our results demonstrate that the degree of scalarization can be expressed in an EOS independent way, which is still quite intriguing and potentially very useful. In particular, we demonstrate how one can use the universal relations presented here to infer the various properties (i.e., the moments) of a scalarized neutron star, as well as the parameters of the specific STT (i.e., the value of β), from astrophysical observations.
The algorithm for doing so using LMXBs can be seen as a proof of principle. A more thorough analysis is necessary in order to determine how accurately the various parameters can be measured and what sort of constrains can be set on STTs from observations. Furthermore, it would be worth exploring how the results presented here could be used in other settings, such as the observation of gravitational waves from the inspiral of NS-NS binaries or BH-NS binaries.
STT with a massless scalar is only a first step in studying the 3-hair relations in these theories. The next would be to consider STTs with a massive scalar field. This is more challenging, as the multipole moments cannot be defined in the same way as in the massless case. Nevertheless some cases can be easier to handle than others. For instance, if the scalar field were to be fully confined inside the neutron star on account of its large mass, then in the exterior for all practical purposes one would only have to deal with the spacetime, and the moments would be calculated in the same way as in GR. This is something that we will explore in future work.

ACKNOWLEDGMENTS
The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC Grant Agreement n. 306425 "Challenging General Relativity". GP acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant agreement no. DarkGRA-757480. DD would like to thank the European Social Fund, the Ministry of Science, Research and the Arts Baden-Württemberg for the support. DD is indebted to the Baden-Württemberg Stiftung for the financial support of this research project by the Eliteprogramme for Postdocs. TPS acknowledges partial support from the STFC Consolidated Grant No. ST/P000703/1. SY acknowledges financial support by the Sofia University Grants No 3258/2017 and the Bulgarian NSF Grant DCOST 01/6. Networking support by the COST Actions CA15117, CA16104 and CA16214 is also gratefully acknowledged.

Appendix A: Calculation of the multipole moments
The calculation of the equilibrium neutron star solutions is done using a modification of the RNS code (see [34] for the original GR version of the RNS code while the STT extension can be found in [16]) and that is why we will follow the formalism and notations that are standard for the KEH method [35,36]. The coefficients B 2l , ν 2n,0 , ω 2n−1,0 and Φ 2n are calculated numerically using integrals of the source functions of the field equations (3)-(5), (8). In the present paper we consider the case of a zero mass scalar field and thus these integrals are the same in pure GR and in STTs (when the Einstein frame is employed). 5 First we describe the calculation of the coefficients B 2l . The function B in the metric ansatz (12) where s is the compacted radial coordinate (1 − s) /s = r eq /r with r eq being a characteristic length scale that gives the coordinate equatorial radius of the star, and µ = cos θ.
The source termS γ (s , µ ) is connected to the right hand side of eq. (3) and is given bỹ For simplicity, the expression (A2) can be written as γ(r, µ) = − 2e −γ/2 π Using the orthogonality conditions of the Gegenbauer polynomials to relate the Γ 2n coefficients to the B 2l coefficients, we have that the B 2n coefficients are given as The only term that will survive the integration is the k −1 = n term that will give the result, with the Γ 2l coefficients given by (A5). The calculation of the rest of the expansion coefficients ν 2n,0 , ω 2n−1,0 and Φ 2n is more straightforward. Thus, in terms of the source functions in the field equations (3)-(5) used by the RNS code, these coefficients are given by