A Spectral Approach to the Relativistic Inverse Stellar Structure Problem

# A Spectral Approach to the Relativistic Inverse Stellar Structure Problem

Lee Lindblom and Nathaniel M. Indik Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125
July 12, 2019
###### Abstract

A new method for solving the relativistic inverse stellar structure problem is presented. This method determines a spectral representation of the unknown high density portion of the stellar equation of state from a knowledge of the total masses and radii of the stars. Spectral representations of the equation of state are very efficient, generally requiring only a few spectral parameters to achieve good accuracy. This new method is able, therefore, to determine the high density equation of state quite accurately from only a few accurately measured data points. This method is tested here by determining the equations of state from mock data computed from tabulated “realistic” neutron-star equations of state. The spectral equations of state obtained from these mock data are shown to agree on average with the originals to within a few percent (over the entire high density range of the neutron-star interior) using only two data points. Higher accuracies are achieved when more data are used. The accuracies of the equations of state determined in these examples are shown to be nearly optimal, in the sense that their errors are comparable to the errors of the best-fit spectral representations of these realistic equations of state.

###### pacs:
04.40.Dg, 97.60.Jd, 26.60.Kp, 26.60.Dd

## I Introduction

The standard stellar structure problem consists of determining the structure of a star by solving the coupled gravitational and hydrodynamic equations with an assumed equation of state for the stellar matter. The solutions to the standard problem determine the various observable properties of the stars with a given equation of state like their total masss , their total radii , etc. The inverse stellar structure problem determines what equation of state is required to produce stellar models having a given set of macroscopic obserables. The goal of this paper is to find efficient and robust methods of solving this inverse stellar structure problem. The method developed here is based on the use of spectral expansions to represent the equation of state. The values of the spectral coefficients in these expansions are fixed in this method by matching stellar models based on these equations of state to observed values of the masses and radii of the stars. Once fixed, these coefficients determine the equation of state that represents the (approximate) solution to the inverse stellar structure problem.

For non-rotating stars in general relativity theory, the simplest version of the stellar structure equations were first derived by Oppenheimer and Volkoff Oppenheimer and Volkoff (1939),

 dmdr = 4πr2ϵ, (1) dpdr = −(ϵ+p)m+4πr3pr(r−2m), (2)

where represents the mass contained within a sphere of radius ; is the pressure; and is the total energy density of the fluid. Solving these equations with a given equation of state is the standard relativistic stellar structure problem. Once an equation of state is specified, these equations determine a one parameter family of stellar models, and , where is the value of the pressure at the center of the star . These solutions then determine various macroscopic properties of the stars, including their outer radii where , and their total masses . These macroscopic properties are (at least in principle) observable.

The standard stellar structure problem can be thought of as a map that takes the equation of state [a curve in energy density – pressure space ], into a curve in the space of macroscopic observables, e.g. . The inverse stellar structure problem consists of finding the inverse of this map Lindblom (1992), i.e. determining the equation of state of the stellar matter from a knowledge of some information about the macroscopic structures of the stars (like their masses and radii ). The solution to this problem, like the solutions to many inverse problems, is less straightforward than the solution to the standard problem.

The inverse stellar structure problem is probably more relevant for practical relativistic astrophysics, however, than the standard problem. The highest density part of the equation of state in neutron stars, for example, is not well known. Matter in this state is well beyond the reach of laboratory experiments, so there is no independent way of directly measuring its properties, including its equation of state. Numerous attempts have been made to model this matter theoretically, but even today there is no consensus among theoreticians. Predictions of the energy density for a typical neutron-star central pressure, for example, often differ by an order of magnitude. Therefore, the standard stellar structure problem for neutron stars is not terribly useful. In contrast, the inverse problem may provide an important tool for learning about high density nuclear matter. Numerous high quality measurements of the masses of neutron stars are now available Lattimer and Prakash (2007), and a few (fairly imprecise and model-dependent) radius measurements are starting to become available as well Steiner et al. (2010); Galloway and Lampe (2012); Guver et al. (2012a, b). In principle then, the inverse stellar structure problem should (eventually) allow us to measure the high density equation of state of neutron-star matter directly. This measurement would provide important information about nuclear interactions that can not be obtained at present in any other way.

