Constraining the quantum gravity polymer scale using LIGO data

We present the first empirical constraints on the polymer scale describing polymer quantized GWs propagating on a classical background. These constraints are determined from the polymer-induced deviation from the classically predicted propagation speed of GWs. We leverage posterior information on the propagation speed of GWs from two previously reported sources: (1) inter-detector arrival time delays for signals from the LIGO-Virgo Collaboration’s first gravitational-wave transient catalog, GWTC1, and (2) from arrival time delays between GW signal GW170817 and its associated gamma-ray burst GRB170817A. For pure-GW constraints, we find relatively uninformative combined constraints of ν=0.96+0.15−0.21×10−53kg1/2 and μ=0.94+0.75−0.20×10−48kg1/2⋅s at the 90% credible level for the two polymer quantization schemes, where ν and µ refer to polymer parameters associated to the polymer quantization schemes of propagating gravitational degrees of freedom. For constraints from GW170817/GRB170817A, we report much more stringent constraints of νlow=2.66+0.60−0.10×10−56 , νhigh=2.66+0.45−0.10×10−56 and μlow=2.84+0.64−0.11×10−52 , μhigh=2.76+0.46−0.11×10−52 for both representations of polymer quantization and two choices of spin prior indicated by the subscript. Additionally, we explore the effect of varying the lag between emission of the GW and EM signals in the multimessenger case.

We present the first empirical constraints on the polymer scale describing polymer quantized GWs propagating on a classical background.These constraints are determined from the polymer-induced deviation from the classically predicted propagation speed of GWs.We leverage posterior information on the propagation speed of GWs from two previously reported sources: 1) inter-detector arrival time delays for signals from the LIGO-Virgo Collaboration's first gravitational-wave transient catalog, GWTC1, and 2) from arrival time delays between GW signal GW170817 and its associated gamma-ray burst GRB170817A.For pure-GW constraints, we find relatively uninformative combined constraints of ν = 0.96 +0. 15  −0.21 × 10 −53 kg 1/2 and µ = 0.94 +0.75 −0.20 × 10 −48 kg 1/2 • s at the 90% credible level for the two polymer quantization schemes, where ν and µ refer to polymer parameters associated to the polymer quantization schemes of propagating gravitational degrees of freedom.For constraints from GW170817/GRB170817A, we report much more stringent constraints of ν low = 2.66 +0.60  −0.10 × 10 −56 , ν high = 2.66 +0.45 −0.10 × 10 −56 and µ low = 2.84 +0.64 −0.11 × 10 −52 , µ high = 2.76 +0. 46 −0.11 × 10 −52 for both representations of polymer quantization and two choices of spin prior indicated by the subscript.Additionally, we explore the effect of varying the lag between emission of the GW and EM signals in the multimessenger case.

