Slowly rotating topological neutron stars -- universal relations and epicyclic frequencies

In the modern era of abundant X-ray detections and the increasing momentum of gravitational waves astronomy, tests of General Relativity in strong field regime become increasingly feasible and their importance for probing gravity cannot be understated. To this end, we study the characteristics of slowly rotating topological neutron stars in the tensor-multi-scalar theories of gravity following the static study of this new type of compact objects by two of the authors. We explore the moment of inertia and verify that universal relations known from General Relativity hold for this new class of compact objects. Furthermore, we study the properties of their innermost stable circular orbits (ISCO) and the epicyclic frequencies due to the latter's hinted link to observational quantities such as quasi-periodic X-ray spectrum features.


I. INTRODUCTION
Some of the most promising extensions of General Relativity (GR) are based on the existence of additional fields which may contribute to Gravity either as additional energy-momentum sources or even through direct coupling to it. While weak-field tests place strong restrictions to such extensions, these fields may have a strong impact on the properties of compact objects such as neutron stars and black holes, which places these objects among the best "laboratories" for experimental tests of the strong field regime of gravity [1][2][3]. The Tensor-Multi-Scalar Theories (TMST) are among the most promising and natural such extensions of General Relativity which are mathematically self-consistent and pass all known experimental and observational constraints [4,5]. In this paper we explore some of the most important observational properties of slowly rotating stars in a certain class of TMST. A new type of neutron stars called topological neutron star was shown to exist in their framework [6] (see also [7] for the case of scalarization in such theories). Following the study of the static structure of this new type of compact objects, we further consider their properties under slow rotation in order to extract additional observational signatures.
One of the most important quantities which can be obtained under rotation is the moment of inertia, which can be computed to first order in the star's angular velocity. Due to the largely uncertain equation of state (EOS) for extremely dense matter, it is beneficial to construct universal relations which provide links between the parameters of neutron stars independently of the EOS. We will focus on two types of universal relations connecting the suitably normalized moment of inertia and the stellar compactness. These relations were examined in the GR case in [8,9]. Later they were generalized to the case of f (R) and scalar-tensor theories with a massive scalar field [10][11][12]. We follow their methodology and compare the topological neutron star results with the GR case for several EOS and at several theory parameters for two coupling functions in order to discover possible observational traits of the former.
The modern and future X-ray observatories such as NICER, LOFT and SKA [13,15] can give us invaluable data for neutron stars possessing less compact companions such as white dwarfs or main sequence stars, due to the high likelihood of accretion in such systems. One of the most important parameters that characterize such an accretion is the innermost stable circular orbit (ISCO) since it determines the boundary region where the accreting matter can no longer orbit the compact object under gravity alone. The epicyclic frequencies for a given orbital radius are the characteristic frequencies of oscillation for particles on stable circular orbits undergoing some small perturbations. There are many theories as to the origin of some quasi-periodic oscillations (QPOs) observed in the X-ray light curves of accreting compact objects, but there is no consensus which is the correct one (see [43] for a comprehensive review). The oscillations themselves are in the order of a few hundred Hz to kHz and most models are based either on orbital and epicyclic motion of matter [16,24] or oscillations and instabilities in an accretion disc [25,29] both of which are expected to occur close to the ISCO radius. Nevertheless, almost all QPOs models are related in some way or another to the ISCO radius and the epicyclic frequencies.
Since ISCO is located close (or on) the surface of neutron stars, the QPOs are a promising opportunity to probe the strong regime of gravity and its potential modifications. The QPOs in scalar-tensor type alternative theories of gravity were considered in [30]- [33]. In the present pa-per we will address this problem in the context of TMST and topological neutron stars.
We start with the standard static spherical metric, perturbing it to first order in the angular velocity Ω and solving the resulting system of reduced field equations. The final goal is to extract the moment of inertia for the topological neutron stars and compare two universal (independent of the EOS) relations, known to be valid in GR [8,9] and some alternative theories of gravity [10,12]. Furthermore, we compute the radius of the innermost stable circular orbit (ISCO), the epicyclic frequencies and their dependence on topological neutron star mass and angular velocity for the stable branches of solutions for one of the coupling functions considered.
Inline with the previous work, we consider a gravitational interaction mediated by the spacetime metric g µν and N scalar fields ϕ a with values in a coordinate patch of an N-dimensional Riemannian (target) manifold E N with a positive-definite metric γ ab (ϕ) defined on it. The Einstein frame action of the theory is given by where G * is the bare gravitational constant, R and ∇ µ are the Ricci scalar and the covariant derivative with respect to the Einstein frame metric g µν , and V (ϕ) ≥ 0 is the scalar fields potential. The matter fields are coupled only to the physical Jordan frame metricg µν = A 2 (ϕ)g µν in order for the theory to satisfy the weak equivalence principle. The conformal factor A 2 (ϕ), the target space metric γ ab and the scalar fields potential V (ϕ) specify the TMST.
The variation of (1) with respect to the Einstein frame metric components and the scalar fields gives the field equations of the theory in the Einstein frame as follows where T µν is the Einstein frame energy-momentum tensor of the matter and γ a bc (ϕ) are the Christoffel symbols with respect to the target space metric γ ab , following closely the notation in [6]. The conservation law for the energymomentum tensor obtained from the contracted Bianchi identities and the field equations is where once again the physical energy-momentum tensor in the Jordan frameT µν and the Einstein frame one in (3) are related through the conformal factor as T µν = A 2 (ϕ)T µν . The matter fields are described by a perfect fluid and by virtue of this relation, the corresponding energy density, pressure and 4-velocity transformations between the two frames are given by ε = A 4 (ϕ)ε, p = A 4 (ϕ)p and u µ = A −1 (ϕ)ũ µ .

