Quasi-universal relations for generalized Skyrme stars

First proposed in 2013 by Yagi and Yunes, the quasi-universal \emph{I-Love-Q relations} consist of a set of relations between the moment of inertia, the spin-induced quadrupole moment and the electric quadrupolar tidal deformability of neutron stars which are independent of the Equation of State (EoS) within an accuracy of $\sim1\%$. In this work, we show that these relations hold for different Skyrme-based nuclear matter EoS and also for the star-like solutions of different Einstein-BPS-Skyrme-models, some of which do not even present a barotropic equation of state. Further, other quasi-universal relations are analyzed, and together with recent GW observations, we use them to select the generalized Skyrme model that better reproduces observations. Our results reaffirm both the universality of the \emph{I-Love-Q} relations and the suitability of generalized Skyrme models to describe nuclear matter inside neutron stars.

The Skyrme model [1] and its generalizations [2][3][4][5][6][7] consist in a set of relativistic, effective-field-theoretic models of interacting Goldstone bosons which have been proposed to describe strongly interacting matter in a low energy regime. Indeed, baryons, nucleons and nuclei, whose existence can not be inferred by perturbative QCD methods, are described within the Skyrme models as (topological) solitonic configurations of the underlying bosonic degrees of freedom. During the last few decades, many properties of nucleons [8] and nuclei [9][10][11][12][13][14][15][16] have been reproduced using these models. These results have contributed to establish the Skyrme model approach as a well-motivated proposal for the description of nuclear matter. Furthermore, in recent years there has been a growing interest in obtaining self-gravitating solutions of the Einstein-Skyrme system in order to determine whether the Skyrme model and its generalizations are also able to describe the properties of matter inside neutron stars (NS) [17][18][19][20][21][22].
On the other hand, one of the most outstanding challenges of current astrophysical research is to obtain information about the equation of state of ultra-dense matter from neutron star observations. In particular, apart from their masses and radii, other interesting observable properties of NS are their quadrupole moments, spin angular velocity (angular momentum), and deformability against tidal forces-which is encoded in the so-called Love numbers [23,24]. All these properties can be constrained by their imprints into the waveform of a gravitational wave arXiv:2011.08573v1 [hep-th] 17 Nov 2020 signal emitted by an inspiraling binary neutron star system. Indeed, binary NS systems are one of the most promising sources of gravitational waves (GW) within the detection range for second generation observatories, such as Advanced LIGO, Advanced VIRGO [25,26], or KAGRA [27]. The observation of GW emitted during the coalescence of the stars in such systems-especially in the last part of the merging, in which the stars are subject to large tidal deformations due to the extremely strong gravitational fields involved-will shed light onto the equation of state (EoS) of matter at very high densities, well beyond the nuclear saturation point.
A particularly interesting property of compact stars is the apparently universal relation between the moment of inertia, the Love numbers and the quadrupole moment (I-Love-Q relations) of such stars. These I-Love-Q relations, firstly proposed by K. Yagi and N. Yunes in [28], when applied to NS, allow to break the degeneracy between the quadrupole moment and the NS spins in the gravitational waveforms of inspiraling NS binaries. Therefore, a much more precise determination of the (dimensionless) averaged spin can be reached in such measurements [29].
In this paper we show that star-like solitonic solutions of different Skyrme-type models not only exist, but also reproduce to a good extent some of the currently best known properties of NS-like the typical values of mass, radius, moment of inertia, Love numbers, etccoming from astrophysical measurements, GW observations and/or computer simulations.
We also address the issue of whether the compact star solutions obtained within different Skyrme-based models and the corresponding EoS satisfy the I-Love-Q relations, and find that they indeed do satisfy them, even though the equations of state for different models present big differences. The Skyrme model, being relatively simpler than other phenomenological or first-principle based relativistic field theories describing nuclear matter, therefore not only stands as an excellent candidate to describe nuclear matter at very high densities such as those inside NS. In addition, it provides a simple toolkit for the construction of a wide range of models of nuclear matter and their corresponding EoS which allow to investigate the resulting NS properties and universal relations in environments not considered previously-like, e.g., for nonbarotropic EoS, see below.
The rest of the paper is organized as follows: In the second section, we present the generalized Skyrme model and the submodels that we want to consider, and find static solutions to the Einstein equations for these models coupled to gravity. For all models, we choose the model parameters such that the resulting skyrmionic matter approaches the nuclear saturation density n 0 = 0.16 fm −3 and the energy per baryon E 0 = 923.3 MeV of infinite nuclear matter at saturation in the limit of zero pressure, as in [22]. In sections III and IV, we review the framework for obtaining perturbative solutions to the Einstein equations that represent, respectively, slowly rotating and tidally deformed stars. We find the equations for the metric perturbations up to second order in spin (for the rotating case) and to linear order in the external tidal field in the tidally deformed case, and solve these numerically for the star interior, after which the matching with the analytical exterior solutions is performed in order to obtain the correct values of the first multipoles in the asymptotic expansion of the metric outside the star. Also, in section IV the definition of the Love numbers is given and a procedure to obtain the first (quadrupolar) Love numbers from the perturbative analysis is explained. Finally, in section V, we show different quasiuniversal relations between the dimensionless moment of inertia, quadrupolar moment, electric and magnetic Love numbers and compactness of the stars for all the models at hand, and discuss them. We conclude with a comment in section VI about some constraints on the deformability of NS coming from recent GW observations of binary mergers and on the feasibility of the different Skyrmebased models to describe nuclear matter inside NS, and finish in the last section with a summary of all these results and prospects of future work. In our choice of units, the speed of light is c = 1. For masses (energies) and lengths we use either solar masses M and km-for astrophysical objects, or MeV and fm-for nuclear physics observables.

II. STATIC STARS AND GENERALIZED SKYRME MODEL
A. The Generalized Skyrme model The original Skyrme model is defined by the following Lagrangian Tr{L µ L µ } (1) where the Goldstone bosons associated to chiral symmetry breaking-the lightest degrees of freedom of QCDform the SU (2) matrix Skyrme field U (x). Furthermore, L µ = U † ∂ µ U is the su(2)−valued, left invariant Maurer-Cartan form and U is a non-derivative part of the model, i.e., a potential. This theory possesses only three coupling constants f π , e and µ which are, respectively, the pion decay constant, the Skyrme coupling constant and the µ parameter related with the pion mass via µ = m π f π / √ 8. More precisely, m π is the mass of small perturbations around the vacuum (pions) if the potential tends to the pion potential U π = 1/2Tr(1 − U ) for U → I.
For finite energy solutions, it is necessary to impose constant boundary values of U at |x| → ∞. This implies the appearance of a nontrivial topology. Indeed, the physically relevant matter field configurations define maps which are classified by an integer number or topological degree where B µ is the topological current. Importantly, it can be rigorously proven that the topological charge is just the baryon charge [30]. Due to this equivalence, topological solitons with a non-zero value of the topological charge, typically referred to as Skyrmions, are identified as baryons and atomic nuclei, although the derivation of their properties as quantum systems from the Skyrme model (a non-renormalizable field theory) requires a careful quantization of the zero as well as massive (vibrational) modes of the classical solitons [8,13,16]. It is one of the most attractive features of the model that all these objects are emergent phenomena which arise from a very simple Lagrangian based entirely on pionic degrees of freedom, which contains a very small number of parameters.
Being an effective theory, the Skyrme model can be extended by adding higher order terms to the Lagrangian. It can be shown that the only possible Lorentz-invariant extra term with at most second order time derivatives of the Skyrme field is [6] where λ is an additional coupling parameter related to the ω-vector meson. Indeed, this term can be obtained by integrating out this vector meson from a model that includes both pions and vector mesons [3]. As a result, we get the generalized Skyrme model For some choices of parameters and potentials [31,32], this model maintains the successes of the standard Skyrme model (1) in the description of vibrationalrotational spectra of some light nuclei, but also leads to physical binding energies, which result too large in the standard Skyrme model. An obvious, next step of application of the Skyrme model is to study properties of nuclear matter at extremely high densities, e.g., to describe the equation of state of neutron stars (see [21] for a recent review). To do so, we have to find the lowest energy solutions of the Skyrme model for a topological charge of the same order as the total baryon number of neutron stars, which typically is N ∼ 10 57 . Basically, within the framework of the generalized Skyrme model, there are two qualitatively distinct possibilities.
First of all, it is well known [33,34] that the lowest energy solutions for the standard Skyrme model (1) with an arbitrarily large baryon number, B → ∞, consist of crystalline structures of Skyrmions. The EoS of this Skyrme crystal may then be used as a starting point for the investigation of NS. A second possibility is related to the fact that there is a very special point in the space of the model parameters, resulting in the so-called BPS Skyrme submodel: L BP S = L 6 + L 0 . The name comes from the fact that this Skyrme model supports topological soliton configurations which saturate the BPS energy bound [6], offering a possibility to resolve the problem of the unphysically high binding energies of the standard Skyrme model. What is more important here, this model describes a perfect fluid for any value of the baryon charge. Indeed, the stress-energy tensor reads [19,20] where the four-velocity u µ , pressure p and energy density ρ are (7) As the sextic term L 6 gives the leading contribution to the energy at high pressure/density [35], the fluid behavior is expected to dominate at this regime. This is, of course, consistent with the usual understanding of the inner core of neutron stars as being formed by a fluid of neutron matter. These two states of the Skyrmionic matter should be smoothly joined in the generalized model, suggesting a phase transition as pressure increases. Although the existence and properties of such a phase transition in the full model is still an unsolved problem, the known regimes at low and high pressure have recently led to the proposal of a generalized EoS [22], see below.