I. INTRODUCTION
The growing roster of significant gravitational wave (GW) observations continues to provide invaluable insight into the nature of the cosmos [1][2][3].These signals, produced from the collisions of compact objects such as neutron stars and black holes, have profoundly expanded-and continue to expand-the catalog of astrophysical objects in our Universe as well as the properties that describe them.Now the exciting prospect of probing fundamental physics with these signals is upon us, with even more sensitive GW observatories on the horizon [1][2][3].Perhaps one the most enticing prospects of precision GW observations is experimental evidence of the quantum nature of spacetime.
On the other hand, since the formulation of classical general relativity, most of the theory community have a consensus that gravity, similar to other field of nature, is intrinsically quantum, and thus classical gravity is just a low energy limit of a full theory of quantum gravity.There have been several proposals for such a theory of quantum gravity, none of which are complete as of now.One of the candidates is loop quantum gravity (LQG) [4][5][6] which is a non-perturbative approach proposing quantum states of space are superpositions of gauge-invariant graphs whose edges have labels associated to the gauge group of the theory.The theory is written in terms of a certain connection called the Ashtekar-Barbero connection and the configuration variables are holonomies of this connection over paths in space.
The use of holonomies as configuration variables has inspired another model for quantization of both space(time) and matter fields with finite degrees of freedom called the polymer quantization [7][8][9].This model is closely related to employing Weyl algebra instead of the usual commutation (or Poisson in the classical case) algebra, and hence has at least two representations 1 .Mathematically speaking, this amounts to the use of infinitesimal generators (algebra members) associated to some of the canonical variables, and finite generators (group members) associated to the canonical conjugated variables.Those finite generator mimicking exponentials of infinitesimal generators look very similar to holonomies, and hence this method of quantization looks quite similar to LQG.The dynamics of the polymer approach leads to the quantization/discretization of the canonical variables that are not "exponentiated", since their conjugate variables are holonomized or "exponentiated" and thus only generate finite transformation of their conjugate non-holonomized counterparts.This method can be applied to both matter, e.g., GRBs [10] and spacetime or its perturbations themselves [11,12].
If gravity or spacetime as a whole is quantized, it implies that its propagating perturbations, i.e., GWs, would also be quantized.When modeling quantum extensions to GWs, there are two approaches: 1) seek out a full quantization of spacetime and derive the quantum-corrected gravitational waveforms as a consequence of the new theory, or 2) adopt a semiclassical approach, separating the GW from its background, and quantize either the background or the GW signal.In this study, we take the latter approach, quantizing the classical transverse-traceless tensor perturbations and assuming that the background spacetime evolves classically.In our previous work [11], we derived the equations of motion for polymer-quantized plane wave metric perturbations and provided numerical solutions, as well as approximate analytical solutions, to those equations of motion.Expanding on this work, we have also calculated the alterations to the response of a Michelson-Morley-like GW observatory under the influence of polymer effects [12,13].The polymer quantization model has also been applied to the propagation of gamma-ray bursts [10].
While there have been numerous proposals for tests of the quantum nature of gravity through observations of GWs [14][15][16][17][18], there has yet to be a test of the polymer quantization as a model for the propagation of GWs.In this work, we provide the first constraints on the polymer scale from GW observations from LIGO's first GW transient catalog (GWTC1) and from the joint detection of the first multimessenger GW event, GW170817 [19], and its coincident gamma-ray burst signal, GRB170817A [20].
The present paper is structured as follows.In section II, we review the procedure for polymer quantizing spacetime metric perturbations using the plane wave approximation, as described in [11].We also restate the results for geodesic deviation in orthogonal GW detectors such as LIGO or Virgo.In section III, we describe how the prediction that polymer effects will cause GWs to propagate slower than classically predicted linearized metric perturbations, can be translated into constraints on the scales introduced during polymer quantization.We outline the process for converting two independent constraints on deviations from classical propagation speed: i) Constraints on ∆v g arising from differences in arrival times and distances between detectors.
ii) Constraints on ∆v g from multimessenger events exhibiting time differences between the GW signal and its associated electromagnetic (EM) counterpart.
When applying these methods to events released in the first GW Transient Catalog (GWTC1 [21]), we find that much more informative constraints can be obtained when using multimessenger information on ∆v g , subject to the caveat that such constraints are highly sensitive to systematics.Constraints from inter-detector time delays, while much less informative by several orders of magnitude, may be improved with additional GW event data.Finally, in section IV, we discuss the implications of our results for both approaches and their potential for motivating future studies.

