Finite tidal effects in GW170817: Observational evidence or model assumptions?
Wolfgang Kastaun Frank Ohme











Max Planck Institute for Gravitational Physics
(Albert Einstein Institute), Callinstr. 38, D-30167 Hannover, Germany
Leibniz Universität Hannover, D-30167 Hannover, Germany








After the detection of gravitational waves caused by the coalescence of compact 
objects in the mass range of neutron stars, GW170817, several studies have searched 
for an imprint of tidal effects in the signal, employing different model assumptions. 
One important distinction is whether or not to  
assume that both objects are neutron stars and obey the same equation of state.
Some studies assumed independent properties. Others assume a universal equation
of state, and in addition that the tidal deformability follows certain 
phenomenological relations. An important question is 
whether the gravitational-wave data alone constitute observational evidence for finite tidal effects.
All studies provide Bayesian credible intervals,
often without sufficiently discussing the impact of prior assumptions, especially in the case of lower limits on the neutron-star tidal deformability or radius.
In this article, we scrutinize the implicit and explicit prior assumptions made 
in those studies. Our findings strongly indicate that existing lower credible
bounds are mainly a consequence of prior assumptions combined with information gained about the system's masses. Importantly, those lower bounds are typically not informed by the observation of tidal effects in the gravitational-wave signal. Thus, regarding them as 
observational evidence might be misleading without a more detailed discussion. 
Further, we point out technical problems and conceptual 
inconsistencies in existing studies.
We also assess the limitations due to systematic waveform model uncertainties
in a novel way, demonstrating that differences between existing models are not guaranteed to be small enough for an unbiased estimation of lower bounds on the tidal deformability.
Finally, we propose strategies for gravitational-wave data analysis designed
to avoid some of the problems we uncovered. 




04.25.dk,  
04.30.Db, 
04.40.Dg, 
97.60.Jd, 


Finite tidal effects in GW170817: Observational evidence or model assumptions?
    F. Ohme
    September 2019; updated September 2020
==============================================================================







§ INTRODUCTION


Among known astronomical objects, the highest matter densities occur
within neutron stars (NSs).
The behavior of cold matter at supranuclear 
densities cannot be studied in laboratories on Earth. Therefore, 
fundamental properties such as the equation of state (EOS) that relates 
pressure, temperature, and density, are essentially unconstrained by experiments. 
Theoretical models calibrated by properties of atomic nuclei 
can be used to extrapolate the EOS of neutron-star matter, but face 
considerable theoretical uncertainties.
Astrophysical observations of NSs can be used instead to infer
EOS constraints, based on the EOS impact on NS structure and dynamics. 

Merging NS binaries in particular are an excellent environment to probe the EOS, 
since the gravitational wave (GW) signal emitted before, during, and after merger
carries signatures of the EOS, as do electromagnetic (EM) counterparts
given by kilonova and short gamma ray burst.
The first multimessenger observation  
<cit.> 
of a binary neutron star (BNS) merger event, GW170817,
indeed opened a wealth of opportunities.

The most direct observable signature of the EOS is contained in the GW signal
emitted during the last seconds before the merger.
While a neutron star and another compact object spiral around each other,
the neutron star is deformed
in the tidal field of the companion.
This deformation modifies the dissipation of energy and angular momentum by
GW and thereby affects the decay of the orbit.
Measuring the inspiral frequency evolution therefore carries imprints of the 
EOS <cit.>, in addition to the constituents' masses and 
angular momenta.

For binary neutron stars, the frequencies during this phase fall within the 
sensitive range of current ground-based interferometers, such as the 
Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO)  <cit.>
or the Virgo detector  <cit.>.
In contrast, merger and postmerger GW emission is dominated by 
frequencies >1 kHz, for which the current detection range 
is considerably smaller, making observations of this phase unlikely (compare 
<cit.>).

For the first observation of GWs from a binary NS merger, GW170817 <cit.>, 
several studies have already presented constraints on the tidal parameters that enter the GW 
signal  <cit.>. 
All of those provide robust upper limits on the tidal effects.
Although more difficult to measure, lower bounds are provided as well.

In this paper, we revisit the methods used in  
<cit.> to derive such constraints.
We scrutinize the implicit and explicit assumptions involved, in order to decide
if the published lower limits can be interpreted as direct observational 
evidence of finite tidal effects, or rather constitute inferences based on
model assumptions and observational constraints on the masses.
The latter would be the case when the same (or very similar) limits are obtained regardless 
of whether the real source was a neutron star or black hole (BH) binary.

We emphasize that we do not argue against combining observational data and prior knowledge,
e.g., from nuclear physics, in order to provide tighter constraints on the EOS.
Before interpreting such constraints as observational evidence, however, 
the contribution of model assumptions needs to be accounted for. In any case,
if the model assumptions include ambiguous ad-hoc choices or inconsistencies,
then resulting constraints are only meaningful to the extent that they are dominated by the 
observational data and not the prior knowledge. 
We review prior choices with regard to the above.


Before explaining the main challenges, we briefly review the 
Bayesian framework that is typically used to infer the properties of binary GW 
sources (see <cit.>). The binary's vector of parameters, θ  (composed of the binary NS 
properties such as masses, distance, inclination, 
etc., and also EOS-related parameters) is determined in the form of its posterior 
probability density, 

    p(θ| d, I) = P(d |θ, I)/P(d | I)  p(θ| I),

given the observed data d and any 
additional model assumptions I. Bayes' Law, eq:Bayes_law, relates the 
probability density distribution on the left-hand side to the product 
of the likelihood of the observed the data given the source parameters, 
P(d |θ, I), and the prior probability density of encountering a system with 
those parameters, p(θ| I). The likelihood is calculated by comparing 
gravitational waveform models with the data, assuming stationary Gaussian 
noise. The prior represents any knowledge or belief that existed before the 
data were obtained.

The main conceptual challenges all such Bayesian EOS analyses of binary 
NS GW data face are the following:

 
  * Parameterization of the unknown EOS. GW models depend on two 
parameters related to the tidal deformability of each object in the binary. 
If a common but unknown EOS is assumed, relating those two numbers is nontrivial. 
In Secs. <ref> and <ref>, 
we discuss model assumptions  <cit.> (often 
referred to as “universal relations”) that have been used to correlate 
masses and tidal deformabilities of NS binaries, analyze their range of 
applicability as well as pathologies arising from their use, and suggest a simpler
alternative.

  * Using the posterior to claim observational evidence for tidal effects.
Even for a completely unconstraining measurement, Bayesian methods provide
a posterior distribution, which then only reflects the chosen prior distribution.
To assess the contribution of observational data, 
one has to compare posterior and prior. 
This becomes nontrivial in a multidimensional parameter space.
If the marginalized posterior for tidal effects is constrained in comparison
to the prior, it still does not signify a direct measurement of tidal effects. 
The reason is that the data might only constrain the NS masses. If the prior 
distribution for tidal effects depends on the masses, the posterior is then influenced
even without any constraint of tidal effects. We therefore separate
the indirect impact of observational constraints on quantities that are unrelated to 
tidal effects from the impact of those that are directly related.
For this, we discuss idealized cases such as fixed mass ratios.

  * Ambiguity of the prior probability for the EOS-related quantities. 
The choice of prior is a well-known issue in Bayesian analysis. 
There are no other observations
that could be used as prior knowledge on tidal deformability, and the 
requirement to specify a prior distribution can only be satisfied by ambiguous 
ad-hoc choices.
We examine the shape of priors used for published EOS constraints
in detail and discuss implications.

  * Credible bounds do not imply statements on likelihood of particular models.
Depending on the shape of the posterior,
a specific model with parameters outside the credible bounds can be more likely
than another model within the bounds. Even the maximum likelihood parameter
can be located outside a given credible interval. It is always possible to 
compute a two-sided credible interval, which, on its own, is therefore not 
useful as evidence against specific models. 
For this, Bayesian model selection is more suitable. A recent model selection 
study on GW170817 <cit.> did not find evidence for or against 
finite tidal effects. In this work, we do not compute Bayes factors, but provide 
arguments, based on features of the posterior and prior distributions, 
against attributing the lower credible bounds to observational evidence.


  * Theoretical uncertainty of waveform model.
Computing the likelihood requires a theoretical model of the gravitational 
waveform for given parameters. The systematic error of the latter is currently
not known, but estimated by comparing different theoretical models 
<cit.>.
Further, estimating the resulting systematic error of credible bounds based on 
statistical analysis of a given event is complicated by the presence of 
the waveform-independent statistical errors.
In sec:syserr we provide a new perspective on estimating the potential impact 
of systematic waveform errors. The order of magnitude we find 
for the systematic errors 
indicates that lower bounds on tidal parameters should 
be interpreted with great caution.


