Building post-Newtonian neutron stars

Owed to their compactness, neutron stars involve strong gravity and extreme density physics. Nevertheless, at present, there are a variety of problems where progress (at least conceptually) can be made in the context of weak gravity. Motivated by this we examine how accurately one can model neutron stars using the post-Newtonian approximation to general relativity. In general, we find there is a significant degree of freedom in how the post-Newtonian equations of stellar structure can be formulated. We discuss this flexibility in the formulation and provide examples to demonstrate the impact on stellar models. We also consider the (closely related) problem of building neutron stars using isotropic coordinates. In this context, we provide a new strategy for solving the equations (based on a scaling argument) which significantly simplifies the problem.


I. INTRODUCTION
Neutron stars involve extremes of physics-in terms of density, magnetic fields, temperature and so on-that make them interesting from many different perspectives.Their special place in the pantheon of astrophysics, nuclear and gravitational physics is justified by a rich observational phenomenology and a range of outstanding issues that remain to be explained/understood [1].The recent focus on neutron stars as sources of gravitational waves, both in binary systems (like the celebrated case of GW170817 [2]) and as emitters of continuous waves [3], bring a number of modelling questions into play.The specific technical challenges may differ but the overarching question remains the same: how do we appropriately represent the physics required to describe a given scenario?By necessity, this has become a gradual grind towards (hopefully) eventual enlightenment.
With the many difficult questions in mind, the scope of the present discussion may seem quite modest.We want to understand how accurately-or not-we can describe neutron stars within the post-Newtonian (pN) approximation (formally, a weak-field, slow-motion expansion of general relativity) [4].At first sight, this may seems like a somewhat misguided exercise.After all, due to their large mass M and small radius R, neutron stars are firmly beyond the pN regime.With a compactness of order (where G is Newton's gravitational constant and c is the speed of light) it is not at all clear that second order pN terms (of order 1/c 4 in a formal expansion) can be neglected compared to first order ones (1/c 2 ).This is, of course, well understood and a good reason to use a fully relativistic approach to studies of neutron-star dynamics.Having said that, one may still be motivated to consider the pN problem.First and foremost, it is important to establish how the star responds to the presence of a binary partner and what this implies for the standard pN expansion for the gravitational-wave signal from an inspiralling binary.This problem is typically expressed as an expansion for small orbital frequencies [4], while the star's response to the tidal interaction is more immediately connected to C and/or the characteristic frequencies of neutron-star dynamics.This "mismatch" is well illustrated by the role of the (static) tidal deformability, which formally enters at a high (orbital) pN order but which nevertheless is expected to be detectable with advanced instruments [5].A closely related issue that springs to mind is the mode-sum approach to the dynamical tides in neutron-star binaries, see [6] for a recent discussion.In the Newtonian case, the solution to this problem relies on the orthogonality and completeness of the star's oscillation modes [7,8].Both properties are complicated in general relativity-in fact, the modes are not going to be complete due to the existence of the late-time power-law tail in the response to a perturbing agent [9,10]-but will remain so in the pN description (at least up to, and including, 2pN terms, which is further than we intend to proceed here).The mathematics will be messier than in the Newtonian case, as exemplified by [11], but so be it.The main question becomes one of precision-are pN model "accurate" enough to be useful for quantitative models of (say) neutron-star tides?This is not a new question (see, e.g., [12], although their focus on alternative theories of gravity is rather different from what we consider here) but we hope to demonstrate that it is worth revisiting.Immediate motivation is provided by the alternative pN descriptions in the literature-from Chandrasekhar's classic work on pN hydrodynamics (which strictly truncates the equations at 1pN order) [11] to the modern description in the monograph from Poisson and Will (which, as we will discuss, retains some 2pN terms in the first order equations) [4] and the formulation of Blanchet, Damour and Schaefer which aims to stay close to the variables used in numerical relativity (and which also mixes in 2pN terms) [13].

arXiv:2209.05871v2 [gr-qc] 1 Mar 2023
If we assume that any 2pN terms may be included without penalty (which is true in the regime where the pN approach truly applies) then the different formulations are all equivalent.However, for actual neutron-star models the inclusion of 2pN terms makes a numerical difference.Hence, it is worth asking if one of the pN prescriptions fares (numerically) better than the others when we employ a realistic matter description, or are they all equally "bad"?This is the question we explore in this paper.It may not seem particularly interesting at first, but a number of subtle issues come into play, making the detailed analysis worthwhile.
The paper is organised as follows.We begin in section II with a summary of the traditional approaches to solving the stellar structure equations in Newtonian and relativistic gravity.This motivates the truncated form of the pN equations.In section III, we formulate the fully relativistic problem (and its corresponding truncated pN form) in isotropic coordinates.This may be a brief aside, but it is important as the pN formulation is commonly based on the use of isotropic coordinates.In section IV, following the standard pN formulation of fluid dynamics, we discuss the (perhaps surprising) flexibility that is afforded in pN theory.We demonstrate how this flexbility impacts on constructed neutron-star models.Finally, we conclude in section V.

II. HOW TO BUILD A STAR
The issue we want to explore involves solving the equations of hydrostatic equilibrium for a self-gravitating body.This is (obviously) a textbook problem so it might be tempting to bypass the introductory aspects and cut straight to the chase.However, in order to understand the issues that arise later it is helpful to start with a few comments on what may be very familiar results.

A. The Newtonian case
Not wanting to miss an opportunity to set the bar as low as possible, let us start with the Newtonian gravity problem.Then we have the mass density ρ, which leads to the mass M N through For reasons that will become clear later, we are using ξ to represent the radial coordinate here.We also need the pressure p, which in equilibrium follows from and the gravitational potential U N , obtained from These three equations form a system that, when complemented by a matter equation of state p = p(ρ), can be solved to obtain a stellar configuration.The surface of the star-at radius ξ = R s -is simply obtained from p(R s ) = 0.The salient point (for future reference) is that we can solve for the mass and radius of a star without involving the gravitational potential: the equation for U N decouples from the other two.This is convenient because the potential involves the a priori unknown value at the centre of the star, U 0 = U N (ξ = 0).After integration, the value of U 0 is fixed by matching the solution to the external potential So far, so good.Unfortunately, Newtonian gravity does not provide an accurate representation of neutron stars.We need to use a realistic equation of state (obtained from nuclear physics) and it is well known that the associated internal energy has a significant impact on the star's mass and radius.Given this, let us consider the (similarly textbook) problem in general relativity.

B. The TOV equations
The equations that describe relativistic hydrostatic equilibrium are known as the Tolman-Oppenheimer-Volkoff (TOV) equations.They are typically expressed in terms of Schwarzschild coordinates associated with the (general spherically symmetric) line element where we are (again) using ξ to represent the (now circumferential) radial coordinate.The problem involves two metric potentials, ν and λ, both functions only of ξ.The standard approach is to introduce a mass M g such that where ε is the total energy density.It then follows that which matches the Schwarzschild exterior if we identify M g (R s ) as the gravitational mass, in turn defined by the asymptotic behaviour of the spacetime metric.
The pressure now follows from where the total energy can be written with Π the (specific) internal energy.Finally, the remaining metric potential is determined by integrating (2.10) The pattern is familiar from the Newtonian case.The potential ν decouples from the equations for the mass and the pressure.We can solve the first two equations for the mass and radius of a star.This is useful because the third equation involves the initially unknown value of the potential at the centre, ν 0 = ν(ξ = 0).This value has to (again) be determined by matching to the external gravitational field at the star's surface.
The solution to the TOV equations is quite straightforward and it may seem self-evident that it should be so.However, there are subtleties here.For example, it is worth noting that e −λ → 1 as ξ → 0 (2.11) ensures that the solution is spatially flat at the star's centre.Also, it is not immediately obvious that that the mass function M g represents the gravitational mass: Intuitively, the total energy contained inside the given volume.In order to argue that this is, indeed, the case it is natural to take a detour and consider the Newtonian virial theorem.

C. The (Newtonian) virial theorem
As will become evident later, one may want to consider different mass functions.One way to do this, certainly within the pN framework, is to make use of the standard virial theorem.Starting from the gravitational potential energy Ω, we have Alternatively, we may use the Poisson equation to get (noting the limits of integration) In summary, we have The key point is that, while the integrands may not be identical for a given value of ξ, when integrated over the star these relations should hold.We can use the virial argument to shed light on the mass used in the TOV equations (following an argument from Box 23.1 in [14]).In Schwarzschild coordinates we have (2.17) Making use of the virial relations-at the integrand level and discarding higher order pN terms-we have From this we identify the total energy as the sum of the mass, the internal energy and the gravitational potential energy.As these are the expected contributions to the total energy for a static star, it makes sense to conclude that the mass M g represents the gravitational mass.The argument is, of course, only valid at 1pN order, but this is good enough for our purposes.

D. Truncating the TOV equations
As we are interested in making contact with the pN expansion, we can truncate the TOV equations at some given order of 1/c 2 [12,15].As a first step in this direction we may expand the right-hand side of (2.8) as This truncated TOV problem was explored by Shinkai [15], who introduced ε = ρ t c 2 to get and then studied the convergence of the expansion as higher order 1/c 2 terms were added in.We are not going to follow precisely this strategy because, in the pN approach it seems more natural to identify the rest-mass density ρ = m B n, where n is the baryon number density1 and m B the baryon mass, and the internal energy density ρΠ as in (2.9).This then leads to (see the post-TOV discussion by Glampedakis et al. [12]) Notably, this equation is not "consistent" in the pN sense as the mass M g itself contains order 1/c 2 contributions, thereby introducing 2pN terms.As an illustration of the accuracy of the truncated model, consider the results in figure 1, obtained for the BSk24 equation of state [16].Although we focus on a single equation of state for the nuclear matter-the same in all examples we provide-we have considered a wide range of matter models to confirm that the qualitative results remain the same.Evidently, the truncated model is good for stars below about 0.5M or so.
As we increase the density the mass-radius curve diverges from the TOV result also demonstrating the absence of a maximum mass configuration.This is as expected for truncated/pN models.For comparison, we also show the mass-radius curve obtained from the Newtonian equations.As an approximation, the Newtonian result clearly breaks down for even lower masses.This is as it should be.There are no real surprises here.
We have already touched upon the notion of different mass functions, hinting that we may profitably work with new-perhaps less intuitive-variables.As an example, let us introduce motivated by the right-hand side of (2.21).This then leads to a different truncated model where Incidentally, this form of the equations was first derived by Wagoner and Malone [17] within the parameterised pN formalism (we will return to this derivation later).Clearly, these equations also include higher order pN terms, although they are not the same as in (2.21): the two models differ at order 1/c 4 .At the surface of the star, we now extract the gravitational mass from (2.25) Does this new mass function make a difference?Evidently so.The example provided in figure 1 shows that (2.23) provides much more accurate stellar models than the alternative truncation (2.21).This is interesting, but at the same time somewhat disturbing.If we can use the freedom to work with different mass functions, basically differing at the higher-order level, to our advantage then how do we know that we are making the optimal choice?

III. ISOTROPIC COORDINATES
Turning to a different aspect of the standard pN formulations, but for the moment staying in the comfortable environment of general relativity, let us see what happens if we try to solve for hydrostatic equilibrium using isotropic coordinates (a problem also considered when constructing initial data for numerical relativity simulations, see e.g.[18]).One reason for doing this is that the results will later provide a direct comparison with the pN results (rather than relying on a "translation" from Schwarzschild coordinates).We will also find the exercise instructive in its own right.

A. Decoupling the metric potentials
In isotropic coordinates, the spatial part of the metric is taken to be conformally flat.That is, we have Again, the gravitational mass is inferred from matching to the vacuum result.In the vacuum exterior, we now have where M g is the gravitational mass (as before) and we see that the transformation to Schwarzschild coordinates (with radial coordinate ξ, as before) is . (3.3) In the weak field limit we get which nicely connects with the static pN metric explored later.
Let us now work out the equations for hydrostatic equilibrium.After a bit of algebra (not provided here as the steps are standard) we arrive at an equation for the first metric potential where we have defined the mass M through The second potential follows from where we have defined a different mass function, M , through with Finally, the pressure follows from where we have defined In these equations we have tried to use variables that make the results appear as "neat" as possible.This evidently includes two different masses, M and M .It is also notable that, in this version of the equations neither of the metric potentials decouple from the pressure.This makes the isotropic problem-at least at first glance-rather different from the standard TOV equations.If we want to solve the equations as given above, then we need to provide values for both potentials at the centre of the star, µ(r = 0) = µ 0 and ν 0 .In effect, the solution involves a two-dimensional root-search, matching the integrated variables to the external metric (see [18] for relevant comments).
At the surface of the star, r = R, we need both µ and its derivative to match continuously to the exterior (continuity of the first and second fundamental forms).We want to ensure that and Starting from the latter, for a chosen value of µ 0 , we arrive at the interior solution The matching at the surface then leads to and it follows from (3.12) that This provides the condition needed to (iteratively) determine the value for µ 0 .However, there is a snag.In order to solve the equations we need the energy ε(r), which follows from the solution for the pressure (and the assumed equation of state).Moreover, the location of the surface follows from p(r = R) = 0.This then brings in the second, a priori unknown, central value ν 0 .It would clearly be better if we could decouple at least one of the potentials from the equation for the pressure (as in the TOV case).One way to go about doing this is to define yet another mass function.With this in mind, we introduce and Working with M we decouple the second potential, ν, from the pressure equation.This then allows us to solve the problem in two stages, where the first step only involves determining the appropriate value for µ 0 .Once this is done, the second potential follows from (3.7)-(3.9).Finally, the value of ν 0 needs to be such that Equivalently, this can be written which leads to Making use of the interior equations, this condition corresponds to At the end of the day, a consistent solution requires The mass functions may differ inside the star, but they all lead to the same result at the surface-evidently in the spirit of the Newtonian virial relations.

B. Scaling argument
The arguments we have provided simplifies the strategy for solving for hydrostatic equilibrium in isotropic coordinates.Essentially, we need two root searches, first for µ 0 and then for ν 0 .This is easier to implement than the original two-dimensional problem.However, the logic is still more involved than we are used to from the TOV problem.It is natural to ask if we can do better.It turns out that we can.
First, introduce and rewrite (3.5) as The other equations we need now take the form Second, ask if the equations are invariant under some suitable rescaling.The argument proceeds as follows.Scale the variables in such a way that where β, α and η are all constant.From the first two equations, we then see that we need to have leading to Next, in order for the pressure equation to be invariant, we must have which is consistent with the final equation, as well.
In essence, once we have a solution to the equations we can find other solutions by scaling.We can use this "trick" to ensure the matching at the star's surface.Starting from an arbitrary value for µ 0 , integrate the equations up to the point where p(R * ) = 0 and M (R * ) = M * .We also obtain the potential ψ * .
At the surface of the star the true values of the variables should be such that we see that we must have which also leaves the equation for N invariant.The surface matching, adding a suitable central value for ν 0 to the previously calculated model, now requires and the scaled solution leads to which fixes the value of ζ and completes the stellar model.Alternatively, we can use the fact that the final result should be consistent with M (R) = M (R), which leads to C. Truncated isotropic models The scaling argument leads to an elegant prescription for building stars using isotropic coordinates, notably not involving a root search at all.It also helps us understand what happens if we expand the equations to some power of 1/c 2 , as in the case of the truncated TOV equations from section II D. This turns out to be an instructive exercisebringing out issues relevant also for the pN discussion to follow.
Focussing on the mass-radius relation, we first note that we can leave out the step that determines the second metric potential N = e ν/2 .That is, we focus on the system of equations (3.5)-(3.6)and (3.17)- (3.18).We do, however, need to commit to an expansion for the potential ψ = e µ/4 .Intuitively, as it connects with the pN logic, it is natural to assume that where Ū plays the role of the gravitational potential.
In order to check to what extent additional assumptions make a difference, we will compare two slightly different models.First, suppose we do not separate the energy into mass and internal energy parts, but leave ε and the combination with the pressure as they are.Then we need to keep in mind that ε/c 2 = O(1) to get and Once these equations are integrated to the surface p = p(R) = 0, we need to match to the external potential at the surface.We then need where we have to keep in mind that we do not yet know M g , but It follows that the central value of the potential ( Ū0 ) is fixed by the condition At this point-and for future reference-it is important to note that the expansion breaks the scaling argument from the previous section.Basically, we now have the freedom to change the central value of the potential Ū by adding a constant.The upshot of this is that the correct central value Ū0 must be iteratively determined by a root search.
In the second model, we expand the energy in the usual way.Using (2.9) we have ) and The matching condition at the surface remains as before.
There should be nothing particularly controversial here.The steps we have taken seem natural.However, when we solve either set of truncated equations-(3.45)-(3.48)or (3.52)-(3.55)-theresult is surprising.Basically, we find that there is a maximum density beyond which no stellar model can be found.Typical results (again for the BSk24 equation of state) are provided in figures 2 and 3.The behaviour is very different from the solution to the full isotropic version of the TOV equations (which leads to a mass-radius curve identical to the original TOV solution from figure 1 once we translate the coordinate value of the radius).What is happening here is that the expansion in 1/c 2 has made the eigenvalue problem for Ū0 nonlinear.At low densities there are two roots.One root brings us back to the Newtonian limit and the other leads to very massive stars with a low value of the central density.As the density increases, these two sequences approach each other until-at a specific maximum value of the density-they merge.Beyond this density there are no solutions.In figures 2 and 3, we have ended our calculation when the two roots are closer than a set numerical resolution.This leaves a gap between the two sequences.While we could make this gap smaller by pushing the calculation to higher precision, we have left it like this to illustrate the main take-home message.
On the one hand, one may argue that the result is what it is.We should not be surprised that the truncated solution differs from the full "nonlinear" one when higher order terms in the 1/c 2 expansion play a prominent role.The result is just a warning that we should not use the truncated models in the regime where they do not apply.The conclusions are nevertheless disturbing, especially as we are about to turn to the pN problem, which is similar in spirit.There are many situations where one might want to use an approximate neutron-star model, even if it may not be particularly precise.What we are learning is that such approximations may have less attractive features.Consider the simple fact that, for a given equation of state we can build neutron-star models up to the maximum mass (and beyond) using the TOV equations (in Schwarzschild or isotropic coordinates, whichever may be preferred).Similarly, we can find solutions to the Newtonian equations for any given central density.These solutions may not provide a very accurate representation of neutron stars (there is no maximum mass, etcetera), but they exist.The "back-bending" of the curves in figures 2 and 3 beyond some density seems an unintuitive and not particularly welcome feature of the truncated problem.Having said that, if all we want is a reasonably good approximation then it is worth noting that we can use the truncated models to build canonical neutron stars with a mass of 1.4M with errors in the radius less than 1 km.This level of error is notably similar to the current observational constraints, obtained from multimessenger data (including the NICER radius results for the low-mass pulsar J0030+0451 and the high-mass system J0740+6620) [19].This suggest that, while the pN models may be viewed as "good enough" at the moment, we have to do better in the future. .For reference we also so the results from the isotropic TOV equations (blue dot-dash).Data correspond to the BSk24 equation of state.