II. BRIEF REVIEW OF THE THEORY
We start from Einstein-Hilbert gravitational action where κ 2 ≡ 8πG/c 4 , with a general perturbed metric in which η µν is the Minkowski metric, and h µν denotes GWs as a small perturbation over η µν .We express the perturbation in transverse-traceless gauge as where h = h µ µ = η µν h µν .Given the two polarizations of the GWs and their properties in this gauge, hµν can be expanded explicitly as where h+ , h× are the aforementioned two polarizations of the GW, and (2.5) The effective equations of motion of the independent polarization modes of the waves then reduces to the familiar Klein-Gordon equation, where ȟλ (x) := hλ (x)/2κ and □ := η µν ∂ µ ∂ ν .The conjugate momentum πλ to ȟλ is derived as usual using the formula πλ = ∂L ȟ ∂∂t ȟλ from the Lagrangian density of the perturbations in terms of ȟλ given by which is written up to second order in linear perturbations.
In the previous expressions, we have The classical solutions of the equation (2.6) and their conjugates, in Fourier modes are where the wave vector k = (k 1 , k 2 , k 3 ) ∈ (2πZ/ℓ) 3 spans to a three-dimensional lattice L .The canonically conjugate variables h λ,k and The reality conditions on the fields indicates that not all the modes are independent.To have an independent expansion for each mode and write the Hamiltonian as a set of decoupled harmonic oscillators, we introduce new new variables A λ,k and E λ,k and split the lattice L into positive and negative sectors; for more details see [11,12].In terms of these new variables, the Hamiltonian of the perturbation field reads where k = |k|.
We now proceed to the polymer quantization of the Hamiltonians in (2.9) in order to extract out the effective terms in the classical limit.To do so, let us first provide the main ideas of the polymer quantization in a self consistent description.To begin with, recall that the fundamental observables required for the Dirac quantization scheme are A and E (see [11] for more details), each for a set (λ, k), and satisfying the Poisson bracket where δ λ,λ ′ and δ k,k are Kronecker delta functions and the other brackets are null.These Poisson relations are now used to construct the Weyl algebra whose generators, denoted as W (a, b), satisfy the algebra multiplication A better notion of these generators W (a, b) is obtained once we recall that in the standard representation, they can be written as the operator W (a, b) = e i(aA+bE) and their domain is the entire Hilbert space used in the standard representation.Actually, in the standard representation one can recover the fundamental operators A and E due to the weak continuity condition of the Stone-von Neumann theorem is satisfied.However, polymer quantization violates such a condition and therefore, it is not unitarily equivalent to the standard quantum mechanics [11,12].In this case, the Hilbert space is given by L 2 (R, dx Bohr ), where the configuration space, R, is the Bohr compactification of the real line and dx Bohr is the Haar measure.Depending on which variable we want to discretize, A or E, the Haar measure will have units of E or A, respectively.We call the first case "polymer E" and "polymer A" the second case.In the polymer E case, the representation is given by and for the polymer A the representation takes the form The main feature of these representations is of course the discrete spectra of the operators A and E which removes the possibility of having a representation for its canonical partner.That is to say, if A is discrete, then there is no polymer representation for E but for its "exponential" form W (0, b) and similarly for the case of discrete E in which the canonical variable is given in the "exponential" form W (a, 0).
Depending on the case we are considering, we now impose that every value on the spectra A j (or E j ) for the discrete operators can be written as 14) where µ and ν are the polymer scales of the system.These polymer scales are considered as the fundamental "lengths" for each of the cases.The parameters A (0) j and E (0) j respectively represent the center of the lattice for the polymer states.This can be confirmed by providing an example of polymer state written in the eigenbasis of the discrete operators, also known as almost periodic functions ) (2.17) At this point, we are not interested in the quantum analysis of the Hamiltonians in (2.9) but in their semiclassical versions.To do so, we apply the procedure given in [22] using path integral analysis.The idea is to obtain the effective action associated with the polymer Hamiltonians using the instanton methods developed for quantum chromodynamic models (for broader polymer examples using group averaging techniques see also [23]).The result yields a modification of the classical Hamiltonians in which the kinetic term is modified in the case of discrete A and in the case of discrete E the quadratic harmonic potential.We call these modified Hamiltonians effective Hamiltonians2 .
This results in two polymer effective (non-operator) Hamiltonians.For polymer E case Hamiltonian we obtain and for polymer A case the Hamiltonian becomes In the plane-wave regime, to the leading order in both polymer A and polymer E cases respectively, we obtain the following GW solutions to the equations of motion, and (2.21) with the group velocities (2.22) Here we introduced new polymer parameters ) where μ has the dimension of length, and ν is dimensionless in natural units.Depending on the quantization scheme, the velocity is sensitive to either a characteristic polymer length scale, μ, or momentum scale, ν.The model for the arms of the GW detectors is a system of two free-falling masses.The geodesic separation equation of these masses are sensitive to the metric perturbations A λ,k , i.e., the incident GWs which play the role of a source in this system.The perturbative solutions to the geodesic deviation of the two arms are [11,12] )

III. POLYMER CONSTRAINTS FROM GW ASTRONOMY
In this section, we outline the procedure for leveraging polymer scale-dependent departures from the classical GW propagation speed to place constraints on the polymer scale.The deviation can be inferred from the equations for the group velocity Eq. (2.22), Clearly, Eqs.(3.1) and (3.2) imply any measurements of ∆v g can be used to infer constraints on the polymer scales.Formally, the extracted probability distribution on the propagation speed, denoted as p(∆v g ), is related to the distribution of polymer scales which, we denoted as p {A,E} (U ), via a Jacobian transformation.By defining dimensionless parameters U A := hI ν ℏ and U E (k) := hI μ ℏ k in natural units, the relation between the two probabilities distributions, p {A,E} (U ) and p(∆v g ) becomes Of course, this relation is valid with the corresponding measures (which we omitted for simplicity) and normalizations of these distributions inside the integrals.Invoking Eqs.(3.1) and (3.2) to evaluate Eq. (3.3), we arrive at a set of simple relations between p {A,E} and p(∆v g ): for the probability distribution for two cases of polymerization.
In the following sections, we use constraints on ∆v g from two independent approaches: The first relies on inter-detector arrival time differences for signals detected in multiple GW observatories, while the second compares the arrival time difference between multimessenger GW signals and their electromagntic counterpart.We apply this procedure to event GW170817 and its associated GRB, GRB170817A, as this is so far the only existing confident multimessenger detection.
A. Constraints from inter-network arrival time delays