The issues discussed in this work are also relevant when planning the data 
analysis strategy for future events.
This involves a choice of how many model assumptions to incorporate already during the GW 
data analysis, and which quantities to provide for comparison to nuclear physics models.
In Sec. <ref>, we propose a strategy that keeps the division 
line between data analysis and neutron-star modeling as close to the measurable 
quantities as possible, which simplifies the interpretation of results and
allows one to study the consequences of different model assumptions without repeating 
computationally expensive parameter estimation studies.
Employing a reasonable strategy will be important to facilitate 
the exciting synergy of GW astronomy and nuclear physics as more observations 
become available.






§ BINARY NEUTRON-STAR PARAMETERS





 §.§ Masses


The evolution of the GW signal during the inspiral depends most strongly 
on a combination of the individual gravitational masses
M_1 and M_2 called chirp mass,

    M_c    = (M_1 M_2 )^3/5/( M_1 + M_2 )^1/5.

For a BNS merger, the chirp mass can be determined very accurately by current GW 
detectors. We therefore treat it as fixed in this article. We 
specifically use the value M_c=1.186 inferred for GW170817 <cit.>.  
The mass ratio q = M_2/M_1 is less well constrained. 
For the given chirp mass, we can express the individual NS masses as a function 
of the mass ratio,

    M_1(q)    = M_c ( 1 + q )^1/5 q^-3/5, 
    
    M_2(q)    = M_c ( 1 + q )^1/5 q^2/5.

We focus on the simple case of vanishing NS spins in this work, 
although they can have a significant influence on the GW signal.
The particular challenges that we discuss are not reduced when allowing 
for spins, while ignoring spins considerably simplifies the discussion.
We also make the common assumption of negligible orbital eccentricity.





 §.§ Effective tidal deformability



In the context of GW signal, the difference between NS and BH cases is 
how the objects' quadrupole moments react to the tidal forces caused by their 
companion. This is expressed by the dimensionless tidal deformabilies
Λ_1, Λ_2 of the two objects, which are 0 for BHs (see <cit.> however).
Since we neglect the NS spins, finite size 
effects are given by the tidal deformability as a function of 
gravitational mass alone, Λ(M), which depends on the EOS.
The correction to the point-particle GW signal mainly depends on 
a combination (see <cit.>) called effective tidal deformability,

    Λ̃   = 16/13(1 + 12 q ) Λ_1 + (q + 12 ) q^4 Λ_2/(1 + q )^5.

In this work, we use the rescaled quantity

    Λ̅   = Λ̃ w_f^-1(q),  
    
    w_f(q)    = 16/131 + 12 q + (q + 12 ) q^4/(1 + q )^5.

The interpretation of Λ̅ is convenient because
for a BNS, it results from a linear interpolation of Λ(M) between
the masses M_1 and M_2, evaluated at a mass M̅ defined by

    M̅(q)    = w_1(q) M_1(q) + w_2(q) M_2(q), 
    Λ̅(q)    =  w_1(q) Λ(M_1(q)) + w_2(q) Λ(M_2(q)).

For mass ratios close to unity, Λ̅≈Λ(M̅).
The weights depend only on the mass ratio, as

    w_1(q)    = 1 + 12 q/1 + 12 q + (q + 12 ) q^4, 
    
    w_2(q)    = (q + 12 ) q^4/1 + 12 q + (q + 12 ) q^4.

where w_1 + w_2 = 1. We follow the convention q<1, and therefore 
w_1≥ 1/2 and w_2 ≤ 1/2.
Note M̅(q) is monotonically decreasing, which allows us to express Λ̅ as 
a function of M̅. 


For the BNS case with both NS obeying the same EOS, Λ̅ is shown in the upper panel of 
fig:lambda_m_vs_M_m_all for a representative set of EOS. Note that the function is defined 
only up to a maximum M̅ (minimum q) for which M_1 exceeds the maximum 
allowed mass for a nonrotating NS for the given EOS. For the case of a mixed binary, we can 
compute Λ̅ by setting either Λ_1 or Λ_2 to 0 in eq:lambda_eff.
The resulting differences are shown in the lower panel for fig:lambda_m_vs_M_m_all.




Although only Λ̅ can be confined by the GW signal, the difference of the 
tidal deformabilities of the two NSs is an important property of the EOS as well.
In the following, we express it in terms of a dimensionless parameter

    S̅   = ln(Λ_1) - ln(Λ_2)/ln(M_1) - ln(M_2)
            = ln(Λ_2/Λ_1)/ln(q).

For moderate mass ratios, S̅ has an intuitive interpretation 
in terms of the slope of Λ(M),


    S̅≈dln(Λ)/dln(M)(M̅).


In the context of waveform models, it is more common to use a quantity 
δΛ̃, which enters at higher post-Newtonian order than 
Λ̃ and has minor impact on the waveform <cit.>.
For fixed Λ̅ and q, the quantity δΛ̃ 
becomes a function of S̅, which therefore has minor impact as well. 
Unless noted otherwise, we discuss the idealized case
that the GW signal depends on
M_c, q, Λ̅ but not S̅.





 §.§ Universal relations


In order to make statements about the individual tidal deformabilites,
the measurement of Λ̅ is not sufficient; the slope S̅ is 
required as well. Since it does not enter the GW signal, it cannot be 
measured. One way to infer properties of the individual NSs regardless of this fact
is to use model-dependent assumptions about S̅.

As proposed in <cit.>, one can express the effective tidal 
deformability Λ̃ as an EOS-independent function of the average 
deformability 
and the mass ratio. 
The residual error (computed in <cit.>)
of the expression in comparison 
to the results for a few nuclear physics EOS models
is small enough to motivate the term universal relation.


In detail, the relations are written as Λ_a(Λ_s, q), where
Λ_s = ( Λ_1 + Λ_2 ) / 2 and
Λ_a = ( Λ_2 - Λ_1 ) / 2.
The function Λ_a is expanded in Λ_s and q as 

    Λ_a    = 
      Λ_s F_n(q) 
      1 + ∑_l=1^3 ∑_k=1^2 b_lk q^k Λ_s^-l/5/
            1 + ∑_l=1^3 ∑_k=1^2 c_lk q^k Λ_s^-l/5, 
    
    F_n(q)    = 
      1 - q^10/(3-n)/ 1 + q^10/(3-n) .

Those relations are used in <cit.> for data analysis of GW170817. 
The coefficients n, b_ij, c_ij, given in <cit.>,
were obtained by fitting to a set of nuclear physics EOS. 
We noticed that the coefficients differ from those given in 
<cit.> (after taking into account that the 
definitions differ by inclusion of a scaling factor).
However, the qualitative statements we make on universal 
relation usage are equally valid for both.

The set of EOS used for calibration in 
<cit.> consists 
of SLY <cit.>, LS220<cit.>, 
Shen <cit.> (here:STOS), AP3 and AP4 <cit.>, 
WFF1 and WFF2 <cit.>, ENG <cit.>, 
MPA1 <cit.>, MS1 and MS1b <cit.>.
We note that the AP3, AP4, WFF1, WFF2, and ENG EOS start to violate causality before 
reaching the central density of the maximum mass nonrotating NS solution. 
It is also worth pointing out that the set of EOS models used for calibration has 
maximum NS masses compatible with the observation of PSR J0348+0432 
<cit.> (although
for the WFF1 EOS, causality is violated at this mass). Since low tidal deformabilities
typically occur near the maximum mass model and increase with decreasing mass,
nuclear physics EOS candidates compatible with PSR J0348+0432 do not yield low 
deformabilities in the mass range measured for GW170817 (for the EOS set above,
Λ>61 in the mass range inferred for the low spin prior 
in <cit.>; compare Fig. 1 in <cit.>).
On the other hand, strong phase transitions are not considered and might alter the 
behavior of tidal deformability.


The universal relations Λ_a(Λ_s, q) are shown in fig:Sl_vs_Lambda_m,
in terms of our Λ̅, S̅ parametrization.
The relations are evaluated for our fixed fiducial chirp mass, and various mass ratios in 
the range 0.4… 1.
For comparison, we plot the values Λ̅, S̅ for a representative set of EOS. 
As one can see, none of the EOS reach low values of Λ̅ for the mass range of the 
detected event. Using the universal relations at lower values of Λ̅
is self-contradictory, since the EOS used to derive the relations explicitly exclude
such values. We note a subtle point: the models used in 
<cit.> to fit the universal relations do include lower values 
of Λ̅. The reason is that the mass range of the models for the fit 
is not restricted. One should however keep in mind that the resulting fit around 
any given Λ̅ is only valid within a range for M̅ allowed by the 
EOS models. In other words, there might be EOS leading to low Λ̅ for the 
detected mass range, but these EOS must then be very different from those used in the fit. 
In turn, the error of S̅ from the universal relations is completely unknown in that 
case.