One naive approach to solving the inverse stellar structure problem for neutron stars would simply be to match their observed properties, e.g. their data, with models of those stars based on different micro-physical models of the dense material in their cores. In this approach the model equation of state whose stellar models best matches the data would be declared the observed neutron-star equation of state. This approach would clearly be ideal if there were wide consensus on exactly what the high density core material is, and if there were a reasonably simple model for this material that were known, up to a few undetermined parameters that could be fixed by these observations. Unfortunately the wide diversity of “realistic” neutron-star equations of state in the literature, suggest that (in the near term at least) this approach is not likely to be effective or conclusive.

A more practical variation of this approach uses some knowledge about the expected properties of neutron-star matter in an intermediate range of densities, and a more empirical description of the equation of state for larger densities Steiner et al. (2010, 2012). This approach is more promising, but the proposed model equations of state of this type have many free parameters that must all be fit by the observational data. Since these data are likely to be sparse for some time, we take a different approach here.

Our goal is to find efficient and robust methods for solving the inverse stellar structure problem that use no prior knowledge of the high density micro-physics at all. A (somewhat impractical) method for solving the inverse stellar structure problem that uses no information about the micro-physics of the high density equation of state was given in the literature about 20 years ago Lindblom (1992). This traditional method can be summarized as follows. The total masses and radii of all of the stars associated with a particular equation of state are assumed to be known. The equation of state is also assumed to be known up to some pressure with corresponding energy density . Let and denote the mass and radius of the star whose central pressure is . Now choose another point, , along the mass-radius curve, with a slightly larger central pressure. The outer layers of this new star are composed of low pressure material, , where the equation of state is known. The stellar structure equations, (1) and (2), can therefore be solved in this outer region starting at the surface of the star, , where and , by integrating inward toward . This integration can be continued until the point where and the known equation of state ends. This integration determines the radius and the mass of a “small” core of high pressure material, where the equation of state is not yet known. If this core is small enough, the stellar structure equations can be solved there as a power series expansion about . The coefficients in this expansion are functions of the central density and pressure of this little core. Since the mass and radius of this core are known, and , this power series can be “inverted” to determine and  Lindblom (1992). This new point provides a small extension of the equation of state beyond . Iterating these steps then determines a sequence of closely-spaced points along the high density portion of the equation of state curve.

This traditional solution to the inverse stellar structure problem is unfortunately very impractical. A large number of points are needed from the mass-radius curve to achieve modest accuracy in the calculation of the corresponding points along the equation of state curve. Since points are very difficult to measure (at least for neutron stars) it is unlikely that there will ever be enough data to use this traditional method to determine the unknown high density part of the neutron-star equation of state.

This paper proposes a rather different approach to the solution of the inverse stellar structure problem, an approach that can be very effective even when only a small number of data points are available. The equation of state in this new approach is expressed as a parametric equation, e.g. , instead of a table of values . The parameters are adjusted to give the best-fit approximation to a particular equation of state model. Parametric representations of this sort, based on spectral expansions, have been shown to be extremely efficient at representing the high density portions of “realistic” neutron-star equations of state Lindblom (2010). Only a few non-vanishing are generally needed to achieve 1% accuracy in most cases.

The basic idea of this new method for solving the inverse stellar structure problem is to choose the equation of state parameters by minimizing the differences between the masses and radii of real neutron stars, and , with those based on the parametric model equation of state, and . Once the are fixed, the parametric equation then provides an approximate solution of the inverse stellar structure problem. Spectral expansions typically converge exponentially as the number of terms in the expansion are increased (for smooth functions). These approximate solutions to the inverse stellar structure problem are therefore expected to converge to the exact equation of state as the number of data points and the number of parameters fixed by this method are increased.