Methods
The following section closely follows section II of [24].In Ref. [24], the propagation speed of the GWs is treated as a free model parameter which deviates from its typical treatment where it is fixed to be the speed of light.The canonical procedure for extracting such parameter information from GWs relies on techniques aimed at sampling the Bayesian posterior probability density.
We use statistical methods to obtain a probability distribution for the inter-detector time delays as a function of v g .This distribution enables us to determine the lower and upper bounds on the speed of gravitational waves.Assuming we have a network of m gravitational wave detectors, each separated by a light travel time ∆t c ij (time delay for light between detectors i and j), we use the relation∆t g ij = c/v g ∆t c ij to map light travel time to GW time delay.By considering uniformly distributed sources in the sky and the antenna patterns of the detectors, we can define a distribution of light time delays between every two detectors p(∆t c ij ).We use this distribution to define the likelihood p(∆t g ij |v g ) [24,25].The posterior distribution for v g for one event and between only two detectors follows from p(v g |∆t g ij ) = p(∆t g ij |v g )p(v g )/p(∆t g ij ), by the Bayes' theorem.Here, p(∆t g ij ) and p(v g ) are the normalization factor and prior knowledge about the distribution in v g , respectively.We assume that v g follows a uniform distribution.Assuming the data measured at detectors labeled by index i is composed of a signal and noise, where n i is the noise.The probability distribution on this parameter for a single event in a network of detectors can be computed by Bayes theorem: , known as the evidence, is a normalization factor useful for model selection that is largely irrelevant for our analysis, and so we do not explicitly compute this.Finally, the likelihood can be written in the frequency domain as assuming the noise in each detector, n i (t), is stationary and Gaussian-distributed, and adopt dt as our Fourier convention.The remaining components h i (f |v g ) and S i (f ) are the frequency-domain waveform and power spectral density (PSD) of the noise respectively.
Lastly, to properly leverage the plethora of GW event data available, Eq. (3.6) can be applied in iteration; for n independent GW events each with data e α , α = 1, . . ., n, the joint posterior on v g is proportional to the product of the posteriors from individual events: assuming a flat prior on v g .Finally, samples drawn from p(v g ) can be trivially converted to samples for the posterior distribution on ∆v g = v g − 1 leaving us with p(∆v g |e 1 , e 2 . . ., e n ), which can then be converted to constraints on the polymer scale following the procedure described in section III.

Data provenance and modeling
All posterior samples used in this work were generously provided by the authors of [24] in which they used Markov Chain Monte Carlo (MCMC) with Metropolis-Hastings algorithm to sample from multi-dimensional posterior distributions.All binary black hole (BBH) events were modeled using the IMRPhenomPv2 waveform [26], which models the inspiral, merger, and ringdown phases of the GW, and includes the effects of precession on the GW strain.The singular binary neutron star (BNS) event included, GW170817, was modeled using the TaylorF2 waveform, a post-Newtonian inspiral model that includes tidal distortions of the neutron matter.
To apply our model to a binary system we make two simplifying assumptions, which are forced on us from numerical perspectives while on the other hand captures the order of corrections induced by polymerization.We assume that the source is classical and that it produces GWs with initial amplitude and template for the evolution of the frequency of the inspiral phase and through propagation, the GW waveform receives corrections effectively for each polymer quantization scheme.
Using posterior samples from [24] and extracting posterior information on the polymer scale by Jacobian transformation (3.3), we technically assume the true effective waveform (i.e., including polymer corrections) is given by and that means by analyzing data with classical waveforms such as IMRPhenomPv2 and TaylorF2, we will miss some of the polymer corrections which may induce biases in the recovered posteriors.However, we sidestep these biases since the estimate of ∆v g depends solely on the time at which the GW signal amplitude peaks, a model-independent quantity.In a future study we hope to introduce our waveform corrections to LALSimulation and extract the posterior information for polymer parameters directly with polymer-corrected waveforms but this entails a significant project that we leave as a follow-up study.As the final remark, we should note that waveforms (3.9) are consistent with waveforms employed in templates IMRPhenomPv2 and TaylorF2 while only the speed of GWs is now modified and receives corrections (3.1) and (3.2).