Besides the conceptual problem above, the relations used in <cit.>
exhibit some technical problems.
For low values of Λ_s, the expression eq:univ has 
a pole and ceases to be meaningful. The corresponding slope S̅ diverges.
The divergence leads to gaps in the possible 
values of Λ̅. However, the divergence only affects very small values.
More important is the fact that Λ_1 computed from the universal relation 
crosses 0 (we recall that negative tidal deformabilities are unphysical 
<cit.>). At this point S̅ diverges since it is defined in a way
that does not allow negative deformabilities.
This happens below some critical Λ̅ that depends on the mass 
ratio. Even before Λ_1 crosses 0, the slope S̅ greatly exceeds the
range given by the EOS models, as shown in fig:Sl_vs_M_m.
Assuming the universal relations to hold universally thus excludes
small effective tidal deformabilities. The excluded region is shown in 
fig:lambda_m_vs_M_m_all as a function of M̅.
Although not relevant for the case of GW170817, we note that for mass ratios
q≲ 0.5, low values of Λ̅ are excluded because the universal
relation develops a minimum. For even more extreme mass ratios, a
region at large Λ̅ becomes invalid because Λ_1<0.


One could argue that the excluded regions at small Λ follow from 
the universal relations, which hold to good accuracy for known nuclear 
physics EOS. However, the EOS used to derive the relations already exclude
even larger Λ, as discussed earlier. The exact extent 
of the region excluded by negative Λ_1 is just the result of 
extrapolating a fit and has no physical meaning. 





A similar set of universal relations is used in <cit.>.
The set of EOS used for calibration consists of many piecewise 
polytropic EOS, each of which leads to maximum NS masses above 2 M_⊙.
For the mass range of the detected event,
the tidal deformability of NS following those EOS is approximately
proportional to q^6. This leads to a parametrization
Λ_1 = Λ_0 q^3, Λ_2 =  Λ_0 / q^3, where
Λ_0 is a free parameter.
In our notation, the slope parameter is S̅=6, independent of Λ̅
and q.
The given slope parameter is also shown in fig:Sl_vs_Lambda_m.
The differences to NS models with nuclear physics EOS are comparable to those 
found for the universal relation used in <cit.>.
In contrast to <cit.>, the relation used in 
<cit.> does not exhibit any divergence or excluded regions 
at low values of Λ.




It has become customary to express constraints on the EOS from NS observations 
in a mass-radius diagram. However, the radius does not enter the gravitational waveform
models at all and therefore it cannot be measured.
In <cit.> and  <cit.>, the NS tidal
deformabilities were therefore converted to radii using another type of “universal” relation,
which relates tidal deformability and compactness C=M/R of the NS.

The universal Λ-C relation used in <cit.>
is given in <cit.>, while 
<cit.> introduces a different one.
Both variants are shown in fig:univ_compactness. 
For comparison, we 
plot values for our standard set of EOS at two masses spanning the range of the detected 
event.
The universal relation C_u(Λ) contains the same conceptual trap as the relation 
Λ_a(Λ̅, q): for a given NS mass, it is only meaningful in the range 
covered by the set of EOS 
used for the fit. Outside this range, none of those models has the correct mass. 
For a NS with given mass and low tidal deformability, one cannot derive any expectation 
for the compactness from nuclear physics EOS if the very same set of EOS firmly excludes 
NS models with those Λ, M.
We note that outside the range supported by NS models, the two relations 
differ significantly, since the relation used in <cit.> was only 
designed to fit the mass range of GW170817.


For each of the two masses shown in fig:univ_compactness, 
the range of EOS spans only a limited range in Λ. For the example at mass
1.14  M_⊙ (blue markers), the universal relation prediction for the compactness becomes
increasingly unreliable below
Λ≲ 500, the lower limit of the range allowed by the EOS set.
This corresponds to a radius of 
≈10 
when using the universal relation. 
NSs of same mass but lower tidal deformability
would require an EOS very different from those used to construct the universal
compactness relation. A possible cause for such qualitatively different behavior 
might involve phase transitions or quark matter, none of which are considered in 
the EOS sets used for calibration in <cit.>.

In conclusion, assuming that the EOS of an observed NS falls in the 
range of EOS used to calibrate the universal compactness relations already 
leads to lower limits of radius and tidal deformability based on the measured 
mass range alone.
Inserting any measured tidal deformability below this limit
in the universal relations to obtain a radius is contradictory,
because the existence of such a NS invalidates the foundations of the relation.

A consistent use of universal relations in statistical analysis would require 
specifying a mass-dependent validity range defined by the set of EOS used to calibrate the 
relations, and to limit the prior support to valid regions in parameter space.
Most likely, the resulting lower limits for GW170817 would then be close to 
the model-dependent ones already resulting from the validity region and 
the measured mass range.











§ APPLICATION TO DATA ANALYSIS







 §.§ Ideal measurements
  

In the following we study the constraints given by an exact measurement 
of the waveform. Although real measurements only yield statistical probabilities
due to the noise of the detector, it is instructive to discuss this idealized case 
first, since real measurements can never be more constraining.
We thus consider the case that the inspiral GW signal was 
measured with perfect accuracy and agrees exactly with a waveform model with 
parameters q,Λ̃. 
As always, we assume a fixed chirp mass, zero 
spins, and that S̅ has negligible impact on the waveform.

The exact measurement of q,Λ̃ corresponds to a single point in the 
(Λ̅, M̅) parameter space shown in fig:lambda_m_vs_M_m_all.
A given EOS is compatible with the measurement if the corresponding function 
Λ̅(M̅) passes through the point. Note that for some of the nuclear 
physics EOS, the Λ̅(M̅) curves shown in the upper panel of 
fig:lambda_m_vs_M_m_all cross. 
With a single measurement of a neutron star binary at such an intersection, the 
two EOS would be completely indistinguishable. 
On the other hand, the curves for our set of EOS cross at shallow angles,
such that the Λ̅(M̅) curves remain similar when moving away from
the intersection point. How large the differences of two EOS away from an intersection 
point can become is a highly relevant open question. 

The similarity between intersecting Λ̅(M̅) curves does not extend 
to maximum mass NS, however. EOS models with otherwise very similar curves can have very 
different maximum masses. Counterexamples visible in fig:lambda_m_vs_M_m_all 
show that the maximum NS mass cannot be determined from a single inspiral-phase GW signal 
alone, although M_1 constitutes a lower limit (assuming slow rotation).

In general, the inverse problem of deducing the constraints on the actual EOS function 
P(ρ) from Λ̅ 
measured at one or several M̅ is a very nontrivial problem, partly because
the full parameter space of P(ρ) has infinite dimension. However, it is not necessary 
to solve this problem for GW data analysis. Since the only information the GW signal
carries is encoded in Λ̅, M̅, one could restrict signal analysis 
to this information. Testing model assumptions for compatibility or deriving 
information not in the GW signal by adding model assumptions can both be accomplished
in a step separate from GW data analysis.


For an ideal measurement of Λ̅ at some mass ratio q, the universal relations
Λ_a(Λ_s, q) yield exact values of the individual Λ_1, Λ_2,
unless Λ̅ is inside one of the excluded regions discussed in sec:univ.
In our notation, the universal relation Λ_a(Λ_s, q) can be written 
as S̅_u(Λ̅, q), with a range Λ̅ < Λ̅_e(q) excluded
by construction.
The accuracy of Λ_a(Λ_s, q) needs to be specified as part of the model assumptions. 
Since there is only one true EOS and one universal relation, the error is systematic. 
From the viewpoint of measurements, the model assumption is defined by the allowed
range for Λ_a.
A reasonable range is given by the largest fit residuals of the EOS used for calibrating 
the relations, although it depends on the selection of EOS models, which is part of the assumption.
A comparison of the slope parameter S̅ between universal relations and various EOS
is shown in fig:Sl_vs_Lambda_m. As one can see, the EOS we consider fulfill
4 < -S̅ < 9. This range could serve as a simpler replacement for the use of
universal relations Λ_a(Λ_s, q) together with comparable accuracy assumptions.

Figure <ref> also leads us to the following  observations.
First, the 
more complicated universal relations used in <cit.>
do not reduce the residuals considerably compared to the more robust S̅=6 prescription
used in <cit.>. Further, the spread for the EOS we considered is somewhat larger
than for the smaller set used to calibrate the universal relations in <cit.>.
Finally, S̅_u(Λ̅,q)
diverges for Λ̅ outside the range allowed by the EOS models 
considered here, and the error is essentially unknown. Using the same error
as for Λ̅ values compatible with the EOS models is not a well-motivated
assumption, given that the divergence is caused by extrapolation of a particular 
form of the fitting function.






 §.§ Bayesian analysis


The study of gravitational-wave data relies on Bayesian statistical methods.
The main ingredients are 
[label=(*)]

  * a theoretical model of the signal 
expected for given source and extrinsic parameters;

  * a characterization of the detector noise that allows one to estimate 
how likely a given stretch of data contains a given signal; and 

  * the prior probability distribution for the source 
parameters, which encodes the state of knowledge before any measurement.


The choice of a prior is a very ambiguous ingredient, although
there are some well-motivated constraints to start with.
For example, gravitational masses of NS are positive in general relativity.
Similarly, negative tidal deformabilities are usually excluded (see 
<cit.>).
Another important decision is whether to assume that the EOS of NS is universal, 
i.e. that both NS follow the same Λ(M). 
Further, the use of universal relations is equivalent to a particular choice of prior.