IV. FORMULATING THE PN PROBLEM
At this point, we could continue to explore the weak-field limit of the isotropic TOV equations with the aim to gain further insight into the 1/c 2 expansion.Rather than doing this, we will move straight to the pN version of the problem.The main difference is that, instead of truncating the relativistic equations we will build on the pN formulation of hydrodynamics.In doing this we heavily rely on the exhaustive discussion in the recent monograph by Poisson and Will [4].In this description the baryon number density plays central role, so it makes sense to begin by explaining how this comes about.

A. Baryon number conservation
Suppose we start from the perfect fluid stress-energy tensor where u a is the fluid four velocity, along with the baryon number conservation law Now introduce a different observer such that where γ = u 0 /c follows from the normalisation condition .For reference we also so the results from the isotropic TOV equations (blue dot-dash).Data correspond to the BSk24 equation of state.
We then have (keeping in mind that x 0 = ct) or where ρ = m B n is the baryon mass density (as before) and we have defined So far, the description is fully relativistic so equation (4.6) should hold at all pN orders [4].This makes ρ * an attractive variable to work with.

B. The pN formulation
Now consider a formal expansion in inverse powers of c.We then need the pN form for the metric (there are different choices, depending on the chosen gauge [20]).Here we take equations (8.2) from [4], which are expressed in the harmonic gauge, as our starting point.We then have ) and where we have introduced the potentials2 and where v i is the fluid (three-)velocity.We also have where and the superpotential X is determined by Given these definitions we find that If we also note that we see that and, for future reference, it is worth noting that