The remainder of this paper presents details on how to implement this new spectral approach to the solution of the inverse stellar structure problem, along with practical tests of its accuracy and efficiency. Section II reviews the particular spectral representation of the equation of state used in the solution presented here. Section III describes how the spectral parameters are fixed by matching to the given data points. Section IV presents a series of numerical tests of the accuracy and efficiency of this new method. Mock data computed from a collection of 34 “realistic” neutron-star equations of state are used as input in these tests. These tests show, for example, that the resulting spectral equation of state agrees with the exact to within a few percent (on average) when only two data points are used. Higher accuracies are (generally) achieved when more data points are used. Section V discusses some of the limitations of the numerical tests presented here, and proposes several ways that the basic method developed here might be extended and improved. Some of the more complicated technical details needed to implement this method are described in two Appendices. Appendix A describes how to evaluate the derivatives of and , with respect to the parameters and . Appendix B describes the interpolation method used here to bridge the gaps between points in the exact “realistic” equation of state tables.

## Ii Spectral Representations of the Equation of State

The version of the relativistic stellar structure equations most useful for our analysis here requires that the equation of state be written in a form where the energy density and pressure are given as functions of the relativistic enthalpy, . The usual form of the equation of state, , must therefore be re-written as a pair of equations and , where the enthalpy is defined as

 h(p)=∫p0dp′ϵ(p′)+p′. (3)

The needed expressions, and , can be constructed by inverting from Eq. (3) to obtain , and then by composing the result with the standard form of the equation of state, , to obtain .

The transformations needed to construct and in this way are difficult to perform efficiently and accurately in numerical computations. Therefore it is best to construct a spectral representation of the equation of state that is based directly on . This can be done by introducing an enthalpy based spectral expansion of the adiabatic index  Lindblom (2010):

 Γ(h) ≡ (4) = exp{∑kγk[log(hh0)]k}, (5)

where is the lower bound on the enthalpy, , in the domain where the spectral expansion is to be used. This is a standard spectral expansion of the function in which the are the spectral basis functions and the are the spectral expansion coefficients (or parameters).

The functions and are obtained from by integrating the system of ordinary differential equations,

 dpdh = ϵ+p, (6) dϵdh = (ϵ+p)2pΓ(h), (7)

that follow from the definitions of and in Eqs. (3) and (4). The general solution to these equations can be reduced to quadratures:

 p(h) = p0exp[∫hh0eh′dh′μ(h′)], (8) ϵ(h) = p(h)eh−μ(h)μ(h), (9)

where is defined as.

 μ(h)=p0eh0ϵ0+p0+∫hh0Γ(h′)−1Γ(h′)eh′dh′. (10)

The constants and are defined by and respectively. While these quadratures can not be done analytically for the spectral expansion of given in Eq. (5), they can be done numerically very efficiently and accurately using Gaussian quadratures.111The numerical accuracy (and hence the efficiency) of the numerical integrations in Eqs. (8) and (10) can be improved significantly by changing integration variables from to before performing standard Gaussian quadratures. Our tests achieved accuracies about for and with 10 Gaussian integration points using the variable, compared to about accuracies for the same tests using the variable.

The method of solving the inverse stellar structure problem proposed in Sec. III is based on this spectral representation of the equation of state: and . Any equation of state can be represented approximately in this way by using a finite number of spectral parameters in the expression, Eq. (5), for . In analogy with other spectral expansions, the accuracy of these approximations are expected to converge exponentially (for smooth equations of state) as the number of spectral coefficients is increased. Numerical tests that fit “realistic” neutron-star equations of state using this method Lindblom (2010) are consistent with this expectation about the convergence of these expansions.

## Iii Fixing the Spectral Parameters