There is however no observational or theoretical knowledge that could motivate
a choice for the exact shape of the tidal deformability prior distribution.
In this situation, it is common to express the lack of knowledge by 
assuming a uniform prior within the bounds. However, this depends on the choice of 
parametrization and is therefore still ambiguous. For a NS, both Λ and 
log(Λ) are natural choices. A further complication is
given by the choice of variables for the multidimensional parameter space 
of a binary system, in particular, in connection with universal relations. 
We do not try to argue which prior is “correct” but describe
the various priors used in previous studies and their implications in detail.

The studies <cit.> 
provide results on the tidal deformabilities in 
terms of the posterior distributions, which combine the prior with the 
information added by the observational data. The main results are 
one- or two-sided credible intervals on the tidal deformability or 
derived quantities.
In the following sections, we investigate how much of the posterior 
depends on the prior assumptions, and how much information is added by
the observational data. 

With decreasing deformability, its impact on the waveform diminishes.
For a given signal-to-noise ratio, we expect that there is some scale 
Λ̃_t below which the ability of the detector to 
distinguish tidal effects gradually vanishes.
One should therefore expect that any lower limit 
depends more strongly on the prior the smaller its value, and that upper limits
depend less on the prior than lower limits. 


With this motivation, we scrutinize the lower credible bounds
that are already given by the priors.
To this end, we introduce a thought experiment we dubbed the 
uninformative detector test.
Consider a detector that measures with perfect accuracy the NS masses and whether 
or not the effective tidal deformability Λ̃ is greater than some 
threshold Λ̃_t 
or not, but does not provide the value of Λ̃. 
This is an idealized model of a real detector,
for which the ability to distinguish increases gradually instead of abruptly.
The posterior obtained by a hypothetical measurement with the uninformative 
detector is given by the prior, after cutting it off above Λ̃_t  and
normalizing it to unity (assuming the true value is below Λ̃_t).

Of course one can compute perfectly valid two-sided credible intervals 
(Λ̃_t^0, Λ̃_t^1) for this 
posterior. However, the lower limit Λ̃_t^0 is rather useless 
because it is entirely
given by the prior and the threshold Λ̃_t. By construction, 
the detector cannot distinguish any deformabilities below Λ̃_t, 
including 0. A measurement with such a detector could never be considered as 
observational proof of finite tidal effects (unless the true value is above Λ̃_t).

The lower limit Λ̃_t^0 can serve as a reasonable estimate for the scale 
below which lower limits are likely not informed by the data, but the prior.
As further motivation consider the likelihood P(d|θ,I) in eq:Bayes_law.
For the uninformative detector, the likelihood is constant for 
Λ̃<Λ̃_t and 0 otherwise. 
Next, consider an actual detection with a different likelihood function that is 
featureless at scales Λ̃_t^0.
In order to obtain Λ̃_m^0 < Λ̃_t^0 for the 
actual detection,
the corresponding likelihood below Λ̃_t^0 would need to
increase on average, compared to the uninformative case, or decrease
above Λ̃_t^0. Hence low values become more likely by such 
a measurement, and a two-sided credible interval is probably not meaningful.

In order to provide Λ̃_t^0 for GW170817, we still need an estimate for 
Λ̃_t. 
Here we choose Λ̃_t ≈ 800, motivated by the posterior shown 
in Fig. 11 of <cit.>, which falls off strongly at larger values, 
indicating that tidal effects of this size can in fact be distinguished. 
Note the scale Λ̃_t could also
be lower, which would imply a lower Λ̃_t^0. Since we use 
Λ̃_t^0 as a conservative
sanity check for the interpretation of lower limits, this is not a problem. 
The case Λ̃_m^0 ≤Λ̃_t^0 can be regarded as a strong 
indication, but not conclusive evidence, that Λ̃_m^0 is mainly a 
consequence of the prior assumptions. In this case, providing credible intervals 
is insufficient and a detailed comparison of posterior and prior shape is necessary.







 §.§ Universal relations in Bayesian analysis


We now turn to discuss the use of universal relations in Bayesian analysis
of GW signals. To focus on the essential problems, we consider the simplified case
where chirp mass and mass ratio are known exactly, but not the tidal effects.


In <cit.>, the tidal effects are encoded in one 
parameter Λ_s, for which a flat prior distribution in the range (0,5000) is 
assumed, independent of chirp mass and mass ratio.
The prior for the quantity that is actually measured, Λ̅, is then derived
from the universal relations. The result for the given chirp mass and various mass 
ratios is shown in fig:Prior_univ_Lambda_m. As already mentioned, the 
universal relations have poles that cause gaps and peaks at very small values of
Λ̅. More importantly, the strict prescription of the slope causes
gaps extending to larger values of Λ̅ at certain mass ratios, 
because otherwise Λ_1 becomes negative. 
This makes the analysis unsuitable for providing lower limits on Λ̅ based
on measurement. 

Further, the upper limit for Λ̅ depends on the mass ratio. The prior 
density at fixed Λ̅ is therefore reduced with increasing mass ratio,
as shown in fig:Prior_univ_Lambda_m (the prior for Λ̃ shows
a similar behavior).
This is not a physical assumption, but an artifact of the functional form of the universal
relation in conjunction with the choice of a flat prior with constant cutoff for Λ_s. 
For given Λ̅ or Λ̃, the prior used in <cit.> 
favors unequal mass systems, independent of the chirp mass.
Therefore we expect this issue
to cause a bias away from the equal mass case.
We note that the prior for the masses itself already favors unequal mass systems,
again not due to physical considerations but due to fixed limits prescribed for the 
chirp mass.

We conclude that the use of universal relations in Bayesian analysis 
—in contrast to an ideal measurement discussed in sec:ideal—
can also affect the measured quantity Λ̅, via the prior.
However, the prior used in <cit.> does not exhibit this problem,
thanks to an additional cutoff in Λ̃ and a simpler choice of
universal relations. 
The prior for Λ̃ from <cit.> is flat, 
without gaps, and independent of the mass ratio.


The universal relations are used in <cit.> to derive posteriors
for Λ_1 and Λ_2. We note that, for given q, Λ_1 and Λ_2 are
uniquely determined by the parameters Λ̅ and S̅. While Λ̅
enters the GW waveform model, S̅ does not and can therefore not be measured 
at all. All information on S̅ is derived from the universal relations. We
thus treat it as a statistical variable derived from Λ̅ based on a 
model (the universal relations).

The resulting prior distribution for S̅ is shown in fig:Prior_univ_Sl. 
For any mass ratio, the distributions are strongly peaked and asymmetric. The 
shape of the prior is given entirely by the functional form of the universal relation.
It is not related to the distribution of S̅ resulting for the EOS 
models (which would not be meaningful anyway since the EOS selection is arbitrary).
For the universal relations in <cit.>, the prior for 
S̅ is even more restrictive, given by a delta function, which is shown 
as a vertical line in fig:Prior_univ_Sl.

We note that the prior resulting from the universal relations is not used
directly in <cit.>. Instead, the prior for Λ_a
is broadened in order to account for the error of the universal 
relations, 
which is estimated from the residuals of the
EOS set used for calibration. The broadening of Λ_a corresponds
mainly to a broadening of S̅, but to some extent also of Λ̅.
The sharp features visible in Figs. <ref> 
and <ref> are smoothed in the final prior used in 
<cit.>, and the gaps in the prior of Λ̅ 
gain some partial support. 
This might alleviate, but not solve, the 
issues discussed here. 
For the rest of this article, we therefore ignore
the effects of the additional broadening.

We also note that deriving the prior of S̅ from a flat
prior of Λ_a and the universal relations and then broadening it to account for 
errors of the latter is not the only possible approach. 
Further, <cit.> uses error estimates based
on treating an arbitrary selection of EOS models as if it were statistical sampling 
of the true EOS, arriving at a complicated model for the error that depends on 
Λ_s and q.
A cleaner and simpler approach would be to directly prescribe a prior in S̅,  
e.g., a Gaussian, with mean value and broadness defining the model assumptions.
That way, one retains control over the shape of the prior distribution and considerably
simplifies the interpretation of the assumptions.



The gaps in the priors of Λ̅ used in <cit.>
in conjunction with the slope S̅ 
restricted by the universal relations
have an important consequence for the prior of Λ_2, which is shown
in fig:Prior_univ_Lambda2.
As one can see, it has even larger gaps than the prior for Λ̅.
For mass ratio 0.7, the gap extends up to Λ_2 ≈ 100.
The method is therefore particularly unsuitable to measure lower limits on Λ_2. 