C. Static fluid configurations
For a static fluid configuration, we ignore all time derivatives and also set v i = 0. Noting that-this is evident from (4.10)-the pN metric is expressed in isotropic coordinates, letting the corresponding radial coordinate be r (as before) and following the steps in [4] we arrive at the equation for hydrostatic balance at 1pN order: That is, we have where We also note that U i = 0 and Ψ = ψ so we do not have to worry about the superpotential in the following.Conveniently, we may write the equations we need to solve as a coupled set of first order differential equations.First we have with M B representing the baryon mass (recall ρ = m B n).This is evident since where dV is the (three-)volume element for the pN metric.Next it follows that In the vacuum outside the star, we must have-regardless of which mass we use to source the gravitational potential (see later):

.26)
The constant A depends on the source of the potential, but the surface matching condition does not.We only need to make sure that That is, we need to hold for all pN models.We also have the two equations and With these definitions (4.21) takes the form These are the pN equations as given in exercise 8.5 of [4].There are no great revelations here, but three comments are in order.First, the equations we have written down clearly include higher order pN terms.This is not a "problem" as long as we are truly in the pN regime, where higher order terms in 1/c 2 are small, but we have already seen from the example of the truncated TOV equations that such higher order terms may significantly impact on the neutronstar models.Inevitably, the inclusion of higher order terms renders the modelling more ad hoc than one might be comfortable with.Second, as we need to solve for the potential U we need to pay attention to the fact that the central value of the potential is not known ab initio (just as in the isotropic case).It needs to be solved for by matching at the surface.The third point relates to the equation of state.
In equilibrium, neutron-star matter can be described in terms of a single parameter, usually taken to be the baryon number density.Thermodynamics then provides the Gibbs relation or In terms of the specific internal energy per unit mass Π, it follows that This leads us to issue a word of caution regarding the variable ρ * .It may seem tempting, given the form of the pN equations we have written down to think of the matter equation of state as a relation p = p(ρ * ).Mathematically, this might make sense but physically it is problematic.The matter description is developed in a local inertial frame and therefore the equation of state cannot be cognizant of the spacetime geometry.It should not depend on the potential U .
Let us now solve the pN equations in this form.When we do this, we find-perhaps not surprisingly-that the solution has the same behaviour as in the case of the truncated isotropic equations.Beyond some maximum central density there are no solutions.This is illustrated by the results in figure 4 which shows the surface matching condition as function of the central value of the gravitational potential.At lower densities the problem has two roots (and hence there will be two stars with the same central density), but at higher densities there are no solutions.
In order to fix the problem, we can try to remove some of the higher order pN terms.For example, it is easy to see that we may rewrite (4.31) (equivalent at 1pN order) as Results from this equation are slightly better, see figure 5, but the problem of the two roots remains.We also show how the mass depends on the baryon number density in figure 6.Noting that the radius of a 1.4M star is now closer to the correct result, it is evident that the inclusion (or not) of higher order pN terms can make a substantial difference.
The two examples we have provided inevitably prompts us to ask if there is an optimal pN model.This is a tricky issue, given the somewhat arbitrary choices involved.Still, it is relevant to consider another couple of options.
An alternative would be to take the lead from Shibata [21], go back to (4.20) and rewrite it as (again, consistently to order 1/c 2 ) Instead of expanding this equation, we then divide3 by the prefactor on the left-hand side.This then leads to Finally, expanding this result in pN orders brings us back to (4.35).Results obtained from (4.37) are also shown in figure 5.This model is notably better.In fact, as can be seen from the results for polytropes in [21] the model also leads to a maximum-mass configuration.This is an attractive feature, even though the predicted maximum mass is much too high.In principle, the results obtained from (4.37) are attractive but one may have to view this with some caution as it is not clear that a strategy that involves dividing by terms involving pN corrections is "helpful" once we turn to dynamical problems.FIG. 4. Demonstration that there is a maximum central density (n0) beyond which the internal gravitational potential cannot be matched to the exterior in the pN scheme from [4].Below a certain density the problem for the central value U0 has two roots (the solid blue curve crosses zero at two points), but above a critical density there are no roots (orange dotted curve).This accords with the results for the isotropic TOV equations from figure 3. Data correspond to the BSk24 equation of state.