The spectral approach to the inverse stellar structure problem fixes the equation of state parameters by choosing the values whose stellar models best match a collection of points from the exact mass-radius curve. The process of fixing the in this way can be made more efficient numerically by using a somewhat non-standard version of the stellar structure equations. Therefore, we digress briefly here to review this alternate formulation.

The Oppenheimer-Volkoff version of the stellar structure problem, Eqs. (1) and (2), determines and as functions of . That traditional approach has two inconvenient features: First the integration domain, , is only known after the fact when the surface of the star at is found numerically by solving . Second, the equation that defines the surface of the star is somewhat difficult to solve numerically because typically vanishes at . These inconveniences can be avoided by transforming the equations into a form where and are determined as functions of the relativistic enthalpy (see Ref. Lindblom (1992)):

 dmdh = −4πr3ϵ(r−2m)m+4πr3p, (11) drdh = −r(r−2m)m+4πr3p. (12)

Solving the equations in this form begins by specifying “boundary” conditions, , at the center of the star where , and then integrating toward the surface of the star where . The derivative in Eq. (12) is non-zero and bounded at the surface of the star, so this formulation completely eliminates the problems associated with solving to locate the star’s surface. This version of the problem is also easier to implement numerically because it is carried out on the domain , which is fixed before the integration is performed. The total mass and radius of the stellar model are obtained in this formulation simply by evaluating the solutions and at the surface of the star where : and . More details about how to implement this formulation of the stellar structure problem are given in Ref. Lindblom (1992) and in Appendix A of this paper.

This alternate formulation of the stellar structure problem, Eqs. (11) and (12), requires that the equation of state be expressed in terms of the enthalpy, i.e. that and be provided. The spectral representations, and , described in Sec. II are therefore ideal for this purpose. The general solution to this form of the stellar structure problem is a pair of functions of the form and . These solutions are determined uniquely by the parameter , the central enthalpy of the star, and , the spectral parameters that determine the equation of state. The total mass and radius associated with one of these stellar models are determined from these solutions by and .

The new method of solving the inverse stellar structure problem, which we introduce here, fixes the values of the spectral parameters by minimizing the differences between the model mass-radius values and points from an exact mass-radius curve. Thus we fix the values of the by minimizing

 χ2(hjc,γk) = (13) +[log(R(hic,γk)Ri)]2⎫⎬⎭,

with respect to each of the . Each of the stellar models used in this fit has a central enthalpy, , whose value is also needed (along with the ) to construct the mass-radius values . Since the will not be known a priori, these additional parameters must also be determined as part of the fitting process. These parameters are fixed, therefore, by minimizing with respect to variations in each of the .

In summary then, this new method of solving the inverse stellar structure problem determines the equation of state by fixing the spectral parameters in a way that minimizes the differences between the model mass-radius values and values from an exact mass-radius curve. This minimization problem is equivalent to solving the equations

 ∂χ2∂hic=0, (14)

and the (the number of spectral parameters) equations

 ∂χ2∂γk = 0, (15)

for the parameters and . Since the number of independent and data values used in this fitting process is , it follows that the maximum number of spectral parameters that can be fixed in this way is .

A number of numerical methods for solving non-linear least squares problems such as Eqs. (14) and (15) are discussed in the literature. Many of these methods require only that be provided numerically for arbitrary values of the parameters and . Some methods also require in addition that the values of the derivatives and be provided. The numerical tests described in Sec. IV use the Levenberg-Marquardt method for solving these equations, and this method requires that both the values and the derivatives of be provided.

The derivatives of are determined by the derivatives of and with respect to and . These derivatives can be approximated numerically by expressions of the form, . We find it is more efficient (and more accurate) however to evaluate these derivatives by solving an auxiliary system of ordinary differential equations, that are obtained by differentiating Eqs. (11) and (12) with respect to these parameters. This method of evaluating the needed derivatives of and is described in some detail in Appendix A.

## Iv Testing the Spectral Inversion Method