We now apply the uninformative detector test introduced in sec:bayesanity
to the prior used in <cit.>. 
As in <cit.>, we use a flat prior in Λ_s, 
cut off at Λ_s=5000. 
Following <cit.>, we use the universal relations 
Λ_a(Λ_s, q) to compute the Λ_i of 
the two NS, and further the radii R_i = M_i / C_u(Λ_i) (we do
not apply the additional broadening of the prior, though).
We exclude samples from the prior that violate Λ_2 ≥Λ_1 > 0.
The fifth percentiles of the hypothetical posteriors of the radii that follow 
from the uninformative detector test are shown in 
fig:univ_M_R for a range of mass ratios, as a function of M_i(q).
As one can see, the radii are still as large as 
≈8 in the mass range of interest. 


What is the reason for those surprisingly large bounds from our uninformative 
hypothetical measurement? 
One cause is the form of the universal relation C_u(Λ), 
which is a second order polynomial in log(Λ). Considering any flat 
distribution in Λ, it is not surprising that the 5th/95th 
percentiles cover only a small range of radii. 
To some extent, the lower limits are also caused by the functional form of 
the universal relation Λ_a(Λ_s, q), which completely
excludes tidal deformabilites below a certain limit (depending strongly on the mass ratio).
Translating this exclusion zone to the radius of the lighter NS via C_u(Λ),
we find that radii
up to 8 are already excluded by construction. This is shown in
fig:univ_M_R as well. It also shows that for extreme mass ratios, further regions in 
parameter space become unreachable, which is not relevant for the detected event,
however, given that the mass ratio is constrained.


The lower limits published in <cit.> for the measured
posteriors are shown as well in fig:univ_M_R. As one can see,
the limits are quite similar to those obtained with our hypothetical
detector. Comparing the radii derived from the universal relations is 
not ideal to assess how much the observational data contributed to the 
posterior. It would be easier to interpret results for the quantity that
directly affects the waveform, Λ̃. 
Although <cit.> do not quote limits for Λ̃,
they provide 90% two-sided credible intervals Λ_1.4=190^+390_-120
for the deformability of a NS with mass M=1.4 M_⊙, obtained by assuming
a linear expansion of Λ(M)M^5 around M=1.4 M_⊙.
For comparison, we computed the fifth percentiles of Λ̅, which increase 
from 44 at mass ratio q=1 to 58 at q=0.7, and Λ̅≈50
at M̅=1.4 M_⊙.
The limits from the prior are roughly comparable to the published lower credible 
bound of 70.
Our findings indicate that the lower radius limits, as obtained in 
<cit.> using the universal relations prior, might be largely 
based on the prior assumptions and less on the observational data.
We stress that our qualitative discussion is not sufficient to entirely 
dismiss the contribution of the data to the lower bounds. 
Nevertheless, the results denoted in <cit.> as measurements 
of lower limits for the radii should not be interpreted as observational evidence 
for finite size effects yet, even under the assumption that both stars were 
neutron stars following the same EOS.


A more basic test than our uninformative detector is to visually
compare the prior and posterior distributions.
Such a comparison is shown in the top-left panel of Fig. 3 of
<cit.>.
The prior for the radius of the lighter NS is much less restrictive
than the posterior, also towards small radii.
This might be interpreted as confirmation that the lower bound
is indeed based on data, contrary to our previous discussion. 
However, there is a subtle but profound
pitfall: the priors for the radii shown in <cit.> 
are not restricted to the mass range where the posterior has support. 
In particular, the full prior contains a range of chirp masses.
The mass enters directly in the inferred radius when
using the compactness-deformability universal relation.
In other words, any measurement of the mass already restricts the radius 
range, even without measuring the deformability.

Comparing the radius posterior to the full prior is therefore not a valid
test to assess if the measurement of the deformability is informative.
Showing only the full prior of the radii obfuscates the fact that
the priors restricted to the almost exactly 
measured chirp mass are much more restrictive for the inferred radius.
In order to judge if the lower credible bound on the radius 
corresponds to a measurement of finite tidal effects or merely
reflects a measurement of the mass range in conjunction with the 
model assumptions,
one should instead compare the radius posterior to the radius prior 
restricted to measured chirp mass and mass ratio range.
A very similar problem was already discussed in the context of spin precession 
in binary black hole mergers <cit.>.

We carried out the same analysis also for the universal relations
from <cit.>. Similar to the previous discussion,
we consider fixed masses and a flat prior
for the overall tidal deformability, this time in terms of the parameter 
Λ_0 used in <cit.>, with the same cutoff at
Λ_0=5000. Since <cit.> remove tidal deformabilities 
violating limits from causality considerations <cit.>,
we apply a cut Λ̃≤ 51, derived in <cit.>
for the chirp mass of GW170817 and assuming a maximum NS mass ≥ 2 M_⊙.

As before, we apply our uninformative detector test
to obtain a hypothetical posterior. We find that the 
fifth percentile for the radii are very close to the lower bounds 
inferred in <cit.> for the actual detection. 
This is shown in fig:De_univ_M_R, together with the 
limits obtained for causal neutron-star models.

Figure 3 in <cit.> shows the posterior of the effective 
tidal deformability. The distribution
drops sharply to 0 a bit below the lower credible bound. Note,
however, that low values of Λ̃ incompatible with causality 
have been removed from the posterior shown in <cit.>,
which likely causes the drop. The limit Λ̃≤ 51 from 
<cit.> is comparable to the size of the gap.
The original posterior and the prior 
are not shown in comparison, but <cit.> provide
limits both before and after removing causality-violating points from 
the distribution. The removal increases the lower bounds for Λ̃ 
by ≈40. When using a prior without
the causality cut, both the uninformative detector and the analysis of actual data 
in <cit.> result in lower limits reduced by a similar amount.

Our findings indicate that the lower credible limits inferred by <cit.> 
are almost entirely due to the prior assumptions, while the measurement did not add
significant support for finite tidal deformabilities. 
This calls into question the claim of “first evidence for finite size effects using 
gravitational-wave observations” <cit.>.
Although <cit.> compute Bayes factors around 10 in disfavor of 
Λ̃<100 compared to the full range,
an erratum <cit.> reports that the employed thermodynamic integration is 
unsuitable for computing Bayes factors. The impact on the Bayes factor discussed 
above is not provided, however.


Besides the problems mentioned before, we also find that the lower limits on radii 
or tidal deformabilites published in <cit.>
are well below the validity range of the universal relations used in their
deviation. Figures <ref> and <ref> show the 
limits in comparison to radii obtained for a number of EOS, including those used to calibrate the 
relations used in <cit.>. 
The lower bounds for the radii are the result of evaluating phenomenological expressions 
in a range explicitly forbidden by the models used to derive the very same 
expressions (see sec:univ).










 §.§ EOS-agnostic Bayesian analysis


In this section, we discuss the analysis of tidal effects in <cit.>,
with a focus on the interpretation of lower limits. 
As in previous  sections, we use the idealized case of known masses and zero spins to
highlight problems that are not obvious in the full statistical analysis.

In <cit.> it is not assumed that both stars obey the same EOS, or that they
are neutron stars, and no universal relations are employed.
In detail, <cit.> employs a flat prior
for the tidal deformabilities Λ_1, Λ_2, cut at a maximum value
0 ≤Λ_1,Λ_2<5000. This implies that the tidal deformabilities
are not only allowed to differ for equal mass systems, but also that they are 
completely uncorrelated, i.e. it is not assumed that stars with similar mass 
are more likely to have similar deformabilities. 


The assumptions of uncorrelated deformabilities have an important consequence
for the prior of Λ̃. Because of the rectangular uniform support in
(Λ_1, Λ_2) space, 
it linearly approaches 0 near Λ̃=0,
with a slope that depends on the mass ratio.
Figure <ref> shows the prior resulting for Λ̃
restricted to various mass ratios.
Small values of Λ̃ are already disfavored by the prior distribution. 
Obviously, this needs to be considered before interpreting lower credible bounds
as a proof of finite size effects.




In <cit.>, the raw posterior of Λ̃ is divided 
by the marginalized prior of Λ̃ before computing credible intervals 
on Λ̃. 
This approach effectively corresponds to using a prior in Λ̃
that is flat after marginalization over mass ratio.
The motivation given in <cit.> is to address the issue above, 
although the physical motivation for the modified prior is
not explained. Since the behavior of the 
original prior is a direct consequence of the assumption of uncorrelated deformabilities,
scaling it to obtain a flat prior in Λ̃ implies that the posterior
obtained for Λ̃ is based on contrary assumptions.
We note that the new prior does not correspond to the assumption of a common EOS, since 
the NSs can still have the same masses but different deformabilities.
However, in this work we are not arguing what prior should be used, but how much
the ones used determine the results. 

First, we investigate more carefully the technical implications of the rescaling.
One problem is that scaling in a mass ratio independent way is not the only choice to
obtain a flat marginalized prior. This already introduces an ambiguity that has 
not been discussed.
More importantly, the resulting prior restricted to mass ratios relevant for the 
detected event is not flat at all. The reason is that the prior for the mass ratio 
strongly favors unequal mass systems (mainly due to the limits imposed on the chirp
mass) for which the prior of Λ̃ has increasingly large slope near 
Λ̃=0. 
For moderate mass ratios, the reweighted prior used in  <cit.> 
disfavors low values of Λ̃, despite the reweighting procedure.
This is shown in the top panel of fig:Prior_flat12_Lambda_eff.