D. Different mass functions
A useful lesson-apparent already from the discussion of the truncated TOV equations in section II D, see also [12]-is that we may use different mass functions to our advantage.Our goal is to explore, given the freedom afforded by the pN approximation, how different definitions at 1pN enable us to construct more, or (indeed) less, accurate neutron-star models.With this in mind, we note that the pN variables we have used may not be the most "intuitive".For example, comparing to the relativistic result one would not expect the gravitational potential to be sourced by the baryon mass.Given this, let us introduce Using this mass to source the potential, we have The new mass M is identical to the ADM mass (by definition equal to the gravitational mass) at 1pN order.To demonstrate this we may use, for example, equation (8.136) in [4] (see also [21]), from which we have We may also derive the result from full relativity.The starting point would be the Tolman mass [22,23], defined as where we have used the virial relations from section II C.This argument shows that M obtained from (4.39) is, indeed, the 1pN version of the gravitational mass.Using the virial result, it is also easy see that (4.39) leads to the same expression for M .The different results are consistent.The main insight at this stage is that the two masses M B and M have clear physical meaning.We now get  we arrive at concise relation Results obtained from (4.43) are provided in figure 7.This model does not seem to be an improvement on the versions considered in figure 5.
Finally, working with M we may also use Mass-radius results obtained from this model are provided in figure 7.This model performs rather well compared to the alternatives we have considered.As it does not suffer from the maximum density issue, this is an attractive (and comparatively simple) formulation of the pN problem.Finally, we may ask if it is possible to take the notion of different masses one steps further.We have shown that one of the problematic features of the pN schemes relates to the a priori unknown central value of the gravitational potential.Given this, it is natural to ask if it is possible to decouple the potential from the equation for the pressure (as in the TOV equations).It turns out that we can, indeed, do this and the argument ultimately leads to us back to (2.23).As we have already considered this model, the step-by-step derivation is mainly relevant for completeness.The interested reader will find the derivation in the Appendix.