B. Static NS solutions
Here and in the following sections, we will obtain solutions to the Einstein equations that describe NS within the different Skyrme models presented above. As a first step, we will consider static, spherically symmetric configurations, which is usually done following the Tolman-Oppenheimer-Volkoff (TOV) approach, in which the Einstein equations are solved using the stress-energy tensor of a perfect fluid. Thus, we suppose the spherically symmetric (Schwarzschild) ansatz for the metric, ds 2 = −e α(r) dt 2 + e β(r) dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ). (8) We extract from the Einstein equations and the conservation of the stress-energy tensor of the perfect fluid type (∇ µ T µ ν = 0) the following system of ODEs, also known as the TOV system, where we have made the usual definition exp(−β) = 1 − 2M/r (11) so that the value M * = M (R * ) of the function M = M (r) coincides with the (static) ADM mass of the star when evaluated at its radius r = R * . To close the system (10), we have to know the relations between the pressure and the energy density, i.e., an EoS. It is at this point where the classical Skyrmion solutions with a very large value of the topological charge become relevant.

C. Skyrme neutron stars
Next, we briefly review the current status of the description of static properties of neutron stars from the Skyrme model perspective.

The LSK Skyrme neutron stars
The usual Skyrme model is a field theory whose energymomentum tensor does not have a perfect fluid form. Therefore, a suitable mean-field approximation has to be performed. In practice, it means a spatial averaging. The ground state is a crystal with a given lattice structure and lattice spacing l 0 (we assume the isotropic case). Obviously, the energy per baryon E(l) has a minimum at l = l 0 . This solution is also a zero-pressure (equilibrium) solution, because where V = l 3 is the volume of the cell. Diminishing the lattice spacing l is equivalent to imposing a nonzero pressure. Finally, as the pressure and the energy density are both functions of l, we can find the corresponding EoS, ρ SK = ρ SK (p). If inserted into the TOV system, the crystal EoS amounts to neutron stars with rather small maximal masses, significantly below the observed NS masses. For example, for the cubic, face-centered lattice of B = 4 Skyrmions (α particles) M max 1.49M [17,21]. The corresponding mass-radius curve is presented in Fig. 1, the pink dots.

The LBP S Skyrme neutron stars
In the case of the BPS Skyrme submodel L BP S , which is a genuine perfect fluid theory for any potential U, one can find lowest energy Skyrmions for any value of the topological charge B in an exact form. There are, in fact, infinitely many solutions for a given B related via SDif f diffeomorphisms, which corresponds very well with the fluid nature of the BPS Skyrmions. Interestingly, the perfect fluid form of the action allows to obtain the mean field EoS in an exact form without solving the field equations [20,36]. This occurs because the pressure enters as an integration constant into the generalized Bogomolny equation. As a consequence, the pressure dependence of both the energy E(p) and the volume V (p) of BPS Skyrmions can be found as target space integrals (averages). The details of the resulting EoS obviously depend on the particular choice of the potential (but, of course, do not depend on a particular solution). On the other hand, since the sextic term provides the leading behavior in the high pressure limit, the EoS tends to the maximally stiff equation of state as the pressure increases As a consequence of this stiffness, it is not surprising that the neutron stars provided by the BPS Skyrme model have rather big maximal masses, easily exceeding 3Msee Fig. 1, black, green, purple and blue dots, which correspond to the four different potentials introduced in [20,36], namely the θ-potential U Θ = Θ(Tr{1 − U }), the standard pion-mass potential U π = 1/2Tr{1 − U } = 2χ(r), the pion-mass squared potential U 2 π = 4χ(r) 2 and the partially flat potential Owing to its perfect fluid nature, the BPS model offers the possibility to close the TOV system without any mean-field approximation. In this case, referred to as the exact case, the pressure and energy densities ρ, p can already be read from the stress-energy tensor (6). Furthermore, they are related in a non-algebraic way, by construction (7). This also means that the obtained matter is an example of a non-barotropic fluid where constant pressure do not correspond with constant energy density. Hence, this exact approach may serve as a laboratory where the impact of non-barotropic EoS on properties of NS can be studied. Further, the different BPS models provide a wealth of new and different EoS which will allow us to test the universal, EoS-independent character of certain relations, like the I-Love-Q relations, in new environments not considered previously.
More precisely, in the exact case the Skyrme field U enters in the Einstein equations as an additional degree of freedom, so that we have to obtain its own differential equation in order to close the TOV system. To do this, we choose the hedgehog ansatz for the Skyrme field, U (x) = e iξ(r)n(θ,φ)· σ , n(θ, φ) = (sin θ cos (Bφ), sin θ sin (Bφ), cos θ), which is compatible with the chosen ansatz of the metric, since it yields a spherically symmetric energy density, which is relevant for static NS. Here σ are the Pauli matrices and (r, θ, φ) are spherical coordinates. The only degree of freedom in this ansatz corresponds to the radial profile ξ(r), and inserting the hedgehog ansatz into the definition of p it can be shown that this function satisfies the differential equation which is added to (10) to close the system. For simplicity, when solving the TOV system we will define the new variable χ := sin 2 (ξ/2), which satisfies Once the system of ODEs is closed, only a set of initial conditions are needed as an input in order to obtain a particular solution. However, in the exact case, the baryon number B of the star is an additional input parameter, and the value of the pressure at the center of the star (p 0 ) that yields the input value must be found via a shooting method, with initial conditions requiring that the pressure vanishes at some finite value p(r = R * ) = 0. This value R * is precisely the radius of the star. The value of α 0 is not needed to solve the system. However, only one value is correct, and it can be obtained by imposing continuity of the metric at the radius of the star, R * , for which, and onwards, the metric is given by the Schwarzschild solution: Also, the central value of the energy density ρ 0 is determined by the BPS EoS (7).
On the other hand, in the Skyrme crystal and the mean-field version of the BPS submodels, we do have a barotropic EoS ρ(p), so that the energy density only depends on the pressure. In these cases, the equation (17) is no longer needed and the input parameter is the pressure in the center of the star p 0 , along with the rest of initial conditions for α and M . The system of differential equations is then solved up to the star radius (R * ), in that point the static ADM mass of the star M * = M (R * ) is also obtained.
In Fig. 1 mass-radius curves for the exact case are presented -see green, blue and purple stars. For the θpotential the MF and exact computations obviously coincide. Therefore, for relatively flat potentials (e.g., the pion-mass and the partially flat potential) the difference between the MF and exact approach is rather small, while it strongly increases for more peaked potentials (e.g., the pion-mass potential squared).