To asses the contribution of the prior to the lower limits
we now apply the uninformative detector test introduced in sec:real
to the reweighted prior used in <cit.>.
We obtain a lower fifth percentile of
Λ̃=128 for the posterior, for all mass ratios q ≳ 0.6.
This should be compared to the fifth percentile of the (reweighted) posterior inferred
from the actual data in <cit.> for the low-spin prior, 
which is Λ̃=110 (the 90% credible bound of the mass ratio 
is q > 0.73). Under the assumption of low spins, <cit.> 
interprets the credible lower bound as evidence for finite tidal effects. 
In contrast, our discussion strongly indicates that  
the lower bound should be regarded mainly as a consequence of prior assumptions
in conjunction with the observational constraints on the mass ratio.

In the remainder of this section, we compare the EOS agnostic approach to the
studies based on universal relations.
Surprisingly, the results in <cit.> can be reinterpreted
under the assumption that both objects are neutron stars obeying the same EOS,
in the sense that the tidal deformability is a unique but unknown function of mass.
For this, one only has to exclude systems with mass ratio exactly 1,
which makes no difference because this set has infinitesimal support in parameter space.
Further, one has to allow any functional form of the tidal deformability Λ(M), 
including arbitrary steep gradients. In our notation,
this implies an unbound range -∞≤S̅≤∞.
Note all EOS considered here fulfill S̅<0. 


Figure <ref> illustrates the the prior distribution of the slope 
S̅ for fixed chirp mass and various fixed mass ratios.
As expected, the support increases towards positively and negatively infinite 
slopes when approaching equal mass. We recall that eq:def_sl then implies
dΛ / dM →±∞.
At the same time, the prior probability for any finite
slope range approaches 0 in the equal mass limit. 
For unequal masses, the prior support for a given range can also increase with
increasing q. In any case, there is a strong dependence on mass ratio.

Since the signal is (almost) independent of the 
slope, the slope posterior is given by the prior. Compared to a slope prior 
bounded to some model-dependent range, the posterior for quantities derived 
from the slope, such as the individual tidal deformabilities, is broadened.
Compared to a more plausible prior S̅<0, on can also expect
a strong bias of Λ_1 and Λ_2 towards larger and smaller values, 
respectively. 

For comparison, we also plot the range of S̅ obtained from our standard set of EOS.
Comparing also to fig:Prior_univ_Sl, one can see that the EOS-agnostic study and
the studies based on universal relations occupy opposing ends of a spectrum of assumptions.
The EOS-agnostic study is overly inclusive regarding the slope, 
compared to the range of theoretical EOS models, while the 
universal relations strongly restrict the slope.
The most restricted slope is given by the fixed value used in <cit.>,
while <cit.> further broadens the raw prior distribution resulting
from their universal relations to account for a certain amount of systematic error.
Given that the universal relation assumptions in <cit.> strongly 
restrain S̅, it is hardly surprising that the credible regions reported occupy 
a smaller volume in parameter space.









 §.§ Parametrized EOS approach




A complementary approach <cit.> to EOS-agnostic and universal 
relation methods is based on parametrizing the EOS and including 
the corresponding parameters in the statistical analysis. This approach is also slightly
problematic, as we explain in the following.
Imagine an observation where chirp mass and mass ratio are measured exactly, but for
Λ̅ only a probability distribution is known. Further, consider
a hypothetical EOS parametrization where one parameter P_u uniquely determines Λ̅ 
for the measured masses,
while the remaining parameters P_i have no influence on Λ̅, but do have an impact on 
the individual NS properties. 
Since the parameters P_i have no influence on the measured quantity, their 
posterior distribution is given by the prior assumed for P_i. 
Therefore, the posteriors for individual NS properties depend on the priors
for P_i.
For the parametrization 
used in <cit.>, the parameters are not as cleanly separated, 
which makes it difficult to assess the importance of the above. 
Note that the implicitly chosen prior in the parameters P_i takes over the role of
the analytic prescriptions for S_m in the universal relation method when
it comes to inferring individual NS properties. 
Knowledge of the priors for a parametrization P_i, P_u instead of the ones 
reported in <cit.>
is required to interpret the significance of single-NS tidal deformabilities 
and radii inferred with the parametrized EOS method.

Since there is no natural choice for the exact parametrization and no 
other experimental data were used, the shape 
of the prior distribution should be regarded as an arbitrary choice.
It is therefore important to assess how much the prior contributes to any results
derived from the posterior. 

The prior and posterior
for the radii are compared in the top-right panel of Fig. 3 in <cit.>.
The credible bounds for the prior distribution are not provided.
Visually, it seems that the prior already has little support below the published 
lower limits of the posterior. In addition, the comparison is complicated by 
the fact that the full prior is shown, not restricted to the measured mass range
(compare the discussion in sec:real).
It is therefore impossible to judge how much of the lower limit is based on a 
measurement of the masses in conjunction with the prescribed prior for the EOS 
parameters, and how much on an actual imprint of tidal effects in the data.

We note that instead of parametrizing the EOS, one could directly parametrize the 
function Λ(M), for example, using a polynomial. One can then
perform Bayesian parameter estimation of the corresponding coefficients.
This approach was suggested in an earlier work <cit.>.
When considering a single ideal measurement, an obvious problem arises.
As discussed earlier, a single ideal detection provides one data point 
Λ̅, M̅.  This single point cannot determine more
than one coefficient, or combination of coefficients. 
For example, if q≈ 1, then Λ̅
determines Λ(M) at M=M̅. The slope of
Λ(M) at M=M̅ on the other hand depends on S̅, which cannot be measured.
Higher derivatives also have no impact on the measured signal. For unequal 
mass systems, the system is equally underdetermined and one can only measure
one combination of coefficients. 
Bayesian parameter estimation nevertheless yields posteriors
for all the coefficients, and some might appear to be constrained by measurement.
However, the signal still only carries one degree of freedom regarding the coefficients.
We recall that the shape of the posterior distribution in Λ̅,M̅ 
around the true value has nothing to do with the behavior of the Λ(M)
curve, but rather reflects the detector noise and the relative impact of mass ratio
and deformability. It cannot be used to constrain the slope of Λ(M), for example.

The approach is also instructive regarding the parametrized EOS method, which in essence
also contains a parametrization of Λ(M), but adds further degrees of freedom
that leave the curve Λ(M) invariant. We note, however, that the requirement 
of causality imposes constraints on the EOS. Assuming there is only 
one branch of solutions, and not a transition to quark stars, for example, 
this implies constraints on Λ(M). Results about Λ(M)
that have no impact on the information in the signal, such as S̅,
are not completely arbitrary. It is however very difficult to disentangle the meaningful
constraints from the results of Bayesian analysis employing the parametrized EOS.





 §.§ Systematic errors of EOS constraints



In the following, we discuss how systematic errors of the models
for gravitational signals affect constraints of the tidal deformability 
that would result from a perfect detection (in the sense of arbitrary large 
signal-to-noise ratio) of the inspiral phase.
Naturally, the impact of waveform systematics in statistical analysis of 
actual measurements with noisy data cannot be any smaller.

In the sensitive frequency range of current detectors, gravitational waveforms of 
coalescing NSs can only be computed using semianalytic approximations. 
Although some models are calibrated with numerical relativity simulations,
we note that such simulations cover only the upper end of the sensitive 
frequency range of current detectors (in contrast to the BBH case, where the 
merger typically occurs at lower frequencies).

Currently, the error of these models can only be estimated by comparing
waveforms based on different approximations. 
In Bayesian statistical data analysis using matched filtering,
the systematic error of the waveform model causes a corresponding 
bias of the posterior probability distributions.
It is worth pointing out that the bias of tidal effects
is not simply determined by the accuracy with which the tidal effects 
are modeled, but is also sensitive to the accuracy of the underlying
waveform model for zero deformability (binary black hole case).

A central measure in GW analysis is the mismatch between waveforms.
The likelihood P(d|θ,I) in Eq. (<ref>)
depends directly on the mismatch between waveform model and data.
It expresses how much waveforms differ from the perspective of a given
detector, weighting different frequency components by detector sensitivity:

    F(h_1, h_2) 
       = 1 - max_ϕ_c, t_c(h_1^ϕ_c, t_c, h_2)/√((h_1,h_1)(h_2,h_2))
 
where h_1^ϕ_c, t_c denotes waveform h_1 shifted in time and phase by 
t_c and ϕ_c, and

    (h_1, h_2) 
       = 4(∫_f_min^f_maxh̃_1(f) h̃^*_2(f)/S_n(f)df ),

where a tilde denotes Fourier transform and S_n is the detector 
spectral noise density. 
In the following we use the advanced LIGO zero-detuned high-power 
design sensitivity curve (provided by the LALSuite 
<cit.> software library, see also <cit.>) and apply a high-frequency
cutoff at f_max = 1 when computing waveform mismatches.