In this section we test the spectral method of solving the inverse stellar structure problem by analyzing sets of mock data points from mass-radius curves based on known “realistic” neutron-star equations of state. We use these mock data to construct best-fit values for the spectral equation of state parameters , using the least-squares method outlined in Sec. III. Then we compare the equation of state constructed in this way with the exact equation of state used to compute the mock data points. We perform these comparisons using different numbers of data points to determine how the accuracy of the approximate equation of state improves as the number of data points is increased. We construct the mock data points using 34 different realistic neutron-star equations of state to explore how well the method works for a fairly wide variety of different equations of state.

The 34 equations of state used to construct the data points used in these tests are the same ones used by Read, Lackey, Owen and Friedman Read et al. (2009) in their study of approximate piecewise polytropic fits to the equation of state. These 34 realistic equations of state are based on a variety of different models for the composition of neutron-star matter, and a variety of different models for the interactions between the particle species present in the model material. Descriptions of these realistic equation of state models, and references to the original publications on each of these equations of state are given in Ref. Read et al. (2009), and are not repeated here. The individual equations of state are referred to here using the abbreviations used in Ref. Read et al. (2009), e.g. PAL6, APR3, BGN1H1, etc. The list of these equations of state are given in the first column of Table III of Ref. Read et al. (2009) as well as the first column of Table 1 here. Spectral fits have already been shown to provide excellent approximations to these 34 equations of state in Ref. Lindblom (2010). Only two or three spectral coefficients are needed to achieve accuracies at the few percent level. The new question being studied here, therefore, is whether the spectral parameters can be determined by fitting data instead of fitting directly to the equation of state itself.

These 34 exact realistic equations of state are provided to us as tables of energy density and pressure points . We compute the enthalpy values corresponding to the points in these tables, interpolate between these tabulated points whenever necessary, and construct complete enthalpy based equations of state from these tabulated data using the methods described in Appendix B. Whenever we refer to one of these exact realistic equations of state, we mean the one constructed by interpolating the tabulated equation of state data in this way.

We construct sets of mock data points by solving the standard stellar structure problem using each of the exact realistic equations of state described above. Figure 1 illustrates these exact mass-radius data for three of the realistic equations of state: PAL6, APR3, and BGN1H1. We select subsets of these points for each of our tests of the spectral inversion method. We limit the points used for our tests to small numbers of models (since we do not anticipate that large numbers of observations are likely to be available any time soon) that fall within the astrophysically relevant range of masses: , where is the maximum mass star allowed for a particular equation of state. We choose models for our tests that are (approximately) evenly spaced in mass within this range:

 Mi≈1.2M⊙Nstars−iNstars−1+Mmaxi−1Nstars−1, (16)

for . The large dots on each curve in Fig. 1 illustrate the points in the data sets with from three of these equations of state.

We have used these spectral methods to find solutions to the inverse stellar structure problem with mass-radius data sets having , 3, 4 and 5 stellar models. In each case, we have constructed approximate equations of state with non-vanishing spectral parameters, for , …, . We use the Levenberg-Marquardt algorithm for solving the non-linear least squares problem, Eqs. (14) and (15), as described in Ref. Press et al. (1992). We find that this method is extremely efficient at finding solutions that minimize the defined in Eq. (13), given a reasonably accurate initial guess for the values of the parameters and . Our purpose here is to explore the overall accuracy of the spectral inversion method, and is not (at this point) directly aimed at producing an optimal robust method for analyzing real neutron-star observational data. Therefore we use our knowledge of the exact equation of state to provide the initial guesses for and in the test solutions reported here. In particular we use the values of for the stellar models computed from the exact equation of state as initial guesses for . Similarly we use the values of obtained by directly fitting the exact equation of state data as the initial guesses for these quantities. While these initial guesses are not precisely the solutions to the non-linear least squares problem, Eqs. (14) and (15), they are close enough that the Levenberg-Marquardt algorithm easily converges.