Neutron stars and the generalized EoS
As we see, the usual Skyrme model crystal and the BPS Skyrme fluid result in too small or too large maximal masses of neutron stars, respectively. It can be expected that these two extremal cases can be balanced in the full generalized Skyrme model. While the EoS for the generalized Skyrme model is not currently available, it motivates the following generalized Skyrme EoS which interpolates between the crystal and fluid phases [22] ρ Gen (p) = (1 − α(p))ρ SK + α(p)(p + ρ SK (p P T )) (20) where the interpolating function tends from 0 for p/p P T → 0 to 1 for p/p P T → ∞. The parameter p P T can be identified with the position of the crystal/fluid phase transition, whereas β measures how rapid the transition occurs. Specifically, we assumed a rather gradual phase transition (β = 0.9) located at p P T ∈ (25, 50)MeV/fm 3 in [22]. We remark that the value of p P T strongly affects the maximal mass.
In Fig. 1 we show the mass-radius curve for the generalized EoS with p P T = 25MeV/fm3 3 -see yellow squares. As expected, the maximal mass of NS is between the two previously discussed versions of the Skyrme model and reads M max 2.55M .

Neutron star crusts and the hybrid EoS
By construction, the generalized Skyrme model contains only pionic degrees of freedom (with some other heavier mesons effectively also taken into account). This means that it is relevant for describing nuclear matter above the saturation density. For lower densities, the electromagnetic interaction starts to have a nontrivial impact on the properties of nuclear matter, leading to the appearance of inhomogeneous phases (such as "nuclear pasta" phases [37]). Although the Skyrme model can be coupled with the electromagnetic U (1) gauge field, which in principle may allow to study such phases within the framework of the Skyrme model, the resulting theory is very complicated and currently no large B Skyrmions are known. However, it is possible to take into account this low density regime, relevant for the crust region of NS, by assuming a transition of the generalized Skyrme EoS to a standard nuclear EoS obtained by the usual manybody techniques. Concretely, we choose the EoS ρ BCPM of [38], as we did in [22]. As a consequence, we arrive at a hybrid EoS (22) where now β = 2 and the position of the transition p * ∈ [0.5, 2]MeV/fm 3 . Further, α(p, p * , β) is defined in (21). The resulting mass-radius curve is presented in Fig. 1, olive squares (for p * = 1 MeV fm −3 and p P T = 25 MeV fm −3 ). We want to emphasize that the NS resulting from the hybrid EoS (22) pass all current observational constraints.

III. SLOWLY ROTATING SKYRME STARS
In this section, we will study how the previously obtained spherically symmetric Skyrmion stars-models of NS based on the (generalized) Skyrme models considered in the present paper-behave under small perturbations. In particular, we will analyze their deformation due to (small) rotation and tidal forces, in order to obtain useful relations between their moment of inertia, deformability and Love number-known as I-Love-Q relations, first proposed in [28]-which may help in extracting information about the internal structure of compact stars. Throughout this and subsequent sections, we will largely follow the approach and notation of [28].