Results
Fig. 1 displays the probability distribution function (PDF) for ∆v g in GWTC1 [21].The events with little or no support for negative ∆v g were excluded from the study.They are composed of 8 BBH events as well as GW170817, the lone BNS event.We combine information constraining ∆v g from these 9 events and apply the methods outlined in section III to resample and interpolate the distribution functions once more.Sky localization and high signal-to-noise ratio (SNR) of the events help us better constrain v g , therefore, some events have more sharp peaks.The two highest peaks in Fig. 1 correspond to events GW170817 and GW170814, which have high network SNR and their sky location is well constrained.On the other hand, some events like GW170823 and GW170104 have poorly constrained sky localizations and low SNR, thus their posterior plots of v g seem relatively flat.The resulting combined as well as per-event constraints on the polymer scale for each of the two polymer quantization choices are depicted in Figs.2a and 2b.A hard cut prior is placed on all ∆v g samples less than 0, as these are excluded by our model.To better present combined results of the events and their probability density distribution on polymer scale μ and ν, compute two combined results, in "CombinedBBH" only data from BBHs are analyzed, while in "Combined" case, data from GW170817, the multi-messenger BNS event is also included in the analysis.Fig. 2a shows that including or removing BNS event data from combined posterior analysis does not affect the results too much, only displacing the maximum of the probability peak by small value, while on the other hand, as presented in 2b, BNS data changes the combined analysis considerably.Unsmooth behavior of the combined case for μ is the result of its dependence on the frequency of each signal, we should note that the frequency of BNS event GW170817 is one order larger than the other events.The maximum posterior values of ∆v g and polymer parameters µ and ν are presented in Table I.To find max posterior values for the polymer scales, first, we compute the posterior on ν and μ and then find the maximum of the distribution.In the process, we do not use any of the positive values for ∆v g in our analysis, which is why the y-axis labels on the polymer scale PDFs now read p(μ|∆v g < 0) as we are effectively computing the conditional PDF.To justify this, we also computed the probability that ∆v g is negative for the combined PDFs on ∆v g , but we have not added them here, which shows more that 50% of the events.The required length scale ℓ for the binary system is set to 10 10 m, larger wavelengths are ignored and could be absorbed in the homogeneous background, because we assume our system is localized.
We should note here that we did not find any resources that had values for the strain and frequency at the peak for all the BBH events.We instead tried to simulate the time domain signals for each of the events using the maximum posterior values for each of the model parameters.Then, simulated each of the interferometer detector responses which accounts for the antenna function and approximates the noise characteristics using the published power spectral density of the noise for each of the event/detector pairs and found the maximum strain value in each of the detector responses.An example of the generated waveform is plotted in Fig. 3 for GW150914 event.For the frequency at maximum strain value, we just found the peak just before the merger and did a really rough frequency approximation from the time difference of the two peaks.Our maximum strain value has the same order of magnitude as the few reported values on the available factsheets of the events, but it differs by about 0.8 × 10 −22 .This is within our acceptable range of tolerance, because by assuming waveforms (3.9), we have already accounted for this level of uncertainty.After all, in the combined cases and corresponding values for polymer scales, the effect of these small tolerances will get even smaller.

B. Multi-messenger constraints
Multi-messenger astronomy has developed rapidly over the past years.The channel type of astronomical messengers now includes electromagnetic radiation, gravitational waves, neutrinos and cosmic rays.One of the main multi-messenger sources are binary pairs (BHs and NSs) [27], since their first detections in 2015 by LIGO and VIRGO [28], several techniques of astronomical Figure 1: Posterior density functions on parameter v g estimated for events in the first and second observing run of Advanced LIGO (O1 and O2), from which the events with little or no support for negative ∆v g are removed from the study.