V. SUMMARY
Neutron stars, with their strong gravitational fields, lie firmly in the relativistic regime.However, while a precise description of their structure and dynamics must be relativistic, there are a number of contexts where the modelling is technically challenging and where progress can be made by assuming weak gravity.With this in mind, we have explored the extent to which neutron stars can be described within the pN approximation.Indeed, working to 1pN order, there is considerable freedom in how one formulates the equations (see, e.g., Refs.[4,11,13]).A particularly important aspect relates to the freedom to introduce different mass functions in the stellar interior.We have shown that the choice of mass function may have a significant impact on the accuracy of stellar models (in terms of mass and radius) and can, in some cases, lead to results that are fairly close to the fully relativistic result.In particular, our results show that neutron-star models based on (4.46), which follows from a fairy simple prescription, are significantly more accurate than models built from (4.31), which may be seen as the "standard" approach [4].Here we demonstrated the result for a single matter equation of state, but we have confirmed it for a wide range of matter models.
As a step towards the pN discussion, we explored the use of isotropic coordinates.In contrast to the familiar TOV equations, the relativistic equations in such a coordinate system couple the two metric potentials.The problem also naturally involves two distinct mass functions (both of which match the gravitational mass at the surface).We have devised a new strategy for solving this coupled boundary-value problem, involving a scaling prescription that neatly circumvents the need for a root search.
Motivated by the usual approach in the pN approximation, we explored a truncated version of the isotropic TOV equations.Perhaps surprisingly, this uncovered a pathology in the isotropic formulation of the pN equations: there (typically) exists a maximum central density beyond which one cannot find a numerical model that satisfies the boundary conditions.We showed that this problem occurs because the pN expansion at these densities leads to a non-linear eigenvalue problem for the central value of the gravitational potential.Below the maximum density, there are two roots for a given central density, corresponding to a more and less massive star (only one of which is physical in the low-density limit).
Moving forward, one of the motivations for this work was to study the dynamical tides of neutron stars beyond Newtonian gravity (with a realistic matter model).While it is clear that the pN strategy for building neutronstar models is far from robust (given the sizeable difference in results obtained from models that are-at least formally-consistent to the same pN order), the approach does not suffer from the well-known issues that hamper the development of a fully relativistic description of dynamical tides.In particular, the set of stellar oscillation modes will remain complete at 1pN order, thus providing a potential avenue to develop more accurate descriptions of dynamical neutron-star tides using the mode-sum approach [6].Work in this direction will be relevant for future, sensitive gravitational-wave observations of compact binaries.
The purpose of this Appendix is to show that, if we set out to define a new mass function for which the gravitational potential is decoupled from the equation for the pressure, then we are naturally led back to equation (2.23).
The argument we need can be found in work by Wagoner and Malone [17] and Ciufolini and Ruffini [24,25] (see also the more recent discussion by Glampedakis et al. [12]) within the parameterised pN approach [20].This strategy is commonly used to quantify the impact of deviations from general relativity and design tests of the theory.Our interest here is quite different.We want to explore different approximations to general relativity, rather than alternative theories.However, the scope does not matter.The steps we need are the same.
We basically want to ask if we can find a coordinate transformation that decouples the gravitational potential U N from the hydrostatic equations, taking as our starting point the (consistent) 1pN equation where Motivated by the transformation from isotropic to Schwarzschild coordinates from section III (noting that the TOV equations in Schwarzschild coordinates decouple from the metric potential ν), we consider The final step then is to introduce yet another mass function to remove this dependence on the energy from the equation for the pressure.The argument we need is discussed in detail by Glampedakis et al. [12].We also need the gravitational potential energy, now given by Ω = − (5.20) If we now choose A = C = D = 1, B = 0, we retain the expressions from Wagoner and Malone [17], i.e. the pressure is determined from (2.23) which we already derived from the truncated TOV equation in section II D. This makes the argument somewhat "circular" but it is rewarding to see that the expressions can be derived from the 1pN relation (5.1).