We assess the accuracy of our solutions by evaluating the differences between the approximate equation of state produced by the inversion process, and the exact equation of state that was used to construct the data. We measure these differences by constructing the error measure,

 Δ2Nγk = 1NeosNeos∑i=1[log(ϵ(hi,γk)ϵi)]2. (17)

The sum in Eq. (17) is taken over the points, , from the tabulated realistic equation of state. Only the points that lie within the range of densities present in the neutron-star models associated with a particular equation of state are used in this sum. Table 1 lists the values of for each of the 34 realistic equations of state used in our tests. The results reported in Table 1 are for tests that fit the maximum number of spectral parameters, , for each set of data points. Fitting data points gives equations of state approximations with average errors of only a few percent, . Using larger numbers of data points (generally) results in higher accuracy approximations to the equation of state, with average values of , and . Thus, the spectral approach to the inverse stellar structure problem is capable of giving high accuracy measurements of the high density equation of state using only a very small number of data points.

Table 1 also contains two additional measures of the accuracy of our test solutions to the inverse stellar structure problem. One of these,

 ΥN=ΔMRNΔEOSN, (18)

provides another way to measure the error in the approximate spectral equation of state. The quantity in Eq. (18) refers to the error in the approximate equation of state obtained by fitting data, as defined in Eq. (17). The quantity in Eq. (18) refers to the error of the best possible -parameter spectral fit to this particular equation of state. The values of for the equations of state studied here were determined in Ref. Lindblom (2010), and are given in Table II of that reference. The quantity therefore measures the accuracy of the approximate spectral equation of state obtained by solving the inverse stellar structure problem, relative to the accuracy of the best possible approximate -parameter spectral equation of state. Table 1 shows that (almost) all of these measures are of order unity: the approximate equation of state obtained with this spectral inversion method is almost as accurate as the best possible -parameter spectral fit to the equation of state. Table 1 also contains the values of the quantity , defined in Eq. (13), for each of the test solutions found here. These values of are all much less than unity, which shows that the least squares method is doing a good job of minimizing the differences between the model values of and the exact data points .

In addition to the accuracy measures given in Table 1, we have made in depth studies of the errors associated with a few of these solutions to the inverse stellar structure problem. We have selected for closer examination the equation of state whose solution has the smallest error, PAL6 with , the equation of state whose error is the median of the cases studied in these tests, APR3 with , and the equation of state having the largest error, BGN1H1 with . Figure 2 shows the quantity , which measures the fractional difference between the best fit model equation of state and the exact PAL6 equation of state . This figure shows that the errors in these solutions to the inverse stellar structure problem become smaller as the number of data points used in the fits with becomes larger. This figure also shows that the model equations of state do a very good job of approximating the actual equation of state over the entire range of densities that are present in the interiors of neutron stars. Figures 3 and 4 illustrate the analogous error measures for the equations of state obtained from data based on the APR3 and the BGN1H1 equations of state.

These three cases illustrate the range of errors obtained using the spectral inversion method: the best case, a typical average case, and the worst case. We point out that the worst case, BGN1H1, equation of state has a strong phase transition within the neutron-star density range. Non-smooth equations of state of this type are difficult to fit using spectral methods, and convergence in these cases is typically a power law rather than exponential. Many more spectral parameters are therefore needed to approximate these cases accurately. We point out, however, that the values of for the BGN1H1 case are about ten percent, so even in this worst case the spectral inversion method gives a reasonably accurate estimate of the equation of state. The values of the for the BGN1H1 case all have values below 3.5, which shows that while it is difficult to model an equation of state of this type using a spectral fit, the spectral inversion method nevertheless does provide a solution that is comparable to the optimal -parameter spectral fit.