Table I:
Locations of maximum a posteriori values of ∆v g , ν and μ for all the events, and the corresponding calculated polymer parameters in their reduced form ν and µ. "Combined" and "Combined (BBH)" refer to combined events data with/out the BNS event GW170817.Uncertainties listed are calculated to the 90% credible level.To have a better upper bound estimates for the polymer parameters, we use the frequency and strain of the peak of inspiral phase with tolerance about 0.8 × 10 −22 , where we assumed ℓ = 10 10 m for the length scale of the system.p E (μ|∆v (b) Probability density on polymer scale μ given ∆v g < 0 from gravitational-wave detections from LIGO's first and second observing runs.The combined posterior for μ has a maximum a posteriori value of μMP = 0.27 +0.05 −0.17 × 10 −12 kg • m 2 .

Figure 2:
Posterior probability density functions on polymer parameters extracted from conditional PDF of ∆v g .All events used for analysis are from the LVC's first gravitational-wave transient catalog paper (GWTC1).observations have been emerged.Observation of the first multi-messenger transient GW170817 [19,29] has raised interests to study the details of the physical processes in their sources from different perspectives.The gamma ray burst GRB 170817A was detected by the Fermi Gamma-ray Space Telescope and INTEGRAL 1.7 seconds after the gravitational wave signal GW170817, which was detected by the LIGO/Virgo collaboration in 2017.These signals were produced by the neutron star collision in the galaxy NGC 4993.In the event of an electromagnetic counterpart coincident with a GW detection, constraints on ∆v g can be placed based on the difference in arrival times between the coincident gravitational and electromagnetic signals as well as an estimate of the distance to the source.

Methods
Following the procedure of [20], deviations from the classically predicted group velocity of GWs can be derived from measurements of the time delay between coincident GW and electromagnetic signals, where ∆v = v GW − v EM , ∆t is the time delay between the two signals, and d L is the luminosity distance to the source.The time delay that appears in Eq. 3.10 is assumed to be caused purely from polymer effects.However, the observed time delay will in reality be a sum of the time delay due to polymer effects and any difference in the emission times of the GW and EM signals ∆t obs ≡ ∆t poly + ∆t lag .While Ref.
[20] predicts a 10 s lag time, others propose significantly longer lags [30,31] up to ∼ 1000s.Initially, we take ∆t lag to be perfectly known, but later on we will explore the dependence of the polymer scale measurements on the choice of ∆t lag .We take both d L and ∆t obs to be random variates, where p(d L ) is approximated from the publicly available posterior samples produced from LIGO parameter estimation analysis.The distribution p(∆t) is instead assumed to be normally distributed with expectation value E[∆t] = 1.74 s and standard deviation σ = 0.05 s, again in accordance with [20].For their lower bound estimate, Ref. assumes a 10 s lag between the emission of the source's gravitational radiation and its associated GRB.We adopt an even more conservative lag of 3.48 s, which is equivalent to simultaneous signal emission with polymer effects inducing a 1.74 s lag in the GW arrival time over the distance traveled by both signals.Compact binary coalescences are expected to be strong GW radiators associated with a delayed emission of short gamma-ray up to a few seconds compared to the GW emission, given that the central engine is expected to form within a few seconds during the inspiral phase [32,33].Therefore, an observer in direction of the outflow is expected to observe the GW/GRB signal with a delay up to few seconds for the electromagnetic counterpart emission.Other models suggest significantly longer lags [31,34], but we take this a conservative estimate.To compute the posterior distribution on ∆v g , we first build the probability distribution p d L (d L ) for d L from a public library.We then construct a normal probability distribution p(∆t) for time delay.The posterior distribution for ∆v g is proportional to the product of marginalized posteriors p d L (d L ) and p(∆t), as mentioned earlier, The constructed probability distribution on ∆t is related to the distribution of p(∆v g d L ), via a Jacobian transformation as (3.3) using relation (3.10).The resulting distribution on ∆v can then be computed by integrating over the joint probability density function, where v EM has been set to 1. Then the polymer scale distribution is then calculated from p(∆v g ) according to section III.

Results
The measured time delay between the arrival times of gravitational-wave detection GW170817 and the coincident gamma-ray detection GRB 170817A [20] constrains the deviations of the GW propagation speed from that of the electromagnetic radiation as, assuming a conservative estimate of luminosity distance to the source binary of 26 Mpc, the lower bound of the 90% credible interval.The upper limit is unphysical when interpreting ∆v as a purely polymer-induced effect, so we restrict our estimate to the lower limit.This is equivalent to enforcing that the polymer scale must be real-valued.By inverting Eqs.(3.1) and (3.2), constraint (3.13) can be translated to the following upper bounds on the polymer corrections.
What we have for polymer E is a frequency dependent correction, in line with the argument we presented and lead to equation (3.9), we replace its dependent with the max value, for the sake of simplicity and postpone its dependent to a future work.High spin prior Low spin prior Rewritting these expressions in terms of the μ and ν we obtain In addition to this point-statistic bound, we also estimate the PDF on the luminosity distance to the source of GW170817 from posterior samples provided in the LVC's public data release 3 .Approximating the PDF on the time delay as a Gaussian, we use Eq.(3.12) to compute the PDF on ∆v g , presented in Fig. 4. Finally, the resulting PDF on the polymer parameters under the two polymer quantization schemes are computed following the methods described in sections III and III B 1, depicted in Figs. 5. We calculate the equivalent 90% credible regions for the PDFs on the effective polymer parameters ν and μ and also bare parameters ν and µ report them in Table II.It should be noted that for this event, the choice of spin prior is particularly important.Higher spin values allow for the neutron star to sustain a higher mass, a parameter which is degenerate with luminosity distance at the level of the gravitational waveform.However, highly spinning neutron stars are thought to be rarer than ones with more moderate spins due to a loss of rotational energy through the powering of magnetically-driven plasma winds [35][36][37].Thus, Fig. 5 includes two posterior density functions, corresponding to two choices of spin prior: one which restricts the spin parameters to low values, and one which assumes all values of the spin parameters are equally likely.
The dependence of polymer constraints on the assumption of time delay between GW signal emission and GRB emission is displayed in Figs.6a and  Additionally, all variances for the distributions on {μ, ν} {low,high} ≤ 2.0 × 10 −3 .These plots show that the estimated constraints does show strong dependence on the time delay assumption we made during our analysis.

IV. DISCUSSION AND CONCLUSIONS
As presented in [11], polymer quantization affects the wave-form and the propagation speed of gravitational waves and predicts departures from classical GR.To place constraint on the polymer scale, we use two procedures for leveraging polymer scale-dependent departures.In the first approach, we employed inter-detector arrival time differences for detected GW signals in observatories to find probability distribution function on deviation of the propagation speed of polymer GWs compared to their GR-predicted classical propagation speed, ∆v g .After simulating and extracting the strain and frequency of the signals at the peak, we translate the constraints from ∆v g to polymer parameters ν and µ.The details of posterior data for each signal can be found in Table I.In the "Combined(BBH)" case, which we combined only the binary black hole events, we obtained constraints on polymer parameters, ν = 0.96 +0. 15  −0.21 × 10 −53 kg 1/2 and µ = 0.94 +0.75 −0.20 × 10 −48 kg 1/2 • s.After including the data for the single BNS event GW170817, the "Combined" value for ν and µ changed to 0.90 +0.16  −0.20 ×10 −53 kg 1/2 and 0.25 +0.05 −0.15 ×10 −48 kg 1/2 •s respectively.The listed uncertainties are calculated to the 90% credible level.In the second approach, we tried to find the constraints on the polymer parameters from a different method, by comparing the arrival time difference between multimessenger GW signal GW170817 with its EM counterpart GRB170817A.By assuming a conservative lag of 3.48 s, we extracted polymer scale constraints from max(∆v g ) for two spin priors, low and high.For the low spin case we find constrains ν = 2.66 +0.60 −0.10 × 10 −56 kg 1/2 and µ = 2.84 +0.64 −0.11 × 10 −52 kg 1/2 • s, respectively.For the high spin case the constraints turn out to be ν = 2.66 +0. 45 −0.10 × 10 −56 kg 1/2 and µ = 2.76 +0.46 −0.11 × 10 −52 kg 1/2 • s.In our previous work [12], from a completely different approach, we found bounds 10 −52 < ν < 10 −58 and 10 −44 < µ < 10 −50 for a given length scale ℓ for the binary system.Constraints on ν and µ extracted in the first procedure are within the detection range of LISA, while analysis in the second procedure shows that only quantum effects in the A scheme falls in the LISA range.One might tend to conclude that constrain on µ obtained from the second approach, makes predictions of the quantization scheme E, undetectable in the future observation of LISA, but one should note that these effects are frequency-dependent, and thus by going to higher frequencies, their chance of detection increases.On the other hand, constraints reported in [12], obtained for different length scale ℓ of the system, that means, we assumed different range of frequency for the GWs, by changing its value, different constraint can be extracted.It is not an odd feature, since, extracting any bounds for the polymer scale, is closely connected with the characteristic properties of the underling system, we will elaborate on this point when we compare our results with previously reported constraints.Although the main motivation the polymer quantization comes from quantizing the gravitational degrees of freedom, but most of the previously reported constrains on the polymer scale obtained by considering matter fields in different setups [38][39][40][41][42][43].In almost all of them, the reported bound on the polymer scale was sensitive to the characteristic properties of the setup, for example in [38,39], by changing the number of particles and the characteristic length of the one-dimensional oscillator, a different bound on the polymer scale can be obtained, or in [40,41], different value for the number of particles, size of the system or barrier width, would result in different bound on the polymer parameter.Even in [10], which employs the same procedure for the mode decomposition of the electromagnetic fields, final bounds on the polymer scale depends on the selected value for the size of the decomposition box and the amplitude of the observed GRB, which shows a similar role to the parameter ℓ in our setup.
When constraints on ∆v g are estimated from inter-detector time delays for pure GW signals, the resulting distributions on the polymer scales are relatively uninformative.The uncertainty on the distributions in Figs.2a and 2b are significantly larger than even theoretical constraints, with GW170817 providing the best single-event constraints due to its precise sky location measurement.There is, however, modest improvement when information from multiple GW events is combined.While our analysis only includes events from the LVC's first GW Transient Catalog (GWTC1) [21], the second [44,45] and third [46] catalogs add an additional 79 high-significance candidates to the list of GW detections.While additional constraints on the polymer scale from O(100) GW events may not be sufficient to make robust claims about the existence of polymer effects, next-generation GW detectors such as Cosmic Explorer (CE), and Einstein Telescope (ET) are set to provide O(10 6 ) observations by 2050 [47].Even in the likely case where polymer constraints do not improve with the current typical signal-to-noise ratio (SNR) of GW events (regardless of how many observations are made at that SNR), ET and CE are expected to make O(10 3 ) observations of BBH signals with SNR ≥ 100.
GW events accompanied by an associated electromagnetic signal offer an opportunity to constrain deviations from the classically predicted GW propagation speed from time delays over astrophysical scales.This translates into many orders-of-magnitude tighter constraints on the polymer scale as shown in Fig. 5.However, this approach is susceptible to systematics.Estimates of the time delay between the emission of GWs and GRBs in BNS systems varies widely which has a non-negligible impact on the resulting PDF on the polymer scales.We explore this dependency in Figs.6a and 6b.We find that the uncertainty is largely independent of the lag.While the mean of the PDF varies with the the emission time difference between the two signals, the means are within one order of magnitude of each other despite the emission delay ranging from 0 s to 100 s.This reflects the fact that when the propagation speed deviation is estimated as ∆v g ≈ ∆t/d L , the polymer scales go as U {A,E} ∼ √ ∆t.With better BNS merger modeling in both the GW and EM sectors, this systematic dependence can be ameliorated yielding more trustworthy polymer constraints.
Furthermore, we note that Eq. (3.2) has a spectral dependence.Since GWs have a non-trivial frequency evolution (their amplitude and frequency change over the inspiral, merger and ringdown phases), we should consider ∆v E (k) as a function and µ as a free parameter of the model, and infer its value from the posterior analysis.However, since the posterior for ∆v g was not binned in k-space (or equivalently, frequency space), we treated ∆v E as a free model parameter and mapped the posterior distribution from this parameter to the polymer parameter µ.We also had to choose a constant value for the frequency and set k = k value at the peak , aiming for the most optimistic constraint for µ.Frequency-dependent constraints are something we hope to explore in a future study, which would require having v g binned in the frequency space.
Our model shows that if spacetime is quantum with a minimal length scale, then this should result in the modification to the waveform of the gravitational waves, including their amplitude dispersion relation, and in particular will lead to the dependence of the speed of propagation of gravitational waves on their frequency.If such dependence is actually experimentally established, then our model shows how to get an indirect bound on such minimum scale from the aforementioned dependence of the propagation speed on the frequency of the waves.This, together with our precise prediction to modification of the waveform can lead to two outcomes: either this precise waveform will match near-future precision experiments and the results match our predictions, in which case it would be an strong indication of the quantum nature of spacetime, or in case of disagreement with experiment, this specific polymer model will be refuted.In our opinion, either case would be fruitful results.
Furthermore, although our results show that polymer effects will modify the propagation of GWs and an upper bound for the polymer parameters can be found using our suggested approach, nevertheless, in order to concretely obtain an indication of quantum gravity polymer effects, i.e., a smoking gun result, we need to also find lower bounds to the polymer parameters.A necessary (but probably not sufficient) improvement in this direction is to explore the aforementioned frequency dependency of the polymer parameters together with other wave-like effects, and particularly extending our model to cases where the background, as well as perturbations, is also quantum.
In future investigations we hope to perform a full forecasting study to quantify the level at which additional GW events of a certain SNR improve polymer constraints estimated from inter-detector time delays.Additionally, with the first-order analytic approximations to the full polymer-corrected gravitational waveforms (Eqs.(2.21), (2.20)), it is now feasible to directly constrain the polymer scale by performing Bayesian parameter estimation with waveforms that include polymer effects-a subject of another future study.
Probability density on polymer scale ν given ∆v g < 0 from gravitational-wave detections from LIGO's first and second observing runs.The combined posterior for ν has a maximum a posteriori value of νMP = 0.99 +0.18 −0.22 × 10 −17 kg • m 2 • s −1 .

Figure 3 :
Figure 3: Example of a waveform generated to find stain and frequency at the peak.Vertical lines show two consecutive peaks, which are used to extract frequency of GW at the peak.

Figure 4 :
Figure 4: Constraints on departure from classically predicted propagation speed of GWs, calculated based on estimates of luminosity distance to the source and time delay between GW170817 and GRB170817A.The maximum a posteriori values for ∆v g under the two spin priors are ∆v g,low /v EM = 3.99 +1.99 −0.31 × 10 −16 and ∆v g,high /v EM = 3.94 +1.44 −0.31 × 10 −16 .

Figure 6 :
Figure 6: Dependence of constraints of polymer scales ν and μ on the time delay assumption for the GW and EM signals.
.6) The prior distribution p(v g ) encodes any prior knowledge about what values the parameters can take on before a measurement is made.The denominator, p(d 1 , d 2 , . . ., d m ) = p(v g )p(d 1 , d 2 , . . ., d m |v g )dv g

Table II :
6b for both polymer quantization schemes and both choices for the spin prior on d L .The means of the distributions 90% credible intervals for PDFs on {μ, ν} {low,high} shown in Fig. 5. on ν are all within the range −19.50 ≤ log 10 νlow ≤ −18.62 and −19.51 ≤ log 10 νhigh ≤ −18.63 for low and high values for spin prior.While the means on the distribution on μ are all within the range −15.47 ≤ log 10 μlow ≤ −14.59 and −15.50 ≤ log 10 μhigh ≤ −14.62 for low and high values.