A. Slowly rotating stars: Hartle-Thorne formalism
To analyze the properties of rotating Skyrmion stars, we will make use of the Hartle-Thorne formalism for slowly rotating stars, introduced in [39]. This formalism establishes a perturbative framework which consists in an expansion of the metric in powers of a perturbation parameter-related with the rotational frequencyand solving the Einstein equations order by order in this parameter. This perturbative expansion has proven particularly useful in the literature since it allows to obtain approximate solutions to the Einstein equations both for the interior and exterior of the star, hence, it enables to retrieve information about the equation of state for the matter inside the star from the multipolar expansion of the external solution. We will now review the procedure to obtain the solution in this approximation for the metric in the interior of a compact star, and in the following sections we will do the same for the exterior solution and the matching between both solutions at the star surface.
The starting point of the slow rotation approximation is to consider a static solution for the metric of a nonrotating configuration, and subsequently add perturbation terms up to a given order in a suitable parameter related to the spin of the star. In our case, we will start from the static metric with line element (8) and, as in [40], defining the spin parameter = Ω * /Ω K in terms of Ω * -the angular velocity of the star as measured by an external, static observer located at spatial infinityand the characteristic frequency where M 0 and R 0 are the non-spinning mass and radius of the star. The characteristic frequency Ω k corresponds to the Keplerian orbital period of a test particle at a radius R 0 around a mass M 0 and thus can be thought of as the rotational frequency for which the mass shedding occurs, i.e., an upper limit for the rotational frequency of the star [41]. For spin frequencies much smaller than this characteristic frequency, the parameter serves as a suitable small parameter about which we can expand the metric. On the other hand, for spin frequencies near the Keplerian limit, ∼ 1 and the Hartle-Thorne approximation is no longer valid. Despite the dependence of the Keplerian frequency on the EoS, the slow-rotation approximation is valid for even the most rapidly spinning neutron stars observed to date [40].
Therefore, let us consider the background spacetime whose metric is given by the static line element (8). We now extend this metric by defining a one-parameter family of metrics g( ) whose components may be expanded in powers of , g( ) = g (0) + g (1) + 1 2 2 g (2) + · · · , with g (0) given by (8). Note that this construction introduces an inherent gauge freedom (for details see, for example, [42,43]) Thereby, following [28], up to second order in , we may write the metric of a slowly rotating star in the Regge-Wheeler gauge as: whereω =ω(θ, r),h =h(θ, r),m =m(θ, r),k =k(θ, r), andM (r) is related toβ(r) in the same form as in (11).
Comparing with the general expansion of g( ), we find: g (2) = − 4eᾱh + 2r 2 sin 2 θω 2 dt 2 + + 4eβm r − 2M dr 2 + 4r 2k (dθ 2 + sin 2 θdφ 2 ). (25) Note that the metric perturbation functionω enters at first order in the spin parameter, whereash,m and k correspond to second order perturbations. This can be easily understood with the following argument [39]: a transformation of the metric for a stationary and axially symmetric rotating spacetime of the form Ω → −Ω should be equivalent to t → −t. This, in particular, implies that an expansion of the diagonal components of the metric in powers of must contain only even powers (since they are unchanged under time reversal), whilst an expansion of the g 0,3 term will only contain odd powers of . Furthermore, sinceω corresponds essentially with the g 03 term of the metric, it is responsible for the dragging of inertial frames. In other words, it measures the rate of rotation that a freely falling observer would undergo with respect to a static one (Lense-Thirring effect). Due to these perturbation terms in the spacetime metric, both the Einstein tensor for the metric and the stressenergy tensor for the matter field will develop perturbation terms, as well. Indeed, just as with the metric tensor, we may define the one-parameter families of perturbed quantities G µν ( ) and T µν ( ), expand them in powers of and impose that Einstein equations are satisfied order by order in the expansion parameter. In particular, both the pressure and mass densities of the matter field will be perturbed, acquiring an angular dependence, i.e.
as well as the fluid four-velocity, u( ). For this latter quantity, we further impose the normalization condition g( ) µν u µ ( )u ν ( ) = 1. Also, stationarity, axial symmetry and rigidity of the fluid flow requires u( ) to be proportional to both killing vectors, i.e. u( ) = f 1 ( )(∂ t + f 2 ( )∂ φ ). The f 1 function is obtained by the normalization condition at each order, and, since the background configuration corresponds to a static fluid, thus the constant C corresponds to the angular velocity of the fluid as measured within the inner coordinate system. Note also that only odd powers of enter the expansion of f 2 , for the same symmetry arguments as for ω.
It is important to notice that all these (one-parameter families of) objects so defined are gauge-dependent, although the Einstein equations themselves do not depend on the gauge (i.e, they must be fulfilled in any gauge). We thus may take advantage of this gauge freedom to choose the most convenient form of the metric functions. In particular, we may choose C = Ω K in (28), so that the coordinate system in the interior of the star is taken to be that of a static observer which measures the angular velocity of the fluid to be du φ /du t = εΩ K = Ω * . It can be shown that any other choice of the constant C = C 0 is equivalent to a gauge transformation of the first order metric perturbation defined by the vector V = (Ω K −C 0 )t∂ φ [43].
On the other hand, the coordinate system we have chosen so far is not quite well suited to perform the integration of the Einstein field equations from the inside of the star, for the following reason: in order to find a numerical solution for the interior metric, we will have to solve Einstein equations with a non-vanishing stress-energy tensor up to the surface of the star, which is usually defined by the surface of vanishing pressure. While in the spherically symmetric case the surfaces of constant density (or pressure) are trivially those of constant radial coordinate, this is no longer the case once the second order perturbations of the metric due to rotation are taken into account. Indeed, the (perturbed) pressure and mass densities (27) will depend both on r and θ, so that the surface of the star will be deformed with respect to the static case.
Therefore, we will consider a choice of gauge in which the surfaces of constant pressure (density) of the perturbed configuration are those of constant radial coordinate. This is in fact equivalent to a change of coordinates in the perturbed configuration from the original (background) coordinate system {t, r, θ, φ} to another, {t,r, θ, φ}, in which the new radial coordinate is defined by so that r coincides withr in the background configuration ( = 0), whilst the function ζ(r, θ) measures the deviation from spherical symmetry of the perturbed configurations. The new radial coordinater is defined so that p 0 (r) = const defines the isobaric surfaces of the rotating star.
Strictly speaking, one could think that, in the exact case, also the perturbations of the Skyrme profile function χ must be taken into account. However, these will be by construction directly related to the energy and pressure perturbations, and, since we will get rid of these perturbations by a suitable radial coordinate change, also the perturbation on the radial Skyrme profile will disappear. We have checked that this is in fact the case, and that no extra degrees of freedom appear in the perturbative formalism for the exact BPS Skyrme case up to second order in .
In the new coordinate system, the metric (23) is rewritten, up to second order in : where all the metric functions are written as functions of r (and possibly θ), and the denotes a derivative with respect tor. The metric (30) has a rather complicated form. However, we may simplify it by redefining, the metric func-tions: coincides with (30) up to second order in . Although both metrics (32) and (23) are different, they are related through a gauge transformation, so that both must satisfy Einstein equations, and the gauge-independent results obtained in both approaches must be the same (at least, up to second order in ). Note that these metrics are compatible with the general form for the Hartle-Thorne metric in an arbitrary gauge, obtained in [43], which have two commuting killing vector fields Although a priori the metric perturbation functions can have an arbitrary dependence on r and θ, an expansion of these functions is always possible in spherical harmonics [39]. Moreover, the angular dependence of the perturbation functions may be further reduced by additional arguments. For example, axial and reflection symmetry in the equatorial plane implies that the m (axial) number in the spherical harmonic expansion does not play any role, so that it may be reduced to an expansion in terms of Legendre polynomials. Therefore, we may expand the metric perturbation functions into a series of Legendre polynomials or their derivatives, depending on the parity of the corresponding perturbation function (see eg [44] or appendix A for details). Thus, for the odd parity perturbation function (r, θ), we have whereas for the even parity functions h, m and k: and the same holds for ζ(r, θ). Furthermore, one can show that the requirements of asymptotic flatness and regularity of the metric at the center of the star impose that only the l = 0 term of (33) survives, and similar arguments can be made for the second order perturbation functions, in which case only the l = 0, 2 terms are non vanishing [39]. Thus, the spacetime metric is reduced to (32), where (r, θ) = 1 (r) ≡ (r), h(r, θ) = h 0 (r) + h 2 (r)P 2 (cos θ), and so forth. Also, in the following, it will become useful to work with the shifted function ω(r) defined by Furthermore, we can make use of the residual gauge freedom of reparametrizations of the radial coordinate to set k 0 (r) = 0 in the expansion.
On the other hand, with these gauge choices, the stress-energy tensor of the system will be given by where, from the normalization condition for the fourvelocity, we have u µ = u t (1, 0, 0, Ω * ), and which, up to second order in , reads (38) To sum up, we have described the metric of spacetime associated to a slowly rotating perfect fluid star up to second order in the spin parameter. To do so, a perturbative expansion must be performed from a spherically symmetric, non-rotating metric in terms of a certain set of perturbation functions. We have chosen a particular coordinate system in which the surfaces of constant pressure coincide with those of constant radial coordinate, and written the stress-energy tensor of the rotating fluid in these coordinates. Therefore, we are now ready to obtain the Einstein equations for the system.

Interior Einstein equations
Thus, we now consider the Einstein equations (9), which can be written as E = 0 with E := G − 8πT , for the interior metric (32). The Einstein equations imply different equations for the perturbation functions, at each order in . Indeed, we may write where E (1) = ∂ E| =0 , and so forth, so that all terms in this expansion must vanish. For example, the zerothorder equations correspond to the TOV system of equations (10). At first order in , the only nontrivial equation is obtained from the first-order Einstein equation E (1) = 0, and corresponds to the (t, φ) component, E (1)φ t = 0, which yields The second order Einstein equations are given by E (2) = 0. As we have seen, the second order perturbation functions can be divided into two sectors, corresponding to the l = 0 and l = 2 terms in the Legendre expansion. Furthermore, these sectors appear uncoupled in the Einstein equations, so that we may separate these into different sets of equations for each sector. At quadratic order in , it will also be useful to consider, apart from the Einstein equations, the stress-energy tensor conservation equation. In particular, from the l = 0 sector of ∇ µ T (2)µ r = 0, one finds This equation can be integrated to yield an algebraic equation for ζ 0 in terms of h 0 and its initial condition, h 0 which is in a priori unknown and will be determined once the system is solved by matching with the exterior solutions. Further, from ∇ µ T (2)µ θ = 0, we have where we have used the zeroth order TOV equations. The last two equations are only valid inside the star since we are supposing p, ρ = 0. In particular, as they both correspond to algebraic instead of differential equations, the variable ζ 2 will not appear in the second-order system of differential equations since we can substitute directly by (42), and the same will happen to h 0 . Let us now obtain the differential equations for the rest of metric perturbation functions. The l = 0 contribution of E (2) t t = 0 gives On the other hand, the Einstein equation E (2)r r = 0 yields two independent equations which must be satisfied separately, namely, one for the l = 2 sector (obtained from the terms proportional to P 2 (cos θ)), and other for the l = 0 sector (from the terms independent of the Legendre polynomial). The equation for the l = 2 sector can be written whereas for the l = 0 sector we have Substitution of eq. (41) into (47) yields a differential equation for ζ 0 , which, together with eq. (43), constitutes a system of two ODEs independent of h 0 . Thus, once this system is solved, h 0 can be found algebraically using the integrated version of (41) up to an arbitrary constant.

Exterior equations and solutions
Following [28,39], we may take (23) as an ansatz for the metric of spacetime in the star exterior. We can indeed do this, since ζ(r, θ) is defined only inside the star, and taken to be constant outside. This means that the exterior metric in terms of r andr will be the same, where now the radial coordinate r goes from a finite value in the star surface R * , -corresponding to the star radius at zeroth-order-to infinity. Following the same steps as in the previous section, the Einstein equations in the exterior of the star can be obtained at each order in simply by setting ρ = p = 0 in Equations (40) and (44) to (46), with the exterior solution of the zeroth-order equations (TOV system) corresponding to the Schwarzschild solution by virtue of Birkhoff's theorem. Hence, to first order in , we have: ω ext = K 1 − K 2 /r 3 where K 1 and K 2 are two integration constants which can be related to the total spin velocity and angular momentum of the star. Indeed, at r → ∞, the metric function ω ext must approach the angular velocity of the star as measured by an static observer, so that K 1 = Ω * . On the other hand, we may calculate the conserved total angular momentum J of the star by integrating the angular momentum density current J µ = T µ ν k ν (φ) over a spacelike hypersurface Σ: from where it is straightforward to see that On the other hand, at second order in , the system given by eqs. (44) to (46) must be solved with vanishing ρ and p.
Using the expressions for the exterior solution of the first and zeroth-order metric functions, and imposing asymptotic flatness of the metric, one finds the analytic expressions [28,39]: being f (r) = (1 − 2M * /r), and where δM and A are integration constants. As we will see, δM corresponds to the correction of the gravitational mass, whereas A will be associated to the Love numbers.