Given an approximate spectral equation of state, , we can use it to compute the complete mass-radius curve for the full range of central enthalpies, . This model mass-radius curve should agree with the exact curve very well, at least at the points used in the inversion process. However, they will not agree everywhere, and the size of the differences is another measure of how well the approximate equation of state agrees with the exact. Figure 5 illustrates the differences between the model masses and the exact masses for the PAL6 equation of state.222Note that the masses in Figs. 5– 7 are compared between models having the same central enthalpy . The central enthalpy of the exact model with need not be exactly equal to the central enthalpy of the best fit model, . Therefore these curves need not have zeros at those points. Figures 6 and 7 make similar comparisons for the APR3 and the BGN1H1 cases. We note that the error measures, , shown in Figs. 57, are somewhat smaller in size than the error measures, , shown in Figs. 24. These error measures provide one more piece of evidence that the spectral solutions to the inverse stellar structure problem are quite accurate.

## V Discussion

In summary, we have developed a new method for solving the relativistic inverse stellar structure problem based on the construction of a spectral expansion of the unknown high density part of the equation of state of the star. The results of our numerical tests of this new method, described in Sec. IV, are quite impressive. Using only two data points, this new method can determine the entire high density part of the neutron-star equation of state with errors (on average) of just a few percent. The addition of more data points (generally) results in higher accuracy approximations. We also show that -parameter spectral approximations to the equation of state determined in this way are almost as accurate as the best possible -parameter spectral approximations. This is quite remarkable. It shows that macroscopic mass-radius measurements are strongly correlated to the properties of the equation of state, and such measurements should therefore allow us (eventually) to measure the high density part of the neutron-star equation of state with great precision.

A close inspection of the results from the various tests summarized in Table 1 reveals a number of anomalies that merit further study. For example, the error measure , defined in Eq. (17), is expected to decrease as the number of spectral parameters is increased, i.e. that . This seems to be true for most of our tests, but there are also a number of exceptions in Table 1. The equation of state FPS, for example, has , , , and . What is going on? Such a sequence of errors would be consistent, for example, with the idea that this particular equation of state is not well represented by these low order spectral expansions, i.e. that these expansions in this case are not yet in the convergent regime. This does not seem to be the case however since the optimal spectral fits to the FPS equation of state do appear to be convergent with these same numbers of spectral parameters, cf. Table II of Ref. Lindblom (2010). Another (more likely) explanation of the anomalous results found in Table 1 is that the minima of found by the Levenberg-Marquardt algorithm for these cases are just local minima and not the desired global minima. An interesting area for further research on this problem, therefore, will be to explore the use of more robust numerical methods for finding global minima of complicated non-linear functions like .

Another interesting direction for future research on this problem will be to explore how robust this kind of solution to the inverse stellar structure problem will be when applied to more realistic data sets. The data used here were idealized in two important ways. First, the mock data used here were supplied with very high precision. Real astrophysical measurements of these quantities will have significant errors. How will measurement errors influence the accuracy of the equation of state that is constructed by these techniques? Second, the mock data used here were chosen to cover uniformly the astrophysically relevant range of neutron-star masses. Real astrophysical measurements will not be distributed in such an orderly way. How will the accuracy of the implied equation of state be affected by different, presumably less ideal, data distributions?

The version of the inverse stellar structure problem studied here is based on the use of mass and radius measurements to determine the high density part of the equation of state. These are not the only macroscopic properties of neutron stars that could potentially be measured. It is not too difficult to imagine for example that the moment of inertias or the tidal Love numbers might be more easily observable for some types of observations. Another interesting direction for future study will therefore be to explore the use of other measurement data, say the mass and Love number (which could be measured using gravitational wave observations of neutron-star mergers), as input for solving the inverse stellar structure problem using the spectral methods developed here.

###### Acknowledgements.
We thank John Friedman, Benjamin Lackey, and Benjamin Owen for helpful discussions about this work, and J.F. and B.L. for providing tables of the various realistic equations of state. A portion of this research was carried out during the time L.L. was a visitor at the Leonard E. Parker Center for Gravitation, Cosmology and Astrophysics, University of Wisconsin at Milwaukee. This research was supported in part by a grant from the Sherman Fairchild Foundation, by NSF grants PHY1005655 and DMS1065438, and by NASA grant NNX09AF97G.