N = e ν/ 2 (c 2
correct central value for ψ (and hence µ).The other variables for the stellar model follow from the appropriate scaling.A similar argument applies to the second potential.Introduce (ε + 3p) N ψ 6 .(3.38) Assuming a scaling such that N → ζN (3.39)

2 FIG. 2 .
FIG.2.Mass-radius curves for the two first order truncations of the equations for hydrostatic equilibrium in isotropic coordinates, (3.45)-(3.48)(orange dotted) and (3.52)-(3.55)(green solid).For reference we also so the results from the isotropic TOV equations (blue dot-dash).Data correspond to the BSk24 equation of state.

2 FIG. 6 .
FIG.6.Mass against central number density (n0) curves for the pN schemes (4.31) (orange dotted) and (4.35) (green solid).Also shown are results obtained from (4.37) (red dashed).For comparions we also show the results from the isotropic TOV equations (blue dot-dash).Data correspond to the BSk24 equation of state.
ρ −U N − ξ dU N d ξ + Π + 3 p ρ d ξ + Π + p ρ + U N − ξ dU N dξ ξ 0 4π ξ2 ρd ξ .(5.8)It is now natural to define a mass function m byd m dξ = 4πξ 2 ρ ,(5.9)such that M N = m + O(c −2 ).At this point, we need to do some integration by parts (making use of the relations from the virial theorem).First, we have what we wanted: an equation for hydrostatic equilibrium that does not explicitly involve the potential U N .The final expression is, however, not quite satisfactory as it requires the integral for the internal energy E = ξ 0 4π ξ2 ρΠd ξ .(5.16)