Numerical solution for the interior and matching
Once we have obtained the system of differential equations for the metric functions, we need now the initial conditions in order to solve it. In this section we will explain how to obtain them and also how to solve the system numerically.
At this point there are no differences between how to solve the exact case and the mean-field case since the shooting method for the exact case is required only for the zeroth-order equations and those have already been solved. Thus we do know which value of the pressure in the center of the star corresponds to a given baryon number.
To obtain the initial conditions, as before, we expand our metric functions in powers of the radial coordinate and insert them in the differential equations to obtain the relations between the coefficients. In the zeroth-order (non-rotating) problem it is enough to expand until the zeroth order coefficient (in powers ofr), however when dealing with the second-order functions we need to expand them to the first nontrivial order (with nonvanishing coefficients). The reason is that the metric functions h 2 and k 2 vanish at the center of the star. Furthermore, for the next term of the expansions we find that they are equal and opposite, thus cancelling each other when substituting into their equations. This implies that: 1. We need a really good accuracy in the step of the numerical integration and 2. We cannot obtain the value of the first nontrivial coefficient of h 2 , in its expansion in powers ofr. To solve both problems, we follow [45] and start the integration in some small radius R (instead of r = 0) such that the expansions (63) are sufficiently accurate and the integration does not depend on the value of R . The expansions of the metric functions, with the nontrivial coefficients expressed in terms of the functions inr = 0 are where ω 0 , as p 0 , is an input parameter when it comes to solve the system. This parameter will determine the angular velocity of the star Ω * , as can be seen from the matching condition of ω with the exterior solution, ω ext , obtained in the previous section. This matching condition is simply given by imposing that the metric function ω and its first derivative are continuous throughout the star surface [28], i.e. ω(R * ) = ω ext (R * ), ω (R * ) = ω ext (R * ).
Therefore, in the rotating case, the stars are identified by a two-parameter family (ω 0 , p 0 ). The values of ρ 2 and N 2 are easily obtained from the EoS (ρ(p), n(p)), and h (2) is obtained from (17). The functions k 2 and m 2 are found to satisfy k 2 , around the center. As we have said, the odd powers inr of almost all the metric functions are null, however the definitions of M (11), m and ζ 0 (32) lead to the expansions given in (63). Now we start the integration with a non-zero, but still unknown, seed for the second-order functions h 2 and k 2 . To solve the unknown initial condition issue we will follow the approach given in [39,45]. First we must obtain a particular solution for h 2 and k 2 (h p , k p ) by solving the equations (46) and (45) for an arbitrary initial value (that must satisfy the regularity conditions given in (63)). Next, we obtain a homogeneus solution (h h , k h ) again for an arbitrary initial condition, using the same equations but with vanishing source terms. With these two functions we can construct the solution h 2 (r) = h p (r) + Bh h (r), k 2 (r) = k p (r) + Bk h (r). (65) In these expressions B is a constant that can be obtained by matching the functions h 2 and k 2 at the surface of the star with their corresponding exterior solutions. This matching condition is simply given by continuity of both functions at R * , i.e.
By doing this we are introducing the integration constant that appears in (54), hence we have an algebraic system of two equations that can be solved for A and B.
On the other hand, to solve the l = 0 sector of the second order system, we first solve the coupled ODEs for ζ 0 and m 0 as explained in the previous section, and then we obtain the solution for h 0 up to a constant h c 0 whose value is determined from the matching conditions where the constant term in (67a) is due to a nonvanishing energy density at the surface of the star, as pointed out in [43,46]. To obtain this constant term we can integrate (87) in the interval [R * − , R * + ] and take → 0. By doing this we have that all terms in the surface of the star vanish but the term dρ/dr, which is unbounded at R * and contributes with a constant term. From (67a) we can obtain the value of δM , which reads We would also like to remark a subtle detail concerning the second-order equations. When solving the TOV numerically, the metric function α is not fixed to its correct initial value since it does not affect the observables of the star. However, for solving the second order problem it is necessary to find the correct initial value of this function since the second order perturbation functions depend directly on α(0) and an incorrect value will affect the values of the quadrupole moment and gravitational mass correction of the star. This can be done by first solving the TOV system, finding the correct initial value of α using the matching condition (19) and then solving both zeroth and second order systems.

B. Global properties of compact stars
A key feature of the Hartle-Thorne perturbative formalism is that it allows us to obtain the values of these observable parameters from the solutions of the Einstein equations for the interior of the star at each order in the expansion parameter. Indeed, once these solutions have been obtained, they can be matched to the exterior solutions, from which observational parameters such as the quadrupole moment can be obtained systematically.
Take for example the moment of inertia I, which is defined as the quantity measuring how fast a star can spin given a fixed spin angular momentum J, and is given by To obtain the value of I for a given (interior) solution of the second-order Hartle-Thorne equations is straightforward: we simply obtain J from (49), by matching the exterior solution to the interior one at R * (64) and dividing by Ω * . It will be convenient also to define the dimensionless moment of inertia as On the other hand, the metric generated by an isolated, static gravitating body at a given point sufficiently far from the source may be written using a multipolar expansion in a system of Asymptotically Cartesian and Mass Centered coordinates [47,48], whose (0,0) component will be of the form where M is the gravitational mass of the star [49], and Q ij is the (traceless) quadrupolar tensor. The induced quadrupolar deformation of the star can be described in terms of the star's l = 2 sector perturbation functions in spherical coordinates. Indeed, defining x i = rn i (θ, φ), where n i is the unit three-vector in spherical coordinates, we may write: (where Y 2m are the l = 2 spherical harmonics). We find, in the case of an axially symmetric deformation, that the expansion (71) reduces to which defines the quadrupole moment of the metric, Q. Thus, we may perform an asymptotic expansion of the Hartle-Thorne perturbative solution for the exterior spacetime metric and identify the gravitational mass and quadrupole moment as the coefficients proportional to 2/r and the P 2 (cos θ)/r 3 term, respectively. Clearly, these quantities get corrections due to the star rotation.
Indeed, for example, the gravitational mass of the star, up to second order in , receives a correction which can be obtained from the expansion of the h ext 0 perturbation function. Furthermore, taking into account the asymptotic expansion for large r of h ext 2 and ω ext one finds that the spin-induced quadrupole moment of the star, up to second order in the spin parameter, is given by For later convenience we also define the dimensionless rotationally-induced quadrupole moment as Dropping the staticity assumption, nontrivial current multipole moments may appear in the expansion of the (0, j) components of the metric, the first term corresponding to a non-vanishing angular momentum. Finally, another interesting property that can be obtained from the solutions is the binding energy, which physically corresponds to the amount of energy that keeps all the particles (baryons) in the star from dispersing to infinity. It is defined as where M is the gravitational mass (in the static case, M = M * ) and M b is the baryon mass of the star. The binding energy so defined includes both the gravitational binding energy and the nuclear binding energy. However, we will be mostly interested in the gravitational contribution to the total binding energy, i.e. the gravitational binding energy, since it contains EoS-independent information about the mass distribution of the star [50]. The gravitational binding energy is defined as E g = M − M p , being M p the proper mass, given by the proper energymomentum density, P µ = T µν u ν , integrated on a spacelike hypersurface with volume form dS µ : In a stationary spacetime, this integral does not depend on the chosen hypersurface, so we may take dS µ = n µ d 3 S, where d 3 S = √ γd 3 x is the volume element of the spacelike hypersurfaces defined by t = const, γ is the determinant of the three-metric associated with these hypersurfaces and n µ = ∇ µ t/ (∇ ν t∇ ν t) is the corresponding normal vector, so that, for the static case, In the slowly rotating case, the perturbed proper mass M p ( ) will also get corrections. Expanding both γ and the product u µ n µ in powers of , we have, up to second order, M p ( ) = M * p + 2 δM p , where ω 2 e −α/2 dr.
(80) Hence, it is straightforward to obtain the second order perturbation to the gravitational binding energy E g ( ) = E g + 2 δE g , with δE g = δM − δM p .

