High-Order Multipole and Binary Love Number Universal Relations

Using a data set of approximately 2 million phenomenological equations of state consistent with observational constraints, we construct new equation-of-state-insensitive universal relations that exist between the multipolar tidal deformability parameters of neutron stars, $\Lambda_l$, for several high-order multipoles ($l = 5,6,7,8$). We confirm the existence of a universal relation between the radius of the $1.4 M_\odot$ NS, $R_{1.4}$ and the reduced tidal parameter of the binary, $\tilde{\Lambda}$, and the chirp mass. We extend this relation to a large number of chirp masses and to the radii of isolated NSs of different mass $M$, $R_M$. We find that there is an optimal value of $M$ for every $\mathcal{M}$ such that the uncertainty in the estimate of $R_M$ is minimized when using the relation. We discuss the utility and implications of these relations for the upcoming LIGO O4 run and third-generation detectors.


Introduction
Due to the constraints imposed by general relativity and causality, there exist quasiuniversal relations between various bulk physical properties of neutron stars (NSs) that are mostly insensitive to the actual equation of state (EOS) of nuclear matter [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Since the nuclear EOS in the high-density regime of NSs is still unknown, these universal relations are a great utility for gravitational wave (GW) astronomy. Universal relations reduce a group of several seemingly independent physical properties to a family characterized by only a few parameters. Ideally, this allows one to break the degeneracies between parameters in the analysis of GW data as well as in waveform modelling.
A robust set of universal relations (called multipole Love relations) holds between the l-th order dimensionless gravitoelectric tidal deformability coefficients of NSs [12], Λ l , which are defined by where C = M/R is the compactness of the NS (here we take G = c = 1) and k l is its l-th order gravitoelectric tidal Love number [15]. The GW waveform of a binary NS (BNS) merger is, quite understandably, highly sensitive to these tidal parameters. How deformable a NS is in a tidal potential affects how its mass ultimately gets distributed during the inspiral of a merger, which, in turn, shapes the GW waveform, especially during the late stages of the inspiral [15][16][17][18][19]. The tidal parameters enter into the waveform at different post-Newtonian orders; however, they are degenerate in the signal [12]. The multipole relations allow this degeneracy to be broken by reducing all of the tidal deformabilities to a family determined by a single parameter. This parameter is always chosen to be the quadrupolar tidal deformability Λ 2 , which is the source of the leading-order finite-size effect in the GW signal and, consequently, is the easiest to measure [12]. Thus, higher-order (l > 2) tidal deformabilities can be expressed through the multipole Love relations as functions of Λ 2 .
The authors and others have demonstrated that the improvements to the accuracy of tidal deformability measurements, to parameter estimation, and to GW modelling offered by the multipole Love relations are significant [12,14,15,20,21] and will become particularly important with the increased sensitivity of upcoming third-generation GW detectors like LIGO III, the Einstein Telescope, and Cosmic Explorer [12,15,[22][23][24][25]. Motivated by these potential improvements, we present entirely new fits to several previously un-fitted high-order multipole Love relations, specifically for l = 5, 6, 7, and 8. Though the finite-size effects of these orders of tidal parameters are currently smaller than measurement error, they will become more measurable with increased sensitivity; hence, faithful GW waveform modelling will need to incorporate them. Previosu studies, such as Flanagan and Hinderer [26] and Damour et al. [27], have discussed the finite-size effects of the l ≤ 4 multipoles.
Zhao and Lattimer [19], De et al. [28] have demonstrated the existence of an intriguing EOS-insensitive relation for BNSs between the radius of the 1.4M NS, R 1.4 and the reduced tidal deformability (also called the binary tidal deformability),Λ. The quadrupolar deformabilities of the individual NSs enter into the GW signal of the merger viaΛ, which is defined asΛ where Λ 2,1 and Λ 2,2 are the deformabilities of the primary and the secondary stars respectively. The quadrupolar tidal Love number k 2 is known to scale roughly as C −1 independently of the EOS [15,29]. According to Eq. (1), this means Λ 2 scales approximately as C −6 . In an apparently analogous fashion,Λ seems to go as (M/R 1.4 ) −6 , where M is the chirp mass of the BNS given by Combining this observation with the definition ofΛ in the manner done by Zhao and Lattimer [19] yields a mostly EOS-insensitive estimate of R 1.4 in terms ofΛ and M that is also mostly insensitive to the binary mass ratio q: The immediate utility of this relation is the ability to produce an EOS-agnostic estimate of R 1.4 from just tidal parameter measurements. This is an alternative to the more involved method of using the universal relation for binaries between the symmetric and antisymmetric combinations of Λ 2,1 and Λ 2,2 [14,30] combined with the relation for individual NSs between Λ 2 and the compactness C [14,31] (a relation which intuitively follows from the definition of Λ l in Eq. (1)). One would first use the symmetric-antisymmetric relation to break the degeneracy between Λ 2,1 and Λ 2,2 and estimate them individually fromΛ, and then use the Λ 2 -C relation and the masses of the binary to extract the radii of both stars. The LIGO/VIRGO analysis of GW170817 is an example of this latter approach [28,32,33]. One need not appeal to universal relations to estimate stellar radii, however. Instead, one could perform an inference of the EOS directly using a parametric representation of the EOS, as was also done in the LIGO/VIRGO analysis [32], or using a much more sophisticated nonparametric representation, as described in Essick et al. [34].
It is an appealing question, then, whether this relation can be extended using the radius of a NS with a generic mass M, R(M) = R M . A R M -Λ relation would allow one to use measurements of tidal parameters and M to place robust constrains on R M directly without the need for a more complicated procedure. Hence, our motivation in this work is to provide a phenomenological study of the R M -Λ relation. We look at the relation for several values of M. For a given M, we compute fits to the relation for twelve fixed values of M between 0.9M and 1.4M . We then generalize the fit for all M ∈ [0.9M , 1.4M ] by interpolating the fitting parameters as functions of M. Fitting to the relation from across a vast set of phenomenological EOSs incorporates the effects of higher-order terms that are dropped when one analytically derives the expression in Eq. (4) as was done in [19]. Eq. ( 4) assumes that the R 1.4 -Λ relation (and, by extension, the R M -Λ relation) is only linearly dependent on M, i.e. that M simply scales the relation but does not change its dependence onΛ. A phenomenological study permits us to observe directly the effect changing M has on the relation.
The outline of this paper is as follows. In Sec. 2, we describe the parameterization scheme and algorithm by which we generate our phenomenological EOS data and the statistics for our analyses. In Sec. 3, we present the fitting parameters of the high-order multipole Love relations, followed in Sec. 4 by an phenomenological analysis of and fits to the R 1.4 -Λ relation as well as to the general R M -Λ relation. We also discuss the implications of these new fits to GW waveform analysis for the LIGO O4 run. A concluding summary is given in Sec. 5.

Methods
We parameterize the space of all possible EOSs consistent with theoretical calculations and astronomical observations using the piecewise polytropic interpolation developed in Read et al. [35], with the only modification being that we allow the transition densities ρ 1 and ρ 2 to vary. We then generate random piecewise EOSs using a Markov chain Monte Carlo (MCMC) algorithm, with the basic summary as follows. For a given candidate EOS, the algorithm first computes a series of solutions to the Tolman-Oppenheimer-Volkoff (TOV) equation using the publically available TOVL code described in Bernuzzi and Nagar [36] and Damour and Nagar [17], and then accepts the EOS if and only if it satisfies three weak physical constraints: causality of the maximum mass NS is preserved (i.e. the maximum sound speed c s is less than the speed of light c below the maximum stable central density), 2.
the maximum stable mass of a non-rotating NS, M max , is greater than 1.97M , and 3.
The full details of parameterization and the MCMC algorithm can be found in Godzieba et al. [14]. With this scheme, we generate a set of 1,966,225 phenomenological EOSs.
To study the multipole Love relations, for each EOS in our data set, we solve the TOV equation for sixteen evenly spaced central densities between ρ c = 3.09 × 10 14 g/cm 3 and the maximum stable central density of that EOS, and then extract Λ l for l = 2 through l = 8 from each solution.
To study the R M -Λ relation, we follow a similar procedure. First, we choose a fixed value of M. Next, for each EOS in a random sample of a quarter of all EOS in the data set, we generate twenty random binary NSs (BNSs). We uniformly sample the binary mass ratio q = m 2 /m 1 (where m 2 ≤ m 1 ) on the interval 1/2 ≤ q ≤ 1. This range is not intended to represent the complete range of values that q could take in Nature, but rather simply to capture the general behavior of q based on observational and theoretical considerations. Observations of the most massive known pulsars indicate that M max 2M [37][38][39][40][41][42], and the analysis of GW170817 suggests that M max 2.3M [43][44][45][46][47]; though, as we await upcoming precision measurements of millisecond pulsar radii by NICER, we cannot as of yet categorically rule out the possibility of extreme EOSs with M max < 2.5M [48]. Meanwhile, the least massive known pulsar has a mass of 1.17M [49], and, depending on the true nuclear EOS, the minimum stable gravitational mass, M min , could be as low as 1.15M [50]. Hence, q ≥ M min /M max ≈ 1/2, and the vast majority of BNSs, being far from either mass extreme, will fall well within this range.
Each q is then converted into the actual binary masses m 1 and m 2 using the value of M and Eq. (3): The TOV equation is then solved with the corresponding EOS for NSs with these two masses, and Λ 2 is extracted from both solutions to computeΛ. We apply this procedure for twelve different values of M between 0.9M and 1.4M . TheΛ values are then plotted versus R M for eight different values of M. These R M values are pulled from our EOS data set.

High-Order Multipole Relations
From our phenomenological EOS data set, we compute 21,994,104 valid individual NS solutions to the TOV equation as our statistics for analyzing the multipole Love relations. Λ 3 , Λ 4 , Λ 5 , Λ 6 , Λ 7 , and Λ 8 are plotted against Λ 2 in Fig. 1, and one can appreciate the universality of each relation across a vast range of scales. (We observe that all the intersections between the six curves lie around Λ 2 ∼ 100, but we are not sure why this is the case.) As in the authors' previous work [14], we employ a fitting function of the form which is an extended version of the fitting function originally used by Yagi and Yunes [11]. The fitting parameters a = {a k } for each relation are given in Table 1.
In Fig. 2, we show the 68%, 95%, and 99.7% relative error of each fit. For each line in the error plot, the corresponding percentage of data points lie below it. We restrict our attention to the domain 1 < Λ 2 < 10 4 , as this is the range of Λ 2 most relevant to current LIGO measurements. The estimate error of each Λ l over this range stays mostly flat with a slight downward trend. (The small ripples that can be seen in the error plots over this range are simply artifacts of how the distribution of NS solutions was computed.) The universality of the multipole relations weaken gradually as l increases, as can be seen in the increasing thickness of the distributions in Fig. 1. This then increases the maximum estimate error of Λ l for larger l despite the faithfulness of each fit to the shape of the corresponding relation (see Fig. 1). While 95% of estimate errors are smaller than ∼7% for the Λ 3 -Λ 2 relation, 95% are only smaller than ∼50% for Λ 8 -Λ 2 .
The phase of a GW in waveform modelling is affected by the highest order out to which one carries finite-size corrections. 1 We demonstrate this with a baseline model of a binary with m 1 = m 2 = 2.7M and Λ 1 = Λ 2 = 1000 using the spin-aligned effective-one-body waveform model TEOBResumS [18]. Often when universal relations are not employed, all finite-size effects are dropped except for the leading-order (l = 2) effect. In the baseline model, just the l = 2 correction alone contributes a phase difference of 36.7 radians compared to a waveform model with no tidal corrections. Further corrections from the l = 3 and l = 4 effects using the Λ 3 -Λ 2 and Λ 4 -Λ 2 relations respectively incur an additional 2.89 radians. Finally, including the l = 5, 6, 7, and 8 corrections using the relations given in this work adds 0.02 radians of dephasing on top of that. (The dephasing between the l ≤ 8 waveform model and models with fewer corrections is plotted in Fig. 3 as a function of time. For all models, most of the dephasing is accumulated in the last 5 milliseconds before the merger.) Combined, the l > 2 corrections contribute 2.91 radians of dephasing. This demonstrates the importance of the multipole Love relations for faithful waveform modelling.
The dephasing of the l > 4 corrections are currently smaller than GW detector uncertainties, but this could only have been known after fitting to the l > 4 multipole relations. Additionally, with the greater sensitivity of future detectors, the l > 4 finite-size effects will start to come into view. The order out to which one should carry finite-size corrections in the waveform analysis of actual GW data is dependent on several factors (the EOS model, the signal-to-noise ratio of the merger, etc.); however, in general it is recommended that corrections up to l = 4 be included in the analysis of data from current detectors [12,14,27].

R M -Λ Relation
We analyze the R 1.4 -Λ relation at twelve different fixed values of the chirp mass M, which are given in Table 2. We compute between 750,000 and 1,000,000 valid individual binaries for each value of M. Several example plots of the relation are shown in Fig. 4. The relation's dependence on the binary mass ratio q is illustrated by the coloring of the points in these plots. Each point in the plot represents a BNS. Points with smaller values of q are plotted on top. An important conclusion to draw from these plots is that the relation does not depend upon both stars having the same radius [19]. For M 1.25M , the relation  remains fairly tight for all values of q. Further, for M 1.1M 2 , the relation actually becomes tighter as q decreases (i.e. as the radii of the two stars differ more and more), which can be understood by considering the definition ofΛ in Eq. (2). For fixed R 1.4 , the range of possible valuesΛ can take is constrained by q. When q = 1, the masses of the binary can span the range from the minimum to the maximum mass, M min ≤ m 1 = m 2 ≤ M max . Hence, min (Λ 2 ) ≤ Λ 2,1 = Λ 2,2 ≤ max (Λ 2 ), and min (Λ 2 ) ≤Λ ≤ max (Λ 2 ). As q decreases, the bounds for both m 1 and m 2 shrink and no longer overlap, causing the same to happen for Λ 2,1 and Λ 2,2 . This, as we see, also shrinks the bounds onΛ. Thus, we expect the relation to tighten as q decreases.
We construct a fitting function for the R 1.4 -Λ relation by considering a slightly generalized form of Eq. (4): where, in this case, M = 1.4M . Here the proportionality constant α and the inverse exponent β are the fitting parameters and, consequently, will be dependent on M. These fits are also shown in Fig. 4. The fitting parameters for all values of M are given in Table 2 Table 2, are shown in Fig. 6. The fitting parameters p = {p k } and q = {q k } for α(x) and β(x) are given in Table 3. What is interesting is that the inverse exponent β is not monotonic. Rather, it has a minimum at M = 1.1661M . A possible contributor to this effect is the decrease in the variety of possible binaries as M increases. The maximum value M could take for a given EOS is found by letting m 1 = m 2 = M max in Eq. (3), which yields M max = 2 −1/5 M max . For 1/2 ≤ q ≤ 1, m 1 and m 2 are bounded by Hence, as M increases, the relation becomes gradually dominated by (1) EOSs with M max ≥ 2 1/5 M and (2) only those BNSs from each EOS that lie in the increasingly narrow range in Eq. (9). This decrease of BNS variety could play a role in the non-monotonic behavior of β(x). One could, of course, consider more generally the relation betweenΛ and the radius of a NS with some mass M, R M . We pursue this thought by looking at R M -Λ for M/M = 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, and 2.14. Just as for the R 1.4 -Λ relation, we utilize the fitting function in Eq. (4), and then find α and β as functions of x using Eq. (8). In Table 3 Fig. 7. This reveals that there is an optimal M for each M such that R M is maximally constrained by the R M -Λ relation at that M. Thus, for example, a chirp mass of M ≈ 1.05M would yield the best estimate of R 1.4 , while a chirp mass of M ≈ 1.4M would yield the best estimate of R 1.8 . Further, there appears to be a linear dependence of the optimal M on M; however, a wider range of M would need to be considered to confirm this. The change in the variety of BNSs as M increases, as previously described, may contribute to this. At larger M, the relation becomes dominated by larger mass NSs; thus, the relation may become more sensitive to the radii of larger mass NSs as M increases.
The R M -Λ relation, then, allows one to use any binary to place a robust, EOS-agnostic constraint on R M using justΛ and M. This offers great prospects for the upcoming LIGO O4 run and for third-generation detectors. The O4 run expects to see 10 +52 −10 detections within a search volume of 1.6 × 10 7 Mpc 3 year [51]. Every BNS detection can be transformed into a maximum constrain on some R M . However, even the weaker constraints afforded by R M -Λ are of still great utility. Just 10 weak constraints on R 1.4 using the R 1.4 -Λ relation will yield a reliable value for R 1.4 . Further, a reduction in statistical uncertainty thanks to increased sensitivity improves the effectiveness of universal relations, as, for example, the systematic errors of fits to multipole relations are generally smaller than statistical uncertainty [12].

Conclusions
We supplement the tool set of GW analysis and waveform modelling by presenting entirely new fits to several universal relations between high-multipole-order dimensionless gravitoelectric tidal deformabilities Λ l and to the universal relation for BNS between the radius of the 1.4M NS, R 1.4 , and the reduced tidal deformabilityΛ. We compute these utilizing a data set of nearly two-million phenomenological EOS sampled from across a broad parameter space using an MCMC algorithm.
First, we present fits to the multipole relations. Previous fits [12,14] had been made to just the Λ 3 -Λ 2 and Λ 4 -Λ 2 relations. We extend the library of fits by looking at the Λ 5 -Λ 2 , Λ 6 -Λ 2 , Λ 7 -Λ 2 , and Λ 8 -Λ 2 relations. The tightness of the relations weakens as l increases. Consequently, though the fits are faithful to the shapes of the relations, the maximum estimate error of the fits increases to the order of 50% for Λ 8 . The inclusion of the finite-size effects of the l < 4 multipoles in waveform analysis can incur as much as 0.02 radians of dephasing compared to including only the l ≤ 4 effects. Collectively, the l > 2 effects contribute as much as 2.91 radians of dephasing, and it is recommended that finite-size corrections for l > 2 multipoles be included in the analysis of GW data wherever they are at least comparable to detector uncertainties. The full usefulness of these l > 4 relations in GW data analysis will be realized with the increased sensitivity of the upcoming thirdgeneration GW detectors like LIGO III [22], the Einstein Telescope [23,24], and Cosmic Explorer [25], as the finite-size effects of these multipole orders are currently smaller than measurement error.
Next, we analyze the R 1.4 -Λ relation. The original derivation of the relation [19] yields an expression, given in Eq. (4), that is linearly dependent on the chirp mass M of the BNS. Fitting the relation for different fixed values of M reveals any nonlinear dependence the relation may have on M and allows us to compute an expression that more accurately estimates R 1.4 . We do this for twelve different values of M between 0.9M and 1.4M , using the fitting function in Eq. (7). We then interpolate the fitting parameters to all M ∈ [0.9M , 1.4M ] by fitting them as functions of M. The accuracy of the estimate of R 1.4 for any value of M is found to be quite good. 95% of the estimates are within ∼5% of R 1.4 .
We then consider a generalized form of the relation R M -Λ for a generic NS mass M. We perform the same analysis as for the R 1.4 -Λ relation for seven other values of M. We find that the level of uncertainty in the estimate of R M depends on both M and M. There is, in fact, an optimal value of M for each M such that R M is maximally constrained by the relation at that value of M. Therefore, this relation will be an excellent tool for combining the results from multiple GW detections of BNSs into constraints on NS radii.
The parameter space of possible EOS explored by our MCMC algorithm to compute our EOS data set can be further restricted with the inclusion of possible future LIGO/Virgo/KAGRA O4 constraints, laboratory constraints, such as those from heavyion collisions [52] and PREX [53], X-ray burst observations from NICER [54][55][56], and by combining the phenomenological EOSs with results from pQCD calculations [57].
Author Contributions: D.G. generated and analyzed the TOV data, D.G. and D.R. interpreted the results and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.