Although we discuss the case of arbitrary loud signals, the results in this section 
are still specific to LIGO-type detectors. They do not depend on the overall 
sensitivity though, only the relative sensitivity at different frequencies.
Further, the choice of cutoff parameter does have a varying impact on the
estimates we will present. Since the cutoff frequency
in <cit.> is chosen dynamically using heuristic
criteria based on masses and tidal parameters, the exact impact of the
waveform systematics might be somewhat larger or smaller than our estimates for
fixed cutoff. We will not address this additional complication here, but 
note that it deserves consideration before claiming observational evidence
for tidal effects.

Before discussing the impact of waveform systematics in general, we start with a 
more specific question. When computing a waveform with one waveform model
evaluated at some example parameters, and computing the best match with a different
waveform model, how strongly do the parameters differ?
We focus on the idealized case of fixed chirp mass (using the GW170817 value) and zero spin,
working in a parameter space M̅, Λ̅.
Figure <ref> shows the best matches with regards to a few reference 
parameters, comparing the IMRPhenomD_NRTidal <cit.>
and TaylorF2(6PN) <cit.> waveform models.
We note that the two do not share the same point-particle baseline and that 
the TaylorF2 shows relatively large differences to more recent 
models (see also <cit.>). How the true error 
of the most sophisticated models compares to the differences found for this 
example pair is an important open question, but outside the scope of this work.

As one can see, the differences for reference points with zero deformability 
are comparable to the lower credible bounds for the tidal deformability 
given in <cit.> (with narrow spin prior). 
In contrast, the Bayesian lower bounds obtained in 
<cit.> for different waveform models (including 
the two examples used here) agree quite well. This is not necessarily 
a contradiction. As we discussed before, one would expect similar lower bounds
even when low tidal deformabilities cannot be distinguished at all at the given
sensitivity. 
The comparison of statistical parameter 
estimation studies with different waveform models is an important cross-check,
but it should not be used exclusively to estimate the systematic error due to 
the waveform. On the contrary, if lower bounds agree better than the 
uncertainties obtained from the measures discussed here, it should be taken as 
an indication that the lower bounds are not informative.



The comparison of best matches for different waveform models gives a reasonable 
lower limit for the waveform-related systematic parameter uncertainty. We note, 
however, that the parameter difference at the best match could a priori be 0 despite 
a large residual mismatch at the best matching parameter. Further, true waveform 
most likely differs from both models. 

We therefore devised a complementary error
estimate that is not based on the best matching parameters, but on the magnitude
of mismatch.
For this, we consider two waveform models A,B and two points P_1, P_2 in 
parameter space.
The mismatch between waveform model A evaluated at the two points can be used
as a measure of how mismatch changes with distance in parameter space.
The mismatch between waveform models A and B evaluated at the same point P_1 
can be used to give us a scale for mismatches between waveform models.
Comparing the two, we define a distance in parameter space 
that is symmetric in the two waveform models,

    D(P_1, P_2)    = 
    min(F(a_1, a_2),
          F(b_1, b_2))/max(F(a_1, b_1), 
                                   F(a_2, b_2)),

where a_i and b_i are the waveforms obtained with models A, B at points P_i.

If two parameter sets differ by less than D<1, they cannot be distinguished reliably 
by measurement as long it is not known if A or B is the correct model.
In practice, we have to assume that the true systematic error is bounded
below D<E if D is computed for a pair of actual waveform models, 
where E is a guess of how much larger the true error is compared to the differences 
between available models.

The contours E=1,2 are shown in fig:mismatch, again 
for IMRPhenomD_NRTidal and TaylorF2(6PN).
It is worth noting that for reference points with zero deformability,
the differences to the best matching parameters are much smaller
than the extent of the D=1 contours. This might indicate that 
the mismatch changes only slowly for low tidal deformability. The slower
the change in mismatch, the larger the errors of the deformability.
The location of the best match is probably not a good error estimate
in this case. 
This concerns studies like <cit.>, which use
the best match location to estimate systematic errors, and which might 
similarly underestimate the true parameter error.


As a more conservative error estimate, we can use the maximum of the two 
heuristic measures.
From this, we conclude that for mass ratios in the range 0.7≤ q ≤ 1
and fixed chirp masses around the one for GW170817, 
we cannot distinguish a deformability Λ̅≈ 200 from the 
binary black hole case Λ̅=0, unless the waveform model employed 
is known to be more accurate than the difference between
the IMRPhenomD_NRTidal and TaylorF2 models, i.e. E<1.


For comparison, <cit.> claims evidence for finite-size effects
based on lower limits Λ̃ > 38,54,45 obtained
for the three different priors considered therein (those limits are the ones
obtained without applying a manual lower cutoff motivated by causality constraints 
in NS models).
In light of our waveform systematic discussion, this claim implies a great deal of 
trust in the TaylorF2 waveform model used in <cit.>.
Similar arguments apply to <cit.>, which does not 
provide a limit for Λ̃,
but lower limits for the individual radii. Using the same universal compactness relations
as <cit.>, we find
that Λ̅ = 200 corresponds to a NS radius of 10.4  km for 
the equal mass case, which is again larger than the lower limit R>9.1 
given in <cit.>.









 §.§ Improving GW analysis



In the following, we propose an alternative to the universal relation and EOS-agnostic
analysis methods that avoids the shortcomings discussed in Secs. <ref>–<ref>.

The first ingredient is to prescribe the priors for Λ̅ and S̅ directly.
In contrast to the universal-relation-based parametrization used in <cit.>,
one retains direct control over the Λ̅ prior, avoiding also the gaps discussed
in sec:real.
It also avoids correlation with the prior for the mass ratio. 
One can simply prescribe a flat prior for Λ̅. This avoids disfavoring low values 
of Λ̅ prior, contrary to the EOS-agnostic flat Λ_1,Λ_2 prior.
We stress that this flat prior is still an ad-hoc choice not based on any independent observations. 
One should always check how much results depend on such a prior.


When using Λ̅ and S̅ as independent variables, 
the individual deformabilities are given by

    Λ_1    = q^-S̅/2Λ̅_0,    Λ_2    = q^S̅/2Λ̅_0,

where

    Λ̅_0    = Λ̅/w_1(q) q^-S̅/2 + w_2(q) q^S̅/2.

The upper cutoff for the Λ̅ prior can be set to any reasonable cutoff for
the tidal deformability of a single NS. From  Eqs. (<ref>) and (<ref>) 
it is easy to show that 
Λ̅≤max(Λ_1, Λ_2) for all mass ratios.

The main advantage of specifying the prior directly for Λ̅ is that the 
prior for the second parameter (here S̅) does not influence the posterior of Λ̅.
The reason is that Λ̅ is the only relevant parameter for the GW signal,
unless S̅ reaches values that would cause significant corrections in the waveform 
model employed.
Although the corrections may never become significant unless the signal-to-noise ratio
is much better than GW170817, this needs to be validated in any careful analysis. 
We propose performing parameter estimation with a prior in S̅ that is much wider 
than any model assumptions one intends to make. 
Most likely, it will turn out that various model assumptions of interest 
constrain S̅ much more than the observational data (at least for any detection 
comparable to GW170817). 
This simplifies the parameter
estimation task considerably because model assumptions for S̅ do not need to be
included during the costly parameter estimation computations, but can be incorporated in a 
postprocessing step, as we detail later.


The second ingredient is to express EOS constraints inferred from the measurement
in terms of credible regions in the parameter space M_c, M̅, Λ̅.
Compatibility with a given EOS can be tested by comparing to the two-dimensional 
subspace allowed by the EOS. For a detection similar to GW170817, the chirp mass is
known very accurately and can be treated as fixed. In that case, the EOS surface
reduces to one-dimensional curves in (M̅, Λ̅) space, as in 
fig:lambda_m_vs_M_m_all.
This approach also allows comparing to the mixed binary hypothesis with given EOS in a unified way, without
resorting to separate parameter estimation studies with parameter spaces of different 
dimensionality, or to the EOS-agnostic method. 
It also allows one to judge the compatibility of arbitrary EOS models
without new parameter estimation computations.

One should avoid providing results from GW data primarily in terms of Bayesian 
credible regions on the NS mass-radius relations. 
Although it simplifies comparison to astronomical observations of isolated pulsars 
and allows one to use known mass-radius curves for EOS models, it also
requires additional model-dependent assumptions, e.g., universal relations.
As we have seen, it is nontrivial to assess how much the final result depends on
measurement and how much on model assumptions. In contrast, M_c, M̅, and 
Λ̅ are directly constrained by the GW signal and do not require model 
assumptions. Presenting results in this parameter space seems the cleanest line of
separation between the overlapping problems of GW data analysis on one side and 
NS structure models on the other.