IV. TIDALLY DEFORMED STARS AND LOVE
NUMBERS.
Until now we have studied the deformation of stars resulting from their own rotation. However, we can also study (non-rotating) stars which are deformed due to some external tidal force. Tidal forces are one of the principal signatures of the presence of a nontrivial gravitational field in spacetime. Such forces are responsible for relative acceleration among freely falling particles. This acceleration induces, on extended gravitating bodies, a field of strains that causes a deformation, which may be measured. By measuring the deformation response of a body to a tidal gravitational field, we may obtain information about the kind of matter that conforms the body, as well as its equation of state. In particular, in the case of binary systems involving neutron stars, it is very useful to analyze the deformation of the stars due to tidal effects, which may be measured from its gravitational wave spectrum previous to the merging.
On the other hand, as we have previously stated, a spherical body immersed in an external tidal field may deform due to tidal forces. Owing to this deformation, the metric in the exterior spacetime will develop a non trivial multipolar structure. To characterize the tidal field generated by a given source, consider an observer immersed in a tidal field generated by an external source. We may expand the metric of spacetime in a region surrounding the observer's worldline in Fermi normal coordinates, with the (0, 0) and (0, j) component of the metric given by [48,51], where E ij and B ij are the (quadrupolar) tidal multipole moments of electric and magnetic type, respectively. These two are related to the Riemann tensor through E ij = R i0j0 and B i j = 1 2 ijk R 0jkl [48]. The quadrupolar tidal moments are independent of the distance to the source, but may depend on the time coordinate if the source is not stationary. Now, instead of the worldline of an observer, we may consider the worldtube of an extended, gravitating body immersed in an external tidal field. We thus may be able to write the (0,0) component of the metric outside this body by combining both eqs. (71) and (81a), whereas the (0, j) component of the metric will be given by the combination of eqs. (77) and (81b), (83) Note that, by writing the metric as in eqs. (82) and (83), we are assuming that there exists a region of the exterior spacetime, called the buffer region, in which the expansions of eqs. (71), (77) and (81) converge simultaneously. This will be well justified in the limit in which the source of the external tidal field is very far away from the body that gets deformed and does not evolve rapidly with time. It can also be shown that, in this limit, the multipole moments appearing in (82) are defined unambiguously [52].

A. Electric quadrupolar Love number
Since we are considering that the body gets deformed due to the external tidal field, the quadrupole tensor Q ij will be a more or less complicated function of the tidal field E ij . However, working to linear order in the tidal moment, we define the (tidal) electric quadrupolar deformability λ t as Assuming that the terms with non-zero axial number m vanish, we may write (81a) in spherical coordinates as so that the tidal electric Love number can be obtained as the ratio λ t = −Q/E, where Q is the quadrupole moment of the star as defined in (73).
That the deformation of the star resulting from an external tidal field will be well described by its deformability λ is consistent with the assumption that the source of this external field is far from the body, since the tidal field will be weak and the linear approximation will be well justified.
The quadrupolar deformation of the star due to an external tidal field and to a slow rotation can be described by a similar spacetime metric (up to second order) [23,53], hence we can take advantage of the differential equations derived above to obtain the results for a tidally deformed star. Indeed, to describe a tidally deformed star, one introduces the metric perturbation h µν as in appendix A. By direct comparison between the metric (23) and (A6), it is straightforward to see that the l = 2 even perturbation functions H 2 , M 2 and K 2 in the tidally deformed case play a similar role as the functions h 2 , m 2 , k 2 in the slowly rotating case. Indeed, this can be seen by redefining these functions as In order to calculate the quadrupolar deformation of the metric due to an external gravitational field, the odd perturbations to the metric are not needed. Therefore, the metric functions of non-rotating tidally deformed stars can be directly obtained from eqs. (44) to (46) by imposing ω = 0. Then, these equations can be arranged into only one equation for h 2 , dρ dp e β + (α ) 2 h 2 .
(87) This is a second-order differential equation which can be solved as a first-order system as in [23] by defining H ≡ h 2 , while H is given by (87). To do this we need, again, to expand the function h 2 in powers ofr and introduce it in (87) to obtain the initial condition Once more, the value of h (2) can not be found from the limitr → 0 of the field equations. However, we will see that this is unimportant to find the correct value of the tidal deformability (see below), so that we can start the integration with an arbitrary value for h (2) 2 , as long as (88) are satisfied.
Once the internal solution has been found numerically, we can calculate the (tidal) Love number from the external solution of the metric after matching with the internal solution using the matching conditions where c 1,2 are constants that can be determined through the matching conditions in terms of H int (R * ), h int 2 (R * ). Studying the behavior of this function in the buffer zone we can extract the expression for Q and E in terms of these constants c 1 and c 2 in order to obtain the tidal electric deformability. Indeed, in the buffer zone and comparing with eqs. (73) and (85) we have and, defining the tidally induced, quadrupolar electric Love number k E 2 = 3 2 λt R 5 * , we can write where y = R * H ext (R * )/h ext 2 (R * ) and C = M * /R * is the compactness of the zeroth-order solution. From the definition of y it is clear that k (tid) 2 , and hence λ t , does not depend on the value of h (2) chosen in the numerical integration of the interior equation. Again, for later convenience we will define the adimensional tidal deformability asλ