## Appendix A Computing Derivatives of M and R.

This Appendix describes how the derivatives of the total masses and radii are computed with respect to the parameters and . To begin, however, we present a little more detail on how the alternate form of the stellar structure Eqs. (11) and (12) are solved numerically. These equations are,

 dmdh = M(m,r,ϵ,p)≡−4πr3ϵ(r−2m)m+4πr3p, (19) drdh = R(m,r,p)≡−r(r−2m)m+4πr3p, (20)

where the quantities and merely represent the expressions on the right sides. These equations are solved numerically by specifying conditions, , at the center of the star where and then integrating out to the surface of the star where . Like the standard Oppenheimer-Volkoff version of the problem, Eqs. (1) and (2), the right sides of Eqs. (19) and (20), i.e. the functions and , have the ill behaved form there. Consequently it is necessary to start any numerical integration of these equations slightly away from the singular point . The needed starting conditions can be obtained using a power series solution to the equations. The needed power series are given in Eqs. (7) and (8) of Ref. Lindblom (1992), and can be written in the form,

 r(h) = r1(hc−h)1/2+r3(hc−h)3/2 (21) +O(hc−h)5/2, m(h) = m3(hc−h)3/2+m5(hc−h)5/2 (22) +O(hc−h)7/2,

where , , and are given by

 r1 = [32π(ϵc+3pc)]1/2, (23) r3 = −r14(ϵc+3pc)[ϵc−3pc−3(ϵc+pc)25pcΓc], (24) m3 = 4π3ϵcr31, (25) m5 = 4πr31[r3ϵcr1−(ϵc+pc)25pcΓc]. (26)

The quantities , and in these expressions are the energy density, pressure and the adiabatic index evaluated at the center of the star where , , , and .

It will be useful for our least-squares minimization problem to know how the solutions to Eqs. (19) and (20) change as the parameters and are varied. Let denote any one of the parameters: . We wish to derive equations for the derivatives of the solutions to these equations with respect to these parameters: and . It is straightforward to determine the needed auxiliary equations by differentiating, Eqs. (19) and (20) with respect to :

 ddh(∂m∂λ) = ∂M∂m∂m∂λ+∂M∂r∂r∂λ+∂M∂ϵ∂ϵ∂λ (27) +∂M∂p∂p∂λ, ddh(∂r∂λ) = ∂R∂m∂m∂λ+∂R∂r∂r∂λ+∂R∂p∂p∂λ. (28)

The various derivatives , etc. are determined directly from the stellar structure equations, Eqs. (19) and (20):

 ∂M∂m = 8πr3ϵ−Mm+4πr3p, (29) ∂M∂r = −4πr23pM+2ϵ(2r−3m)m+4πr3p, (30) ∂M∂p = −4πr3Mm+4πr3p, (31) ∂M∂ϵ = −4πr3(r−2m)m+4πr3p, (32) ∂R∂m = 2r−Rm+4πr3p, (33) ∂R∂r = −12πr2pR+2(r−m)m+4πr3p, (34) ∂R∂p = −4πr3Rm+4πr3p. (35)

For the case when , the derivatives and are determined from Eqs. (8)–(10). The needed expressions are given by:

 ∂~μ(h)∂γk = ∫hh0[log(h′h0)]keh′dh′Γ(h′), (36) ∂p(h)∂γk = −p(h)∫hh0∂~μ(h′)∂γkeh′dh′[~μ(h′)]2, (37) ∂ϵ(h)∂γk = ∂p(h)∂γkϵ(h)p(h)−∂~μ(h)∂γkehp(h)[~μ(h)]2. (38)

The integrals needed to determine these quantities can be performed numerically with good efficiency and accuracy using Gaussian quadrature. The equation of state does not depend on the parameter , and so . Consequently the equations that determine and in Eqs. (