II. STRUCTURE EQUATIONS FOR THE NEUTRON STARS AND SETTING THE THEORY PARAMETERS
We perturb the static spherical solution by a term scaling as O(Ω) in the Einstein frame where Γ, Λ and ω depend on the radial coordinate r only. Using the standard Hartle procedure [36,37], one can obtain that at spatial infinity ω tends towards where J is the angular momentum of the star. Following the previous work on topological neutron stars, we consider the target space for the TMST to be the round three-dimensional sphere S 3 γ ab dϕ a dϕ b = a 2 dχ 2 + sin 2 χ(dΘ 2 + sin 2 ΘdΦ 2 ) , (6) defined by a radius a > 0 of S 3 and the standard angular coordinates χ, Θ, Φ on S 3 . The motivation behind this choice is the fact that S 3 is among the simplest target spaces that allow the existence of spherically symmetric topological neutron stars. The assumption that the χ field depends only on the radial coordinate (χ = χ(r)) and that the Θ and Φ fields depend only on the corresponding angular coordinates (Θ = Θ(θ) and Φ = Φ(φ)) leads to unique solutions of the latter two, compatible with the spherical symmetry of the structure equations. These solutions are Θ = θ and Φ = φ as shown in [6].

A. Background structure equations
Under these simplifying assumptions, one can readily derive the structure equations for the static configuration. Using the field equations (2) and the conservation law (3), the reduced field equations take the form: 2 where the prime denotes differentiation with respect to the radial coordinate r. Naturally, these equations must be complemented by an appropriate equation of state (EOS) which provides the required dependence between the physical pressure and energy density (p =p(ε)) in order to complete the system. Throughout the work we use a piecewise polytropic approximation of several realistic nuclear matter equations of state [39]. These are the APR4 and Sly, that fit very well to the current observational constraints, as well as several higher/lower stiffness and maximum mass equations of state (MS1, MPA1, APR2, H4) in order to check the universality of the relations. Apart from the standard boundary conditions for the metric functions derived from asymptotic flatness Γ(∞) = Λ(∞) = 0 and regularity at the center Λ(0) = 0, one further requires such regularity from the scalar field equations. This leads to the conditions χ(∞) = kπ and χ(0) = nπ with both k and n integer numbers k, n ∈ Z. The former condition can be set to zero without loss of generality (k = 0). Since χ(∞) = 0, the extension of the map φ : Σ → S 3 is topologically equivalent to the map φ : S 3 → S 3 , and it can be shown that n is its degree [6]. Thus the solutions with n = 0 are topologically nontrivial.
We focus on two different coupling functions characterized by a single dimensionless parameter β. The first one is that originally used in [6] given by while the second one is a monotonic function of χ Solutions of the static topological neutron stars are obtained by setting the appropriate central conditions for Λ, ε and χ and performing a shooting method to determine (dχ/dr)| 0 and Γ(0) from the asymptotics of χ and Γ. For the numerical integration of the differential equations we are using the adaptive Dormand-Prince embedded error estimation method [40]. As shown in the previous work, solutions exist only for certain theory parameters. One of the main constraints is the size of the target space a with topological neutron stars emerging only for small values of it. The radius of the target space throughout this work has been set to a 2 = 10 −3 . There also exists a restricted range of β values at fixed a which depends on the coupling function. Extending the results obtained previously, several different values of β have been used for the function (11). These include β = 0.08 which was already employed in [6,41] as well as β = 0.04, β = −0.04 and β = −0.08. In the case of the coupling function (12), solutions were restricted to an even smaller interval of values for β and there is strong evidence that no stable configurations exist for β > 0 in that case.
The mass-radius relations for topological neutron stars is shown in Fig. 1 for the non-monotonic coupling function (11). In this case stable solutions exist both for n = 1 and n = 2 topological charges (excluding the case β = 0.04 where the n = 2 solutions appear after the branch's maximum mass). Fig. 2 on the other hand presents the mass-radius relations for the monotonic coupling function (12) where stable solutions exist only for n = 1 topological charge and in a much smaller range of β (around β = −0.01). Since we are working in slow rotation approximation keeping only linear terms with respect to the angular velocity Ω, the mass and the radius remain unchanged with respect to the nonrotating case. The results depicted in these two figures are for the stable branch of topological neutron star solutions. It was shown in [6] that other branches of solutions exist as well for a fixed coupling function and values of n and β, but all of them are unstable against radial perturbations [41] and we will not comment them further.
It is easy to notice that in general the n = 1 branches are more massive than their n = 2 counterparts (when the latter exist). Considering Fig.1, the overall effect of negative β values is to increase the masses of the solutions with the same radius (making them more compact) while the effect of positive β values is the opposite. The actual shape of the mass-radius relation, however, is also heavily influenced by the value of the parameter and these effects are not easily predictable or monotonic in nature.
In most cases the "end" of a branch (i.e. the highest pressure for which solutions were found) has been chosen arbitrarily after the corresponding sequence's maximum mass is reached. However, in some cases the branch of solutions not only appears at some minimal central pressure but also disappears at another maximum central pressure. This maximum pressure where solutions disappear The only non-trivial equation to first order in Ω is that for the Ricci tensor components R 03 = R 30 , which after taking into account the static equations (7-10) as equalities, simplifies to where we have definedω = Ω − ω. The central values of the functions areω ′ (0) = 0 andω(0) =ω c with the latter determined through a shooting method in order to obtain the desired angular velocity of the star Ω. Note that Ω is actually the same in both Einstein and Jordan frames. An important characteristic for compact objects is their moment of inertia. Having integrated (13), one can use its asymptotic form to extract the star's angular momentum J = IΩ and thus obtain the moment of inertia along the axis of rotation. Alternatively, one can integrate the two sides of (13) leading to a conserved quantity proportional to the angular momentum. After a division by Ω, the asymptotic and integral definitions are given by whereω = ω/Ω = 1−ω/Ω is the reduced angular velocity. Both definitions are implemented as part of the numerical integration to further verify the fidelity of the results and the obtained values differ from each-other with relative error consistently on the order of 10 −11 for the various branches. Section III outlines our results for the moment of inertia and its normalization conditions which lead to universal relations.

C. ISCO, orbital and epicyclic frequencies
For a general stationary axially-symmetric metric with g µν = g µν (r, θ) and a line element of the form ds 2 = g tt dt 2 +g rr dr 2 +g θθ dθ 2 +2g tφ dtdθ+g φφ dφ 2 , (16) the ISCO and epicyclic frequencies can be found by analysing the orbital motion for massive test particles. It is easy to show that there are two constants of motion generated by the timelike and axial killing vectors E = −u t and L = u φ . Raising these we obtain for the contravariant components of a particle's 4-velocity The 4-velocity's normalization condition for timelike observers g µν u µ u ν = −1, can be rewritten in the form where we have defined the potential using the proper orbital angular momentum l ≡ L/E. In the equatorial plane (θ = π/2) this result is further reduced to an effective 1D problemṙ 2 = V (r) with effective potential For some fixed E and L, the stable circular orbit at some coordinate radius r 0 is given by the conditions V (r 0 ) = V ′ (r 0 ) = 0 and V ′′ (r 0 ) > 0, while the radius of the ISCO is given by the marginal stability condition V ′′ (r 0 ) = 0. The orbital angular velocity of massive particles in the geometry can be found from the geodesic equations written in their Lagrange-Euler form The radial component of these equations (µ = 1) reads which is transformed into a quadratic algebraic equation for the orbital angular velocity Ω p = dφ/dr = u φ /u t . The two solutions are easily found to be corresponding to prograde and retrograde orbits with respect to the star's rotation. The follow-up computations are considered for the prograde case as the retrogade one's angular velocity is always lower and stability is lost further out from the star. For a particle on a stable circular orbit, small radial or angular perturbations will cause periodic oscillations about the potential's minimum (with respect to r or θ respectively). These frequencies are known as the radial and vertical epicyclic frequencies and can be found by investigating time-dependent perturbations of a stable equatorial circular orbit in the form Inserting (25) into the equation of motion (19) and assuming δr ∝ e iωrt , δθ ∝ e iω θ t yields for the radial and vertical angular epicyclic frequencies (ω i = 2πν i ). Given the interpretation of these frequencies and the proportionality between ω r and ∂ 2 r U (r 0 , π/2), it is clear that the radial epicyclic frequency must be zero at the ISCO radius, real for r > r ISCO and imaginary for r < r ISCO . Furthermore, the vertical epicyclic frequency is equal to the orbital one ω θ = Ω p for the static case (Ω = 0).
Of course, realistic neutron stars possess magnetic fields which can be far from negligible on the dynamics of charged particles and recent studies show that those can be more complicated than originally anticipated [42]. Nevertheless, since plasma is electrically neutral on the scale of several Debye radii, each "clump" of plasma would be held together by electrostatic stresses between the constituting particles so even though individual charged particles' trajectories may differ due to the magnetic field, the overall plasma motion is well described by the considered dynamics.
The same general methodology is also followed for TMST with one major difference -all derivatives and quantities in the latter case must be computed using the physical (Jordan) frame metric componentsg µν . This leads to a noteworthy behavior in the orbital angular velocity Ω p as given by (24) which can turn imaginary (i.e. no stable circular orbits exist) in regions outside the star, before the appearance of a true ISCO. This can be traced to the negative term under the square root of eq. (24). In GR the derivative ∂ r g tt is always negative outside of the star (the metric function g tt is monotonically decreasing). However, the sign of the corresponding quantity in the generalized theory ∂ r (A 2 (χ)g tt ) = ∂ rgtt is now determined not only by g tt but also by the coupling function A(χ). It turns out that in the case of the non-monotonic coupling function (11) this leads to regions outside the ISCO where stable circular orbits do not exist and where the accretion disc may be split. This effect presents considerable interest in terms of potential observational traits but requires a more detailed study. For that reason it will be explored in a separate work. In section IV therefore, we outline our results for the radius of ISCO, the orbital frequency at ISCO as well as the maximum value of the radial epicyclic frequency for the n = 1 stable topological neutron stars of the monotonic coupling function (12) at three different β values and three different rotation rates, comparing them to the results for GR.

III. MOMENT OF INERTIA AND UNIVERSAL RELATIONS RESULTS
In this section we present the results obtained for the moment of inertia as well as our study of two relations known to be universal (independent of EOS). These are relations between the dimensionless compactness M/R and two normalizations of the moment of in-ertiaĨ = I/M R 2 andĪ = I/M 3 . The functional forms were fitted with a polynomial though least squares. Several EOS were used and the error of those fits was evaluated in order to verify their validity. As discussed, we focus only on the astrophysically relevant stable branch of topological neutron stars without commenting on the other unstable branches discovered in [6]. A comparison is made between GR and the branches with n = 1 and n = 2 topological charges (where present) using APR4, Sly, MS1, MPA1, APR2 and H4 equations of state for several different theory parameters.
Before considering these results for the different EOS, the following Figs. 3 and 4 show the moment of inertia as a function of the neutron star masses. These results are obtained for APR4 EOS and several combinations of n and β in order to provide some intuition on the effect of the coupling function and the β value. As a matter of fact these are the same branches of solutions depicted in Fig. 1 and Fig. 2. way (for the non-monotonic coupling function (11)). A complete lack of maximum of I(M ) is observed in the case of the monotonic coupling function. It is important to note, that up to first order of the angular velocity, the moment of inertia for a star is independent of Ω as no deformations are taken into consideration. We have confirmed that the value for I is indeed numerically independent of the different Ω values set in our integrator. Let us now turn to the universal relations. While in the case of mass-radius relations the plots are qualitatively different in shape from those of the GR case (see Figs. 1 and 2), the normalized moment of inertia for all theory parameters follows a very similar dependency. The functional shape is very similar with an overall translation depending on the values of the free parameters (with the exception of β = 0.08, n = 1 for the non-monotonic function).
The fourth order polynomial with zero second and third order terms is now a standard fit for the normalized moment of inertiaĨ = I/M R 2 as a function of the compactness M/R [8,9] We will first give plots of the universal relations polynomial fits for different combinations of the free parameters in order to gain intuition on the possible deviations from GR and only afterwards we will comment on the deviation from EOS university. Let us just point out, that this universality is more or less preserved for topological neutron stars as well. The olive curve is that obtained for the monotonic coupling function (12) while the remaining four curves show the qualitative difference in results for positive and negative beta in the case of non-monotonic coupling (11). β = −0.08 and β = 0.08. These fits are obtained using an equal number of points for each of the 6 EOS (APR4, SLy, MS1, MPA1, APR2 and H4) and fitting the functional form (28) through least squares. In Table I the numerical values of the fitting coefficients, as well as the fitting error, for some representative combinations of the parameters are displayed. It is evident both from the table and the plot that there is very little difference between the polynomial fits for GR and the TMST for the monotonic coupling function A mon (χ) given by (12). This is not the case, though, for the non-monotonic function A nmon (χ) given by eq. (11) where significant differences are observed. In this case the normalized moment of inertia is significantly higher or lower as compared to GR both for n = 1 and n = 2. For example in the β = −0.08 case the polynomial fit is consistently around 20 % higher as compared to GR. The n = 1 fit for the positive β = 0.08 has a slightly different shape from GR particularly at lower compactness but would be hard to discern within experiments available today. The n = 2 solutions of the same value β = 0.08, however, lead to values ofĨ consistently lower by about 7%.
The second alternative normalization of the moment of inertia we consider, namelyĪ = I/M 3 , is fitted with a polynomial function of the form [9] The numerical values of the fitting coefficients, as well as the fitting error, can be found once again in Table I  The olive curve is that obtained for the monotonic coupling function (12) while the remaining four curves show the qualitative difference in results for positive and negative beta in the case of non-monotonic coupling (11).
certain representative cases. Fig. 6 compares theĪ fit obtained using the above formula for the same values of the parameters and equations of state as in Fig. 5. Once again, the curve for the coupling function (12) is virtually indistinguishable from the GR one and the most significant deviation is observed for the non-monotonic coupling (11) with n = 1 topological charge at β = −0.08 (about 20 % higher) and with n = 2 topological charge at β = 0.08 (about 7 % lower).
The following two Figs. 7 and 8 show the actual data points for each of the six EOS considered as well as the fits obtained from them. Only the value of β = −0.08 of the non-monotonic coupling function (11) has been shown for the purpose of better visualization. This value of β has been chosen since it presents the highest difference as compared to GR. The n = 1 and n = 2 solutions are depicted with pluses and crosses and follow the colors of the GR solutions for the different EOS. The highest relative error of the fits is displayed in the bottom of the figure and it is no greater than 6% (highest for the H4 EOS). Based on the relative error of the fits, we can conclude that the EOS universality ofĨ as a function of the stellar compactness is fulfilled in TMST at least as well as in GR, and the differences between GR and TMST fits are significant at least for some of the β and n values.
Similarly, as seen on the following Fig. 8, the error between the actual results and the fit for the second normal-izationĪ does not exceed 6% with the highest maximum of the error once again visible for H4 EOS which confirms the universality ofĪ as well for β = −0.08. The fits for the remaining β values and the monotonic coupling function have a relative error of the same order, confirming that the topological neutron stars of both charges n = 1 and n = 2 indeed satisfy the proposed universal relations and in some cases the predicted difference in the obtained fit is significant as compared to GR.
Theory parametersã 0ã1ã2Ĩχ 2ā 1ā2ā3ā4Īχ 2   (28) and (29). We have chosen some representative combinations of parameters for both coupling functions, where A mon denotes the monotonic coupling (12) and A nmon -the non-monotonic coupling (11). Each set of values was obtained by taking an equal number of points for each EOS (in order to obtain the same weight) for mass range starting at 1 M ⊙ and ending at the maximum mass for the corresponding sequence.

General
It is evident from the values ofĨ χ 2 andĪ χ 2 displayed in Table I that the topological neutron stars exhibit a good degree of EOS universality with respect to the two normalizations of the moment of inertiaĨ = I/M R 2 and I = I/M 3 and therefore, these relations are well described by the standard fits (28) and (29). The obtained fits for the monotonic coupling function (12) are practically indistinguishable from GR. The non-monotonic coupling function (11), however, leads to quantitatively significant differences between GR, the n = 1 and the n = 2 topological neutron stars independently of the EOS.

IV. ISCO, ORBITAL AND EPICYCLIC FREQUENCIES RESULTS
In this section we outline our results for the ISCO as well as the orbital and epicyclic frequencies. The results are obtained for static as well as f = 80 Hz and f = 160 Hz rotating neutron stars. Both of these values for the frequency place the stars at the slowly rotating regime. Unlike the moment of inertia, which is not influenced to the first order in Ω, the innermost stable circular orbit (ISCO) as well as the orbital ν p and epicyclic frequencies Stars of n = 1 and n = 2 respectively with the same color scheme for the EOS. Solid, dashed and dash-dotted lines show the corresponding fits following (29).
ν r , ν θ are modified based on the star's rotation. Since ISCO is the limiting radius where particles can orbit stably, it has important applications to determining the inner edge of accretion discs around compact objects. We should note that all the results presented in this section are the corresponding physical Jordan frame quantities. As discussed in section II C, the ISCO computation in the case of a non-monotonic coupling function such as (11) is not trivial and can lead to a splitting of the accretion disc outside of the star. For the present section, therefore, we have considered only the monotonic coupling function (12) while the more complicated case of non-monotonic A(χ) requires much more profound investigations and will be discussed in a future publication. Fig. 9 displays the static ISCO radius (in the physical Jordan frame) for the stable branch of topological neutron stars with n = 1 at different β parameters for the APR4 EOS and compares them to GR. Whenever the ISCO radius is less than the stars' radius, the latter is used (thus the visible discontinuity of the derivatives around 1.25 M ⊙ ). Stable solution for n = 2 do not exist for this coupling function and as previously discussed, the maximum mass that the neutron stars with n = 1 and β = −0.012 can reach, is much lower compared to those for the additional β values. It is obvious that all results follow closely GR and deviations of the ISCO for any physical radii between the considered set of parameters and GR are marginal.
More important than the actual radius of the ISCO  from observational standpoint is the orbital frequency obtained at its radius. Fig. 10 presents the results for the three cases of β studied following the computation of eq. (24) in the Jordan frame of the theory at f = 0 Hz. Once again, deviations from GR are minor. For masses exceeding 1.3 M ⊙ , the dependence for the topological neutron stars is essentially the same as that for the GR case, while for lower masses, particularly around 1 M ⊙ , the frequencies can be up to 3 % higher than those predicted in GR. These results are not surprising since the deviation is in the region where ISCO is at the star's surface and Fig. 2 clearly indicates that the topological branch of neutron stars are more compact than their GR counterparts particularly at lower masses. This small difference, however, does not appear to be experimentally significant.  The radial epicyclic frequency's maximum value outside of the star at f = 0 Hz is given in the following Fig.  11. The highest deviation is once again around 1 M ⊙ reaching as much as 6% and once again becoming negligible at higher masses. Since neutron stars with masses lower than that are not observed and currently considered as non-existing, the radial epicyclic frequency does not offer any significant observational difference between GR and topological neutron stars similarly to the orbital one.
The following three figures display the same quantities in the case of slow rotation. The color schemes are the same as on figures 9, 10 and 11 with the solid lines corresponding to the f = 0 Hz case, while the dashed and dash-dotted lines give the same dependencies at f = 80 Hz and f = 160 Hz. The results are obtained in the slow rotation regime but the majority of observed neutron stars fall within this approximation. There is also reason to believe that results at higher order of Ω are not qualitatively different [11]. The ISCO for prograde orbit in all cases is shown on Fig.12, followed by the orbital frequency and the radial epicyclic frequency on Figs. 13 and 14.
The decrease in ISCO radius due to the star's rotation leads to an overlap between the plots in an even larger region making the GR and the topological cases harder to discern without knowing the rotational frequency precisely. For each of the rotational cases by itself, no significant distinction is possible once again.
As far as the orbital and maximum radial epicyclic frequencies are concerned -no major difference can  Colors scheme follows Fig.10.
be seen between GR and the n = 1 topological neutron star branch with coupling function (12) for masses higher than ∼ = 1.2 − 1.3M ⊙ . Once again, in the region 1 − 1.25M ⊙ minor quantitative differences on the order of 6 % and less can be observed for both but these do not appear to be significant experimentally. It appears that the monotonic coupling function (12)   does not predict significant qualitative or quantitative deviations from GR as far as the X-ray observations are concerned. As discussed, other coupling functions, such as the non-monotonic one given by eq. (11), predict both qualitatively and quantitatively significantly different results from GR. The new effects that appear poses considerable interest and so it will be investigated further in a follow-up work.

V. CONCLUSIONS
In the present work we have further explored the properties of topological neutron stars (TNS) that is a new and very interesting class of compact objects in the Tensor Multi-Scalar Theories (TMST) of gravity. We considered the case of slow rotation and in addition we expanded the study of these objects' dependence on theory parameters, topological charge and coupling functions. The main focus of the paper was on the investigation of their moment of inertia and the quantities related to accretion under slow rotation.
We have confirmed the EOS independence between the suitably normalized moment of inertia for the TNS and their compactness, comparing the obtained results with those for neutron stars in General Relativity. It turns out that for one of the coupling functions, showing a monotonic behavior with respect to the scalar field, these universal relations are almost indistinguishable from GR, making them not only EOS independent, bus also theory of gravity independent up to a large extend. For more complicated non-monotonic functions, though, the devia-tions can be signification allowing not only to distinguish between GR and TMST, but also between solutions with different values of the topological charge. This demonstrates that the future observations can help us test the very interesting hypothesis that neutron stars can posses a new property that is the topological charge.
We have also obtained accretion-relevant properties such as the ISCO radius, orbital and epicyclic frequencies for TNS at different rotation rates, comparing them once again to the GR results. The obtained TNS quantities differ significantly from those in GR only for some coupling functions and theory parameter values which leaves several observational traits to be sought after and explored. In the case of monotonic coupling function there are no qualitative difference between GR and TMST and in addition, the qualitative deviations from GR are very small, practically unmeasurable. Some interesting qualitative differences appear in particular for non-monotonic coupling function for the ISCO-related quantities, where the quantitative deviations can be large as well, that will be further explored in a dedicated work as they require more profound investigations. This provides a firm reason to pursue more accurate observations and systematic studies in order to constraint alternative theory parameters in the strong regime of gravity.
Further study of the higher order corrections to rotational neutron stars will additionally allow us to get their quadrupole moment and explore potential I-love-Q relations [44,45] also opening new prospects of experimental confirmation through the gravitational waves observation data which is promising to become more abundant in the near-future [46,47]