B. Magnetic quadrupolar Love number
While electric-type Love numbers measure the induction of different multipole moments on a star due to an external gravitational field and can also be calculated in the Newtonian limit of general relativity, the current multi-pole moments induced by an external magnetictype tidal field have no analogue in Newtonian gravity, and thus the magnetic tidal Love numbers are a genuine prediction of general relativity. In the simplest (quadrupolar) case, the tidal magnetic deformability, in analogy with the electric case, measures the magnitude of the quadrupolar current S ij induced in the star by an external tidal field of magnetic type, B ij . At the linear level, the relation between both is Therefore, it is interesting to study the response of a NS under a magnetic-type external gravitational field, whose effects may be relevant for such compact objects.
To do so, we consider an axially-symmetric perturbation of the spherical metric. For the calculation of magnetictype Love numbers, only the odd metric perturbations are relevant. The magnetic Love number can be obtained by assuming a perturbation of the static metric of the form g( ) = g (0) + h odd where g (0) is the static spherically symmetric metric (8), here does not have to do anything with rotation, but will play the role of a bookkeeping parameter, and h odd is the odd-parity perturbation: h odd µν dx µ dx ν = 2V (r, θ)drdφ + 2ω(r, θ)dtdφ.
Notice that we have dropped the barred radial coordinate, since the star shape will not be altered by the odd metric perturbations. On the other hand, using the notation of appendix A, we define n i A = ∂ A n i . Then, we can transform the (0, j) components of the metric to spherical coordinates by g 0j → g 0A = rn j A g 0j , and expand into odd-parity vector harmonics (see appendix A), for instance, so that eq. (83) is transformed into In particular, for the simplest case of axially symmetric perturbations, we have, for instance, where · · · in eqs. (99) and (100) denotes the non-leading terms in the expansion at the buffer zone. Hence, the magnetic tidal deformability can be obtained as the ratio where the constants S and B will be determined from the buffer zone expansion of the odd-parity metric perturbation functions (97). To find these functions, the Einstein equations must be solved for the star interior, and matched to a suitable exterior solution, as we have previously done for the even parity case. However, in the case of odd-parity tidal perturbations, the energy-momentum tensor of the fluid will also get perturbed through a perturbation of the 4velocity, u µ ( ) = u µ + δu µ , working to first order in , and δu µ = g (0) µν δu ν + h odd µν u ν . Now, in principle, the four-velocity perturbations are independent of the metric perturbations, and the latter are only related with the former through the perturbative Einstein equations. However, there are two simple cases for which these perturbations are closely related to each other: the static case, in which the fluid remains static -with vanishing spatial four-velocity components i.e. δu µ = 0 -even when the metric perturbations are taken into account, and the irrotational fluid, which is based on the assumption that the fluid perturbations preserve the relativistic circulation theorem [55], and can be shown to be equivalent to the condition of a vorticity-free fluid, i.e. with vanishing vorticity four-vector ω α = 1 2 αβµν u β ∇ µ u ν = 0 which in turn implies δu µ = 0 [56], since the static initial configuration is trivially vorticity-free. The latter assumption is usually considered as more physically relevant, as the static fluid is only adequate for the non-physical case of timeindependent tidal perturbations [57]. Hence in the following we will only consider the case of an irrotational fluid and write the four-velocity perturbation as Substituting this expression into the stress-energy tensor perturbation, we may expand the Einstein equations as in (39) and solve the linearized equations E (1) = 0, which yield V 2 = 0 and the following ODE: for the metric perturbation function ω 2 . The above equation is numerically integrated in the star interior starting with an initial condition where, as in the electric case, the exact value of ω (3) 2 is undetermined from the equations but won't be needed for obtaining the corresponding Love number.
On the other hand, the exterior solution of eq. (104) can be written in terms of the hypergeometric function 2 F 1 (α, β, γ; x) as [58] ω ext 2 (r) =d 1 r 2M The expansion of (106) in powers of r and r −1 in the buffer region reads And, comparing with (100), we have so that the magnetic quadrupolar Love number k M 2 = 48 σt where now y = R * ω 2 (R * )/ω 2 (R * ), and once more it is clear from this expression that the initial condition for ω 2 does not enter in the expression for the magnetic Love number.

V. QUASI-UNIVERSAL RELATIONS
A. I-Love-Q In their original paper [28], Yagi and Yunes present a set of EoS-independent relations between the dimensionless moment of inertia, quadrupole moment and Love numbers of slowly rotating and tidally deformed compact stars, the so-called I-Love-Q relations. Soon after these relations where proposed, in [59] two possible reasons for these relations to exist were given. The first one relies on the fact that these relations depend mostly on the outer core (10 13 ≤ ρ ≤ 5 10 14 g/cm 3 ) of the NS, where all the EoS extracted from the experimental data of the nuclear physics (SLy, APR, WFF1, etc.) follow the same behavior. The second is related to the no-hair conjecture of black holes (BH) since the three parameters (Ī,λ andQ) must approach the limiting values of a BH for stars with large compactness.
In fig. 2 we show that these relations for the BPS (both exact and mean-field limits) also satisfied. We also show the data for the standard Skyrme crystal, generalized and hybrid EoS of [22] (which satisfy these relations as well) and the numerical fit for each of these relations obtained in [28] is plotted with a black line. Although somewhat expected, this result is remarkable at least for the case of exact BPS models, for which the I-Love-Q relations are satisfied even when they present a non-barotropic EoS which varies depending on the chosen potential. Furthermore, the relations are satisfied for these models in the exact and mean-field cases. As we will see, this will not be true anymore for other quasi-universal relations. This points out the universality and EoS independence of the I-Love-Q relations. obtained in [28].
Additionally, more quasi-universal relations have been found between electric, magnetic and higher multipole Love numbers [58]. For example, in 3 we show how there is as well an EoS independent relation between the (dimensionless) electric and magnetic quadrupolar tidal de-formabilities in all models considered. Apart from the I-Love-Q relations, there exists other set of relations between the moment of inertia, the Love numbers and the compactness of neutron stars that share some chatracteristics with the I-Love-Q but are accurate only up to ∼ 10%. These I-Love-C relations were approximately derived analytically in [60], as well as a possible explanation for these relations, in terms of the behavior of the energy density in the star interior. It turns out that these relations, as opposed to the I-Love-Q relations, are not universally satisfied for all the models we have considered. Indeed, from fig. 4 it can be seen that the relation between I, Q and C generally splits into two branches, corresponding to usual neutron stars and incompressible stars. This is consistent with the findings of [60]. However, we also find that, although the mean field version of the BPS models does lie in the incompressible star branch, the exactly solved cases behave quite differently. Whereas the behavior of the Partially-flat and 2χ-BPS models is better adjusted by the NS branch, the 4χ 2 −BPS model does not fit in neither branch. This behavior can be traced back to the radial dependence of the energy density in each model. Indeed, from fig. 6 we can see that the mean field approximation is not good in order to describe the low-density regime of neutron stars within the BPS Skyrme models in general, which translates into very different behaviors of the I-Love-C relations for these models. Indeed, it is clear from this figure that the MF approach overestimates the energy density of the stars in the outer regions.
Furthermore, the energy density profile for the different BPS models highly depends on the chosen potential. For example, while the θ-potential yields almost incompressible stars, the 2χ-potential curve can be well approx- imated by a quadratic function. This quadratic behavior is in fact expected for realistic neutron stars, whilst the behavior of the density profile for the 4χ 2 -model is actually more similar to that of white dwarfs [60]. Indeed, as we have seen in figs. 4 and 5, the (exact) 2χ-BPS model I-Love-C relations are very close to the NS fit from [60]. Finally, it is interesting that the curve for the P.F. potential does not fit to a quadratic curve, yet the I-Love-C relations within this model are still satisfied. Skyrme stars for different models. We include the nuclear based EoS BCPM and AP4, as well as the quadratic curve 1 − r 2 .

C. Gravitational binding energy relations
A different set of quasi-universal relations involving the static gravitational binding energy and other global properties of neutron star solutions have been recently proposed in [50]. For instance, we show in fig. 7 the universal behavior of the static gravitational binding energy normalized to the TOV mass and plotted against the adimensional moment of inertia. From the error plot one can see that all models follow the same universal behavior with a deviation of 5% (but the exact 4χ 2 BPS model, in which case the error is as high as ten percent) with respect to the numerical fit obtained in [50]. Further, the rotation of the star has measurable effects both in the gravitational and proper mass of the star. Indeed, as we have seen, the gravitational mass of the star receives a correction δM , the dimensionless version of which, δM = δM × M 3 /J 2 , was also shown in [46] to satisfy a universal relation when plotted against the (dimensionless) tidal deformability. We show this relation in fig. 8, together with the numerical fit of [46] obtained for the region λ t < 10 3 , at which the deviation for all models is less than ten percent.
On the other hand, the gravitational binding energy will also get a second order correction, namely, δE g . Remarkably, as opposing to its zeroth-order counterpart, the correction to the gravitational binding energy does not seem to follow a simple, quasi-universal relation. Since the correction to the gravitational mass indeed does follow a relation as shown in fig. 8, the non-universal nature of δE g can be traced back to the correction to the proper mass δM p , which involves an integral over the star, see (80). In fig. 9 we show the behavior of universal behavior may be inferred, which corresponds to the numerical fit we have obtained. However, this behavior has not the same universality as others previousy analyzed, as the deviation can grow up to 30% for realistic masses.