When results on individual NS properties are presented, the unavoidable direct dependency 
on model-specific constraints of S̅ needs to be quantified. Otherwise, it is impossible
to judge how much results depend on assumptions and how much on measurement constraints.
We note that current models only yield approximate ranges for S̅, possibly as a function 
M_c, M̅, and Λ̅,
but not a well-motivated shape of the prior distribution. 
As we remarked earlier, S̅ is typically much more constrained
by assumption than by measurement. 
It is therefore not meaningful to repeat expensive 
parameter estimation computations with different priors for S̅. 
Instead, one can apply cuts to the posterior, limiting S̅ to the range allowed
by a given model assumption.

To get an estimate on the impact of model assumptions constraining S̅ to a range,
on can further marginalize the restricted posterior over S̅, and then derive two new 
posteriors for each of the individual deformabilities Λ_1, Λ_2 by setting 
S̅ to the extreme values allowed by the model constraints.
Repeating this for different models provides a computationally inexpensive 
way to separate model uncertainties and measurement error.

What type of credible intervals, e.g., two-sided symmetric, highest posterior density (HPD), 
or upper limits,
are meaningful to provide for the statistical error 
depends on the posterior shape. For the tidal deformability, there likely exists a 
reliable one-sided upper limit, since larger tidal effects are increasingly 
easier to detect. For a detection similar to GW170817, lower limits require more care,
as argued in previous sections. We note that using a HPD interval, as done in 
<cit.>, does not automatically 
yield a meaningful choice. For example, after adding an infinitesimal perturbation
to a flat posterior, the HPD choice can result in anything from one-sided lower to
one-sided upper limits.



In general, one should compare the posterior and prior in detail, also 
distinguishing the impact of mass constraints and deformability constraints.
This is simpler when using a flat prior in Λ̅ as proposed, 
since the latter is not correlated with the masses. 
It is then not necessary to restrict 
the prior to the measured mass range for the comparison, in contrast to the universal
relation method (see discussion in sec:real).

A simple diagnostic check is given by the uninformative detector test introduced 
in sec:real. The cutoff value 800, motivated by GW170817, should of 
course be replaced with the one-sided upper limit for the posterior at hand.
Lower bounds should be interpreted with care if the lower bound 
is not larger than the one given by the uninformative detector test. In such cases, 
lower bounds should not be regarded as observational constraints, at least not 
without a more careful analysis.
Even when passing this test, one should 
still perform Bayesian model selection as in <cit.>, since credible 
intervals do not imply statements on the likelihood of specific 
models (see Sec. <ref>).


The waveform systematics should not be assessed solely by comparing 
parameter estimation results using different waveforms.
It is possible that statistical bounds agree much better than the systematic 
uncertainty, if the bounds are determined mainly by the prior.
One should use mismatch computations as shown in sec:syserr as a guideline,
and not trust lower limits below values for which the waveform differences
are potentially larger.

We conclude with some remarks on incorporating EM counterparts. Improvements in 
numerical modeling of counterparts likely allows one to obtain robust 
predictions of EM counterparts.
However, there is no reason to assume that the parameter
combinations that are most constrained by EM counterparts
will bear any similarity to those combinations constrained by the GW signal.
In particular, it might not be possible to express the impact of the EOS on the 
kilonova signal solely in terms of tidal deformability. In <cit.>, 
a few numerical simulations were used to propose a loose correlation between effective tidal 
deformability and an upper limit for mass ejection. 
A later work <cit.> considered 
different EOS and found counterexamples, also demonstrating a large impact of
the mass ratio. For reliable results, it is likely necessary to consistently 
use a parameter space with enough parameters to uniquely characterize 
both EM and GW observable signals. 
One way to construct such a space is to use
a parametrized EOS both for GW data analysis and numerical simulations predicting
EM counterparts (the latter also requires one to parametrize the thermal part).
Since numerical simulations are computationally expensive,
analytic fits cannot be avoided completely. In order to prevent conceptual 
problems similar to those we discussed for the universal relations, 
they should however be based on data points in the parameter region relevant for a given event.








§ SUMMARY




In this work, we investigated the question whether results presented in 
<cit.> could be interpreted
as observational evidence for tidal effects during the merger that caused 
gravitational-wave event GW170817. 
Such a claim was made in <cit.>, and also in 
<cit.> assuming low NS spins.
Under the assumption that both objects are NSs, the authors of
Ref. <cit.> claim a measurement of lower bounds for the NS radii. 
Since those are derived by feeding the tidal deformabilities into a given monotonic 
analytic expression, this claim implies a measurement of lower limits for tidal effects.

All of those works use Bayesian statistics, but do not satisfactorily address the 
question of how much the results are already determined by the prior assumptions, and 
how much the observational data contribute.
We discussed the prior probability distributions in detail. 
Their features are based largely on arbitrary choices,
which is particularly problematic for lower limits on tidal effects, 
since small tidal effects have a weaker impact on the measured signal.

Since the detector cannot distinguish small tidal effects as well as large ones,
we discussed the idealized case of a detector that cannot distinguish
tidal effects at all if they are below published upper credible bounds.
Using such manifestly uninformative hypothetical measurement with the 
same prior assumptions, we still obtained lower limits comparable 
to the published ones.
Our findings indicate that the Bayesian lower credible bounds in
<cit.>
might largely reflect prior assumptions.
They should not be regarded as observational evidence, at least not without 
a more detailed analysis.


During our study, we also found several technical and conceptual problems. 
The two studies 
<cit.> are using model assumptions based on 
universal relations in a way that is not self-consistent. In short, those 
relations already imply much larger lower limits just based on the measured
mass range, and hence cannot be used to measure finite tidal effects.
While the use of universal relations in <cit.> has no consistent
motivation when applied to low tidal deformabilties for the mass range of GW170817, 
it still allows such values.
The universal relations used in <cit.> on the other hand 
completely exclude values below a limit depending on mass ratio. This limit
is not a physical assumption, but results from extrapolating a complicated 
fit function given by a fraction of polynominals.

Regarding the alternative approach of parametrized EOS discussed in 
<cit.>, we presented arguments suggesting that the impact
of the prior is not sufficiently studied to claim observational evidence.
However, our discussion is not conclusive since we did not quantify the impact.

Our reassessment of the EOS-agnostic study <cit.> also 
revealed some problems. On the conceptual level, we argued that the use of a 
rescaled posterior, which mimics a flat prior in the effective tidal 
deformability, in effect corresponds to a prior that not based on well-defined 
physical assumptions.
Moreover, it is ambiguous how the prior in the multidimensional parameter
space should be rescaled to achieve a given marginalized prior. 
On the technical level, we investigated 
the correlation with the mass ratio and found that the full rescaled prior
from <cit.> is not flat at all when restricted to the 
measured mass range, defeating the apparent purpose of the procedure.

For use in future studies, we sketched a GW data analysis strategy that sidesteps
some of the emerging problems. The main source of complications seems to be
the inclusion of quantities that cannot be measured, necessitating the use of 
model assumptions, and complicating the interpretation of results considerably.
The main idea of our proposal is to restrict the data analysis to quantities that 
actually affect the waveform and can hence be measured. Compatibility with model
assumptions can be tested in a second step, as well as derivation of quantities 
that require further model assumptions. 
Our strategy also simplifies quantifying
the impact of other model assumptions independently.

We demonstrated explicitly how to test compatibility with EOS models.
This also showed that, based on a single event with masses both compatible with 
NS, it is fundamentally impossible to distinguish a BNS merger 
with one EOS from a BH-NS merger with a different EOS, as long as tidal effects enter
the waveform model only via the effective deformability.
We also provided examples demonstrating that the maximum mass of NSs
cannot be inferred from the inspiral GW signal alone.

Another aspect particularly important for measuring finite tidal effects
is the impact of systematic waveform model uncertainty.
This was addressed in <cit.> by repeating the full 
statistical analysis with different waveform models. 
A similar approach
is to perform parameter estimation studies based on 
one waveform model on artificial data constructed using a different waveform
model. We argued that this 
is only a valid error estimate for credible bounds if they depend on the 
observational data, instead of the prior. We derived an independent estimate
by directly computing differences between waveforms, using the mismatch 
measure relevant for GW data analysis for LIGO-type detectors.
The results from this estimate indicate that also previous error estimates
based on best-matching waveforms, e.g., <cit.>, might be too optimistic.
We find that the uncertainties expected from differences of current waveform
models are larger than the lower limits inferred for GW170817. 

Conclusive 
observational evidence of finite tidal effect therefore requires  
better knowledge of the waveform accuracy, as well as another merger event
with either a larger signal-to-noise ratio or involving objects with larger 
deformabilties. For standard NS models, that means
less massive NS. However, the mass range of observed NS is rather limited
so it is not guaranteed that such an event will occur.







We thank Soumi De, Duncan A. Brown, and Katerina Chatziioannou for useful 
discussions, and Philippe Landry for a careful review of the initial draft of 
this paper. This work was supported by the Max Planck Society's Independent 
Research Group Programme. 
In this work, we use a neutron-star model database computed on
the Holodeck cluster at the Albert Einstein Insitute in Hannover. 



apsrev4-1-noeprint