VI. DEFORMABILITY CONSTRAINTS FROM OBSERVATIONS
In addition of constituting an outstanding experimental confirmation of the validity of General Relativity, the direct observation of gravitational waves can be used to place direct constraints on the neutron star EoS (see [61][62][63] for a recent review of nuclear EoS constraints from GW observations). Indeed, the waveform produced by the coalescence of two realistic extended bodies deviates significantly from a point-particle waveform and thus this difference can be observed with Advanced LIGO. The degree of the deviation, in the case of binary neutron star mergers, depends on the underlying EoS. Although the magnitude of the deviation is strongest at later times in the inspiral and during the merger, Flanagan and Hinderer found that the early phase of the inspiral depends mostly on the tidal Love number of the neutron stars, introducing a phase shift with respect to the point-particle waveform [53]. However, the individual Love numbers for each component of the merger cannot be separately distinguished in the observed gravitational waveform. Instead, what can be sharply measured is the so-called effective tidal deformability,Λ a mass-weighted average of the dimensionless deformabilitiesλ 1 andλ 2 of both components (with masses m 1 and m 2 ), given bỹ (111) Similarly, the two component masses are not measured directly. Instead, it is the chirp mass, where q = m 1 /m 2 is the mass ratio, what can actually be constrained. In the case of the GW170817 event, the chirp mass was measured to be 1.188 +0.004 −0.002 at the 90% confidence level. Moreover, within the same confidence level, the mass ratio was found to be in the range 0.7 − 1, and the effective tidal deformability to be smaller than 800 [64].
Such such measurements of the NS properties can be used to further reduce the set of Skyrme models able to reproduce physically realistic NS solutions and impose some constraints on the possible values of the free parameters of these models. Indeed, once the equations for the tidally deformed stars of section IV are solved for a specific model, we may obtain the dimensionless tidal deformability of stars described by this model as a function of their TOV mass, so thatΛ may be seen as a function of both m 1 and m 2 , or, equivalently, of M c and q. On the other hand, since the chirp mass of the binary progenitor of GW170817 is well measured, for any given EOS the effective deformability reduces to a simple EOS-dependent function of the mass ratio.
In fig. 10 we show the effective tidal deformability as a function of the mass ratio for a chirp mass of 1.19 M for different Skyrme models, together with the constraints from the GW170817 event. It is clear from this figure that, as we have already argued, the mean field approximation is not suitable for describing the low energy region of BPS stars, which makes the most relevant contribution to the deformability of the stars. Indeed, this approximation overestimates the values of effective deformability by at least a factor of ∼ 2.
On the other hand, we find that both the generalized and the hybrid EoS provide an excellent description of the tidal deformability [22].
In addition, we see again that the exact BPS Skyrme models present very different behaviors -in this case, different values ofΛ-depending on the chosen potential. For example, the contribution from the Θ potential is clearly too big as compared with the GW observation, which sets an upper value ofΛ ≤ 800, whereas the 2χ and PF potentials yield very largeλ, near the upper value, although still allowed by the bound).  Furthermore, it is likely that additional observations of gravitational waves from binary NS mergers will further constrain the tidal deformability of these compact stars. In particular, some recently observed GW events [65,66] strongly suggest that highly massive NS and compact objects within the NS-Black Hole mass gap (around 2.5 M ) could exist.
However, it is difficult to distinguish between an extremely massive neutron star and a small black hole from the GW waveform alone with first generation GW detectors, since the tidal deformability and quadrupole moment of such massive stars is usually very low due to their high compactness, and almost no realistic EoS is able to produce stars with such big mass.
In fig. 11, we show the dimensionless moment of inertia, quadrupolar moment and tidal deformability of all the Skyrme models as well as for the AP4 and BCPM EoS. These plots show that not only high mass neutron star solutions can be found for any BPS Skyrme model as well as for the generalized and the hybrid EoS. We also find that, depending on the potential, these parameters can acquire sufficiently high values to be able to be measured by current generation GW observatories. Therefore, we conclude that if the tidal deformability of a mass-gap compact object were measured to be non-zero, it is very likely that its EoS will be well approximated by a BPS Skyrme model and, in particular, by the hybrid model which approaches the BPS behavior at high density.

VII. CONCLUSIONS
In this work, we have solved the Einstein equations using the Hartle-Thorne perturbative formalism to find slowly rotating NS solutions with nuclear and Skyrme model-based EoS. Moreover, we have presented perturbative solutions to the Einstein-BPS Skyrme system describing slowly rotating and tidally deformed, selfgravitating solitons which can also be considered as idealized models for Neutron Stars. For all these models, we have computed different global properties of the corresponding star configurations, such as the moments of inertia, quadrupole moments, gravitational masses or binding energies, and checked whether or not all the models satisfy some (quasi-)universal relations previously proposed in the literature. As we have found, the I-Love-Q relations presented in [28] are satisfied up to a ∼ 2% error, even for the exact, non-barotropic BPS Skyrme models, which reaffirms the universality of these relations. Other relations involving the second-order correction to the gravitational mass (including the correction proposed in [46]) and those involving the (gravitational) binding energy are also quite well satisfied for all models at hand.
On the other hand, we have found that while the I-Love-C quasi-universal relations still hold for the meanfield BPS and Skyrme-based EoS, these relations break up for the exact BPS Skyrme models. This fact, as argued, can be traced back to the behavior of the energy density profiles of the solutions for such models, which strongly depends on the particular potential chosen due to the non-barotropic nature of these models. This finding is consistent with the explanation given in [59] about the difference in nature between these relations and the I-Love-Q.
The extension of previous works on Skyrmion stars to include the effects of small rotations and tidal deformations allows to enlarge the set of observable quantities that can be compared to actual measurements. Owing to the increasing number of observed GW events in recent years, those observables that can be inferred from the waveform of a GW produced at a binary NS coalescence are of particular interest. An example is the (effective) tidal deformability of the binary, which, together with the quasi-universal relations, allows to constrain the EoS of strongly interacting matter in the extremely high density regime. In this paper, we have shown that, within the Skyrme model, these universal relations still hold and that the current experimental bounds on NS deformabilities can be well accounted for. Furthermore, a remarkable property of the solutions based on generalized Skyrme models is that very high masses (of approximately 2.5M ) can be reached even for not too large energy densities at the center of the stars. In other words, such massive stars can be produced from mesonic degrees of freedom alone without the need of additional degrees of freedom such as unconfined quarks. This is consistent with the assumption that the Skyrme model is a valid approximation for the description of matter at the core of a NS, which, if true, implies that the pressure and density reached at NS cores are still far from the energy density regimes in which perturbative QCD becomes relevant.
Although some recent GW events can be seen as possible evidence that such massive stars may exist, additional observations are required to further clarify the detailed properties of massive NS cores.
We conclude by summarizing our main results: i) we find that the hybrid EoS, where the EoS of the generalized Skyrme model is complemented by a standard nuclear physics EoS for low densities, is compatible with all observational constraints, both for static NS observables and for observables related to slowly rotating and/or tidally deformed NS, and ii) we verify (quasi-)universal relations, like I-Love-Q, for a broad range of models, based on the minimal Skyrme model, the BPS Skyrme model with a variety of potentials, the generalized Skyrme model, and the hybrid model, respectively. In particular, the BPS Skyrme model also allows for an exact field-theoretic treatment (beyond mean-field theory), because it represents a nonbarotropic perfect fluid. These results contribute to a deeper understanding of the range of validity of these relations, because we investigate them for qualitatively different models not considered previously. In addition, our investigation also provides a better insight into the role played by each component of the full, hybrid Skyrme model. the problem of determining their explicit expression by solving the perturbed Einstein equations. A useful gauge choice is the so-called Regge-Wheeler gauge [44], in which j lm a = G lm = h lm 2 = 0. Furthermore, for stationary, axially symmetric perturbations, we may discard the φ−dependence of the harmonics, and the only nonvanishing contribution will be that of the m = 0 terms. Hence, the general expression for an stationary, axially symmetric metric perturbation with these gauge choices is given by h µν = h even µν + h odd µν , where: where Y l0 = P l (cos θ) and X l0 φ = sin θ∂ θ P l (cos θ), and summation over l is implied.