Multimessenger parameter estimation of GW170817: from jet structure to the Hubble constant
DDraft version September 10, 2020
Typeset using L A TEX default style in AASTeX63
Multi-messenger parameter estimation of GW170817: from jet structure to the Hubble constant
Hao Wang and Dimitrios Giannios Department of Physics and Astronomy, Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47907, USA (Received xxx; Revised yyy; Accepted zzz)
Submitted to ApJABSTRACTThe electromagnetic radiation that followed the neutron-star merger event GW170817 revealed thatgamma-ray burst afterglows from jets misaligned with our line of sight exhibit a light curve with slowlyrising flux. The slope of the rising light curve depends sensitively on the angle of the observer withrespect to the jet axis, which is likely to be perpendicular to the merger plane of the neutron starbinary. Therefore, the afterglow emission can be used to constrain the inclination of the mergingsystem. Here, we calculate the GRB afterglow emission based on the realistic jet structure derivedfrom general-relativistic magnetohydrodynamical simulations of black hole-torus system for the centralengine of the gamma-ray burst. Combined with gravitational wave parameter estimation, we fit themulti-epoch afterglow emission of GW170817. We show that with such a jet model, the observing anglecan be tightly constrained by multi-messenger observations. The best fit observing angle of GW170817is θ v = 0 . ± .
02 rad. With such a constraint, we can break the degeneracy between inclination angleand luminosity distance in gravitational wave parameter estimation, and substantially increase theprecision with which the Hubble constant is constrained by the standard siren method. Our estimationof the distance is D L = 43 . ± . ± − Mpc − .As a result, multi-messenger observations of short-duration gamma-ray bursts, combined with a goodtheoretical understanding of the jet structure, can be powerful probes of cosmological parameters. Keywords: gravitational waves, gamma-ray burst: individual: GRB 170817A, cosmological parameters INTRODUCTIONSince the detection of the neutron star merger event in gravitational waves (GW170817 Abbott et al. (2017a)),gamma-rays (GRB 170817A Goldstein et al. (2017)) and optical (AT2017gfo Tanvir et al. (2017)), short durationgamma-ray bursts became one of the most successful targets of multi-messenger astrophysics. Combining informationfrom both gravitational waves and electromagnetic waves, an impressive breadth of Astrophysical questions can beprobed, such as, for example, the origin of the heavy elements (Kasen et al. 2017), the neutron star equation ofstate (Abbott et al. 2018a) or the value of standard cosmological parameters (Abbott et al. 2017b). Recently theHubble constant has attracted much attention because of the tension between the value measured based on the cosmicmicrowave background and the cosmic distance ladder methods (Verde et al. 2019; Planck Collaboration et al. 2018;Riess et al. 2019). An alternative method to constrain the Hubble constant, the so-called standard siren method,takes advantage of the luminosity distance measured by gravitational wave signal and the cosmological redshift of thesource’s host galaxy (Nissanke et al. 2010). Such a method is independent to the previous methods, and may serve asan additional probe of the expansion rate of the Universe.The standard siren method requires an astrophysical source that provides bright electromagnetic counterparts alongwith a gravitational wave detection. One of the most promising sources is a binary neutron star merger which, not only
Corresponding author: Hao [email protected] a r X i v : . [ a s t r o - ph . H E ] S e p gives a powerful gravitational wave signal, but also expected to launch relativistic and non-relativistic outflows. Theoutflows can power a gamma-ray burst, afterglow emission and a kilonova (Metzger 2019). In this case, gravitational(electromagnetic) waves can be used to derive the luminosity distance (source location), respectively. However, theluminosity distance of a binary neutron star merger, as determined by the gravitational wave signal, can be stronglydegenerate with the orbital inclination of the system, that is to say from the gravitational signal alone, one is hardto distinguish whether a source is further away with the binary orbit facing Earth, or closer but the binary orbitbeen highly inclined to the line of sight. Such degeneracy results in rather large uncertainty of distance measurement,making the standard siren method harder to achieve a precision level comparable to other methods.One way in improving the accuracy of the method is by breaking the distance-inclination degeneracy. This can beachieved by including information from the gamma-ray burst afterglow observations. Several studies suggest that theafterglow emission that followed GRB 170817A originates from the interactions of a structured jet (i.e., a jet with itsproperties like energy and Lorentz factor smoothly varying as a function of polar angle) with the ambient gas (see,e.g., Ghirlanda et al. (2019)). In such a scenario, the slope of the rising flux, observed for ∼ § § §
4, including the implications of some parameters and how observing angleand luminosity distance can be better constrained. In § § AFTERGLOW MODEL2.1.
Jet Structure
The afterglow emission is powered by the shocks that the jet drives into the ambient gas. The afterglow light curvestrongly depends on the jet profile at large distance from the source where these interactions become substantial. Onthe other hand, the jet profile is determined at much more compact scale by the interactions of the jet while it breaks E () E ( ) 10 Figure 1.
The energy and Lorentz factor as a function of polar angle in our jet model. Energy is normalized to its peak value. out from slower gas expelled within a fraction of a second before and after the merger time. In this work we adopta jet structure derived from 3D GRMHD simulation where the jet is launched from a 3 solar mass central black holewith 0.8 spin surrounded by a compact torus of 0.03 solar mass embedded with a poloidal magnetic field. A detaileddescription of the setup can be seen from Fern´andez et al. (2019); Kathirgamaraju et al. (2019). At the initial stages ofthe simulation, the torus is characterized by slower wind-type outflows while mass, and magnetic flux, starts accretinginto the black hole. Once sufficient magnetic flux has accumulated through the black hole ergosphere, a powerful,relativistic jet forms. The jet is surrounded by the wind from the accretion disk; it interacts with it and colimatesout to distances of ∼ a few 1000 gravitational radii. After the jet breaks out, the jet turns conical and its structure isalmost frozen as a function of distance from the central engine. The simulation predicts a tightly collimated relativisticjet with an opening angle of ∼ .
2. The angular distribution of initial Lorentz factor and energy of the jet predictedby that work is shown in Figure 1. In the following section, we describe how one can calculate the afterglow emissiongiven the jet profile. 2.2.
Dynamical Evolution of the External Shock
The relativistically moving jet drives a blastwave into the ambient gas. Here, we work in the thin shell approximation,where we assume that the swept-up ambient matter is compressed in a thin region whose scale is much smaller thanthe radial scale. While the blastwave is initially ultra-relativistic, it gradually slows down and turns trans-relativisticwithin months after the burst. For a misaligned jet, with its core located at an angle θ v > / Γ with the line of sight,the emission from the core is initially beamed away from our line of sight. The peak of the afterglow emission takesplace when the Lorentz factor of the core of the blast slows down to Γ ∼ /θ v which falls into mildly relativistic domainif θ v is relatively large. Since we wish to follow the blast emission up to, and shortly after, its peak, a treatment thatfollows the blast deceleration down to a trans-relativistic speed is required. In this work we use a semi-analyticalmodel following the shock jump conditions derived by an approximate trans-relativistic EoS (Uhm 2011). Such modelworks well at both high and modest Lorentz factors and can smoothly transform to Sedov-Taylor solution when theshock is no longer relativistic.Our model assumes one-dimensional radial solution for each blast segment. Therefore, we neglect sideways (lateral)expansion in our model. The sideways expansion of the blast is important only at later stages of the blast evolutionand it affects the light curve predictions after its peak (Gill & Granot 2018). However, the key factor to constraininclination angle, of interest in this study, is the arising slope before the peak flux, where the jet is still relativistic. Sothe sideways expansion won’t influence our result.To describe the physical quantities more clearly, here we consider two reference frames. One is burst rest frame(hereafter BR frame) which is at rest with the central engine. The other is the shock co-moving frame (hereafter CMframe) which is at rest with respect to the gas just downstream of the forward shock. Primed variables (i.e., denotedwith (cid:48) ) are measured in CM frame, while non-primed one are measured in BR frame.Under the approximation of neglecting sideways expansion, we can regard the jet evolution at each radial direction asbeing independent. It is therefore convenient to describe the evolution at each direction by using isotropic equivalencequantities. Consider an isotropic burst of total energy E iso = Γ M ej c , where M ej is the mass of ejecta and Γ is itsinitial Lorentz factor. The outflow spreads in a cold uniform environment of number density n and forms a relativisticshock. We assume that the blastwave is expanding adiabatically without further energy injection. Energy loss due toradiation is negligible compared with E iso . We define the radius R of the shock by the distance from burst center tothe forward shock. All swept-up mass are compressed in the downstream of forward shock, moving with Lorentz factorΓ measured in BR frame, and its isotropic kinetic energy is E k . We also assume that all ejecta are moving togetherwith shocked ambient at same Lorentz factor behind the contact discontinuity and its total energy is E ej = Γ M ej c .The trans-relativistic shock jump conditions is (Uhm 2011; Ryan et al. 2019) n (cid:48) = 4Γ n (1) e (cid:48) th = 4Γ(Γ − nm p c (2) p (cid:48) = 4(Γ − nm p c / R = 4 βc Γ / (4Γ − , (4)where n (cid:48) , e (cid:48) th , p (cid:48) are the number density, internal energy and pressure of shocked region, respectively. β = (1 − Γ − ) / is the dimensionless speed. c is the speed of light. Because all swept-up matter are gathered in the shell, one can solvethe width of the shell by equating 4 πR n/ πR ∆ R ( n (cid:48) Γ). Applying Eq. 1 we have∆ R = R/ (12Γ ) . (5)The kinetic energy is therefore (Uhm 2011; Ryan et al. 2019; van Eerten 2013) E k = (cid:2) ( ρ (cid:48) c + e (cid:48) th + p (cid:48) )Γ − p (cid:48) − Γ ρ (cid:48) c (cid:3) πR ∆ R = 4 πR nm p β (4Γ − c / , (6)where ρ (cid:48) = n (cid:48) m p is the density of shocked region.As the shock spreads, the ejecta gives its kinetic energy to swept-up matter while keeping the total energy conserved.i.e. E iso = E k + E ej = const. Substituting M ej in terms of E iso and Γ , the total energy is E iso = 11 − Γ / Γ π R nm p β (4Γ − c . (7)The shock dynamical evolution can then be solved by combining Eq 4 and Eq 7. The solution shows three partsof the evolution. Initially the coasting phase where most of the energy is carried by the ejecta and its Lorentz factorremains almost constant. When the accumulated ambient mass is sufficiently large, the shock enters the decelerationphase, where Γ ∝ R − / , while R ∼ ct because the shock is still ultra-relativistic. When the shock is no longerrelativistic, we reproduce that R ∝ t / as expected from the Sedov-Taylor solution.2.3. Synchrotron Radiation
The accelerated, non-thermal electrons in the shock produce synchrotron radiation. In the standard GRB af-terglow model, the non-thermal electron Lorentz factor distribution can be characterized by a power-law function n (cid:48) e = Cn (cid:48) γ − pe , γ e > γ m , where γ m is the minimum electron Lorentz factor and C is a normalization factor. C and γ m can be solved assuming that a fraction (cid:15) e of internal energy is transferred to non-thermal electrons (Sari et al. 1998).In this work we also use an additional free parameter f which describes the fraction of accelerated electrons to thetotal number in the swept-up region, i.e. n (cid:48) e = Cf n (cid:48) γ − pe (van Eerten 2013).With such modification, the electron minimum Lorentz factor is γ m = p − p − (cid:15) e m p f m e (Γ −
1) (8)Assume a fraction (cid:15) B part of internal energy transferred to that of the magnetic field, the magnetic field thereforeis B (cid:48) = (cid:113) πe (cid:48) th (cid:15) B = (cid:113) π Γ(Γ − nm p (cid:15) B c . (9)The cooling Lorentz factor of the electrons is found by equating the synchrotron loss to the expansion time of theblast, resulting in γ c = 6 πm e Γ cσ T B (cid:48) t , (10)where t is the time since the burst measured in BR frame.For each characteristic Lorentz factor γ i = γ m , γ c , the corresponding frequency is ν (cid:48) i = 3 eB (cid:48) γ i πm e c . (11)The synchrotron radiation peak emissivity is (cid:15) (cid:48) P = √ e B (cid:48) f n (cid:48) m e c . (12)The broad band observation of GRB 170817A from radio to x-ray shows a constant spectral index (Mooley et al.2018a; Margutti et al. 2018). Such spectrum is consistent with optically thin synchrotron radiation. . The Spectrumis then (cid:15) (cid:48) ν (cid:48) = (cid:15) (cid:48) P × ( ν (cid:48) /ν (cid:48) m ) / ν (cid:48) < ν (cid:48) m < ν (cid:48) c ( ν (cid:48) /ν (cid:48) m ) − ( p − / ν (cid:48) m < ν (cid:48) < ν (cid:48) c ( ν (cid:48) c /ν (cid:48) m ) − ( p − / ( ν (cid:48) /ν (cid:48) c ) − p/ ν (cid:48) m < ν (cid:48) c < ν (cid:48) ( ν (cid:48) /ν (cid:48) m ) / ν (cid:48) < ν (cid:48) c < ν (cid:48) m ( ν (cid:48) c /ν (cid:48) m ) / ( ν (cid:48) /ν (cid:48) c ) − / ν (cid:48) c < ν (cid:48) < ν (cid:48) m ( ν (cid:48) c /ν (cid:48) m ) / ( ν (cid:48) m /ν (cid:48) c ) − / ( ν (cid:48) /ν (cid:48) m ) − p/ ν (cid:48) c < ν (cid:48) m < ν (cid:48) . (13)2.4. Observed Flux of a structured jet
In order to calculate the total observed flux, we need to integrate over the whole solid angle. It is convenientto consider a spherical coordinate system where the origin is the center of the burst and z axis is the axis of jet.It is also natural to set the observer’s direction as ϕ = 0. In the structured jet case, Lorentz factors, radius andemissivity are functions of θ and t and can be solved given the initial condition at angle θ . The isotropic equivalent For consistency, we have included the synchrotron self-absorption (SSA) in our calculation, but the SSA break is below the observed band,so for simplicity, we do not include SSA here. energy at each direction can be calculated by a free parameter E tot which is the angle integrated jet total energy. i.e. E iso ( θ ) = 4 πE ( θ ) E tot / (cid:82) E ( θ )d θ . The isotropic equivalent radiation power at coordinate ( θ, ϕ ) is then P iso ( θ, ϕ, t ) = (cid:15) (cid:48) ν (cid:48) ( θ, t )4 πR ( θ, t ) ∆ R (cid:48) ( θ, t ) δ ( θ, ϕ, t ) , (14)where ∆ R (cid:48) = ∆ R Γ is the shock width measured in CM frame, δ = [Γ( θ, t )(1 − β ( θ, t ) µ ( θ, ϕ )] − is Doppler factor, µ isthe cosine angle between ( θ, ϕ ) and the observer’s direction µ = cos θ cos θ v + sin θ cos ϕ sin θ v . (15)The observed flux at frequency ν obs and observer’s time T obs since burst is therefore F ν obs ( T obs ) = 1 + z πD L (cid:90) P iso ( θ, ϕ, t )4 π dΩ , (16)where z is redshift and D L is the luminosity distance.The frequency at CM frame can be calculated through Doppler factor ν (cid:48) = ν obs (1 + z ) /δ . The time t at BR frameis numerically solved by equal-arrival-time-surface provided θ , ϕ and T obs t − R ( θ, t ) µ ( θ, ϕ ) /c = T obs / (1 + z ) . (17)2.5. Light Curve Behavior
In Fig. 2, we plot a collection of light curves of the afterglow emission predicted by our model for a jet observedat different angles. At 10 Hz (solid lines), the light curves are characterized by two breaks. The early-time breakrepresents the transition from coasting stage to deceleration stage of the segment of the jet that dominates the early-time emission. After that stage, because of the beaming effect as the jet slows down, observers can see a graduallyincreasing area closer to the core of the jet. The slope of the light curve at this stage strongly depends on θ v , and is akey quantity that we constrain with the fits it in our work. When the observer sees the emission from the jet core, themost powerful segment of the jet contributes, and dominates, the observed flux and the light curves hereafter have asecond break.At 3GHz, the light curves of low θ v at early time show a feature of increasing flux, which is different from onesfrom same angle but higher frequency. Such feature originates from “spectral break”, i.e., the source frequency ν (cid:48) = ν obs (1 + z ) /δ is lower than ν (cid:48) m because of the high Doppler factor. PARAMETER ESTIMATIONIn this section, we apply our model to the radio data of GRB 170817A and perform a Bayesian parameter estimation.The data points that we use for our fits are gathered from a number of publications (Hallinan et al. 2017; Mooley et al.2018a; Margutti et al. 2018; Alexander et al. 2018). Higher energy data like in the optical and X-ray bands followthe power-law extrapolation of the radio data (Margutti et al. 2018). As such, they can help to constrain electronpower-index p , but do not provide additional information to constrain θ v and is unnecessary in our work. We do,however, check for self consistency that our best fit models also account for the X-ray observations (e.g., the coolingfrequency is not crossing the X-ray band during the observational window; see § p = 2 .
17 in our fitting. We also fix the source cosmological redshift to z = 0 . ν obs in our afterglow model are then ( E tot , θ v , n , (cid:15) e , (cid:15) B , D L , f ). We calculate the fitting goodness by χ and the likelihood function is then L = exp( − χ / . We have used the cleaned version where the glitch has been removed. Thealgorithm is described in Ashton et al. (2019) and the priors are chosen to the same as Abbott et al. (2019a), but with T obs / Day10 F o b s / J y v = 0.03 v = 0.1 v = 0.2 v = 0.3 v = 0.4 Figure 2.
The light curves at different observing angles generated by our model. The parameter set is E tot = 10 erg, n = 10 − cm − , (cid:15) e = 0 . (cid:15) B = 10 − , D L = 40Mpc, f = 1. The solid (dashed) lines show the 10 Hz (3GHz) light curves,respectively. Angles are measured in rad. sky localization fixed to it’s host galaxy NGC4993. A fixed coordinate will help reduce luminosity distance uncertainty.The waveform template we used here is IMRPhenomPv2 NRTidal (Dietrich et al. 2017, 2019; Hannam et al. 2014).The posterior distribution of luminosity distance and inclination angle then serve as a prior to the correspondingparameters in our model, where we have assumed that the inclination angle is exactly the observing angle (i.e., thatthe jet is ejected perpendicular to the plane of the merger). To deal with the facing-off cases (i.e., the inclination angleis larger than π/ π − θ v in the afterglow calculation.To show how model parameters work in the radio data fitting and how we combine our model with gravitationalwaves, we did three fitting with different prior sets. They are summarized in Table 1. Below are the description ofthose priors. Prior set 1.:
This is the simplest afterglow fit setting where the fraction of accelerated electrons f is set to 1 anddistance D L is fixed to 40 Mpc. Other parameters are uniformly distributed in bounds. Prior set 2.:
In this prior set we leave f free to show the importance of this parameter in our fit, while the distanceis still fixed. Prior set 3.:
We leave all parameters free, while the prior of observing angle and luminosity distance is a kerneldensity estimation function built from the posterior samples of GW170817 in our analysis above. Applying thisprior is equivalent to a joint parameter estimation of GRB 170817A and GW170817.
Table 1.
The prior distribution in parameter estimation. All the distributions are uniform in the given bound. GW refers tothe kernel density estimation built from posterior samples of gravitational waves.Sampling Parameters Prior Set 1 Prior Set 2 Prior Set 3log( E tot / erg) [47, 55] [47, 55] [47, 55]cos θ v [0, 1] [0, 1] GWlog( n/ cm − ) [-6, 0] [-6, 0] [-6, 0]log( (cid:15) e ) [-3, 0] [-3, 0] [-3, 0]log( (cid:15) B ) [-6, 0] [-6, 0] [-6, 0] D L / Mpc 40 40 GWlog( f ) 0 [-4, 0] [-4, 0] For each prior setup, we implement a nested sampling algorithm to generate posterior samples, employing the openpackage Pymultinest (Buchner et al. 2014). We plot the parameter estimation results of different prior sets in Figure 3,4 and 5 respectively. In addition, we draw the best fitting light curve of prior set 3 together with radio data normalizedto 3GHz in Figure 6. IMPLICATIONS4.1.
Model Parameters
The fitting process provides intriguing constraints to the parameters. All our fittings result in a tight constraint onthe observing angle θ v = 0 . ± .
02 radian, which arises from the sensitive relation between the angle and light curveslope. Our inferred value for θ v is consistent with the angle estimated by the jet’s superluminal motion (Mooley et al.2018b).It is not surprising that we can’t precisely constrain the magnetic field parameter (cid:15) B , because the cooling frequency ν c is not observed in radio band. In fact, a further constraint on ν c can be placed since the spectrum shows a singlepower-law function from the radio to x-ray band, which means the cooling break lies beyond the X-rays (Marguttiet al. 2018). The constraints of the afterglow parameters reported here (and the best fitting parameter in Figure 6)are compatible with these studies. The parameter ν c in our study therefore only plays a role of a scale factor in thecalculation and doesn’t affect our estimation of θ v . However it still helps us to extract some information of electronacceleration process which is implied in Figure 3 and 4. Assuming that (cid:15) B is not too high ( (cid:15) B < .
01, see the X-rayconstraint from Hajela et al. (2019)), the (cid:15) e − (cid:15) B contour plot in Figure 3 shows a very large estimation of (cid:15) e ( > . (cid:15) e results from the early time observation when the frequency is possibly near orlower than ν m (e.g. The radio data of VLA at ∼ . ν m at very early time doesn’t necessarilymean that a high amount of internal energy is transferred to non-thermal electrons. If we relax the assumption thatall electrons in the shock region are accelerated (i.e., setting the fraction parameter f as a free parameter), the fittingresults in Figure 4 show that a high estimation of (cid:15) e is no longer required, although the constraint is much loosenbecause of the strong degeneracy between f and (cid:15) e , as we can see from Eq 8. It can be regarded as an evidence showingthat f = 1 is disfavored. This parameter, however, doesn’t affect the observing angle estimation.The fitting after leaving f free implies an outflow with initial kinetic energy of E tot ≈ − erg in a low densityenvironment of n ≈ − − − cm − . These results are compatible with the constraints placed on the ambient densityby observations of the jet superluminal motion (Mooley et al. 2018b) and X-ray emission in the host galaxy (Hajelaet al. 2019). 4.2. Breaking the D L -inclination degeneracy The tight constraint of θ v implies that, if we assume θ v being the same as the binary inclination (i.e., the jet axis isperpendicular to the merger plane), it can be used to break the degeneracy between luminosity distance and inclination Note that we didn’t include the Synchrotron Self-Compton effect in this work, which may have a modest effect to the electron cooling break(Beniamini et al. 2015). +0.480.43 l o g ( E t o t / e r g ) +0.430.41 . . . v +0.020.02 l o g ( e ) +0.200.21 log( n /cm ) l o g ( B ) log( E tot /erg) . . . v log( e ) log( B )1.79 +1.221.26 Figure 3.
Parameter estimation result with prior set 1, where f = 1 and D L = 40 Mpc. The contours with different colorsrepresent 1 σ − σ uncertainty levels. angle from the gravitational wave parameter estimation. After applying prior set 3, we can see from Figure 5 thatluminosity distance can be tightly constrained to D L = 43 . ± (Abbott et al. 2019b). Benefited from a good estimation of θ v , The uncertainty of D L in ourresult has been reduced by a factor of 7 compared with pure gravitational wave estimation, and by a factor of 4 if https://dcc.ligo.org/LIGO-P1800370/public +0.930.88 l o g ( E t o t / e r g ) +0.940.85 . . . v +0.020.02 l o g ( e ) +0.960.95 l o g ( B ) +2.101.93 log( n /cm ) l o g ( f ) log( E tot /erg) . . . v log( e ) log( B ) log( f )1.28 +0.820.83 Figure 4.
Parameter estimation result with prior set 2, where f is left free while D L is still fixed to 40 Mpc. host galaxy localization is added. Such result shows that taking advantage of multi-messenger parameter estimation, asimulated jet structure model with a small opening angle uncertainty can significantly increase the precision of distancemeasurement. 4.3. Hubble Constant
Our tight constraint of luminosity distance which is independent of the cosmological model can be used to infercosmology parameters. Especially, at the nearby universe, it is possible to set a better constraint of the Hubbleconstant compared with the standard siren method using gravitational waves data alone. For the case of GW170817which lies fairly close, the Hubble constant measurement is insensitive to other cosmological parameters such as Ω m +0.940.89 l o g ( E t o t / e r g ) +0.940.88 . . . v +0.020.02 l o g ( e ) +0.940.96 l o g ( B ) +2.081.99 l o g ( f ) +0.810.84 log( n /cm ) D L / M p c log( E tot /erg) . . . v log( e ) log( B ) log( f ) D L /Mpc43.38 +1.201.15 Figure 5.
Parameter estimation result with prior set 3, where all parameters are left free. The prior of θ v and D L follows Table1. and Ω Λ . It is also unnecessary to distinguish between different cosmological distance scales, since their difference is atthe order of v H /c ∼
1% at such a distance, where v H is the Hubble flow. We then follow the same method describedin Abbott et al. (2017b), where the following expression provides a good approximation v H = H · D. (18)To make a comparison with standard siren method, we use the same Hubble flow v H = 3017 ±
166 km s − (Crooket al. 2007) as in Abbott et al. (2017b), and assume it follows a normal distribution. The result is shown in Figure 8,where we also plot the 1 σ confidence interval of results from Planck Mission (Planck Collaboration et al. 2018) and2 T obs / Day10 F o b s / J y Best Fitting Curve1 uncertainty of v VLA Data (3GHz)
Figure 6.
The light curve fitting of the afterglow of GRB 170817A with 1 σ uncertainty region of θ v . All data points arenormalized to 3GHz. SHOES (Riess et al. 2019). We constrain the Hubble constant to H = 69 . ± − Mpc − and the uncertainty hasreduced by a factor of more than 2 by the inclusion of the radio constraints. It is worth mentioning that half of theuncertainty comes from the uncertainty in the peculiar velocity of the host galaxy in the Hubble flow which potentiallycould be constrained better in future observations. PROSPECTS OF FUTURE DETECTIONSThe prospects of multi-messenger joint fitting of future gravitational wave events strongly depend on whether EMcounterparts can be observed. Especially, in our model, the key factor is the detectability of the afterglow. To quantifythe prospects of afterglow detection of future binary neutron star merger events, we perform a Monte Carlo simulationthat generates a million of candidate events uniformly distributed in the universe, and pick up those who can bedetected by LIGO and Virgo. The method, which for simplicity focuses to the radio afterglow emission, is describedbelow.We generate samples of merger events produced by equal massive binary neutron stars M = 1 . M (cid:12) without spin.The samples are uniformly located in the volume up to a distance of 400 Mpc. Since the signal-to-noise-ratio(SNR) ofa binary neutron star merger is dominated by it’s inspiral stage, for simplicity, we adopt gravitational wave templateTaylorF2 (Vines et al. 2011) which provides good approximation at this stage. The detectors’ sensitivity curve is3 D L = 39.96 +5.399.12 D L = 40.73 +3.855.95 D L = 43.38 +1.21.15 D L / Mpc . . . v . . . vv = 0.58 +0.330.29 v = 0.52 +0.230.22 v = 0.38 +0.020.02 GW OnlyGW + Host GalaxyGW + Host Galaxy + Radio
Figure 7.
Contour plot of posterior probability distribution of θ v and D L of GW170817. The green contour is from puregravitational wave result. The blue contour is also from gravitational wave but with sky localization fixed to it’s host galaxy.The red contour is the joint fit using our model with radio afterglow observation and gravitational waves. It’s the same contouras in Figure 5. Inner and outer contour represents 1 σ and 2 σ uncertainty levels respectively. Note that the EM data greatlyreduce the uncertainty in both θ v and D L .
50 60 70 80 90 100 110 120 130 140 H (km s Mpc )0.000.020.040.060.080.10 P r o b a b ili t y D e n s i t y H = 74.51 +12.927.72 H = 69.48 +4.324.21 GW + Host GalaxyGW + Host Galaxy + RadioPlanckSHOES
Figure 8.
Hubble constant estimation using luminosity distance posteriors from our joint fitting. The green region representsthe measurement from Planck Mission(Planck Collaboration et al. 2018) and the orange one is from SHOES(Riess et al. 2019).All uncertainty levels here are at 1 σ . The uncertainty in reduced by a factor of 2-3 by the inclusion of the radio afterglowconstraints. adopted from the designed sensitivity of LIGO’s public document , and we use it as an approximation to the O4sensitivity. In the simulation an event is regarded as a confirmed detection if the LIGO and Virgo network SNR islarger than 12. We use open source package Bilby (Ashton et al. 2019) to perform signal injection and SNR calculation.To test if those events can produce detectable radio emission, we use our model to generate light curves withparameters same to the best fitting curve in Figure 6, except that distances and observing angles are free. Such setupcan well estimate the flux of GRB 170817A-like events that are randomly distributed and orientated in the universe.The radio telescopes’ detection limit is somewhat hard to estimate due to calibration issues. Considering that the firstfew radio data points of GRB 170817A detected by VLA are about 10 − µ Jy (Hallinan et al. 2017), we set two fluxdensity limit: 10 µ Jy and 20 µ Jy in our calculation. Another limitation to the radio telescopes is their response time,which depends on the progress of searching for host galaxy. Here we set the initial time to be 1 day which is similar tothe case of GW170817. Given the above setup, we can calculate the maximum distance for the detection of an eventas a function of it’s observing angle. Here we assume that a successful detection happens when the peak flux densityexceeds the telescope’s detection limit. https://dcc.ligo.org/LIGO-P1200087-v47/public v (rad)50100150200250300350400 D L ( M p c ) R a d i o D e t e c t a b l e . % v (deg) O4 Events Density10 Jy Limit at 3GHz20 Jy Limit at 3GHzGW170817GW190425 68.3% CI
Figure 9.
Prospects of O4 multi-messenger detection. The simulated gravitational wave events are represented by blue densityplot, where deeper color means high gathering event numbers. The red lines are maximum detection distance assuming 10 µ Jyor 20 µ Jy flux density limit of radio telescopes. Samples located on their low-left side can produce detectable radio afterglowsin our model.
Our result is shown in Figure 9. The simulated events number density at each coordinate point are represented bythe color depth. For those facing-off events (i.e. θ v > π/ π − θ v to match theobserving angle.Our result shows that a considerable amount of events can be detected at larger distance even beyond the averagehorizontal distance (Abbott et al. 2018b) because of their optimal orientation. Our simulation estimates that detectionrate of binary neutron star merger is approximately R GW = 20 (cid:16) R BNS − yr − (cid:17) yr − , where R BNS is the binary neutronstar merger volumetric rate. Among those detections, according to different flux density limit setup, 10% −
15% mayproduce detectable radio counterparts at their peak emission. Using the latest estimation of R BNS (Abbott et al. 2020),those numbers imply that we may have a few electromagnetic counterparts during LIGO O4. However, we should notethat further upgrades of gravitational wave detectors’ sensitivity may not help to substantially increase the number ofelectromagnetic counterparts, because the flux of faraway off-axis GRB afterglows will drop below detection thresholdof any current radio telescope.Radio emission is not the only EM counterpart expected in these sources. Other work has shown that the promptemission of GRBs lying in radio detectable region are also likely detectable (Kathirgamaraju et al. 2018; Beniaminiet al. 2019), just as the case of GRB 170817A. Optical radiation from kilonovae whose radiation is more isotropic alsohelps from localization. The estimated fraction of events with EM counterparts is therefore robust.6Considering the telescope response time caused by searching, we may still not be able to constrain the observing angleas well as GW170817 in events that are more distant. As we can see from Figure 9, most events with EM counterpartshave relatively low observing angle of 0 . − . − CONCLUSION AND DISCUSSIONThe blastwave driven by the GRB jet is highly anisotropic. As a result, the light curves of GRB afterglows stronglydepend on the observing angle and can thus be used to constrain it. Such constraint requires an accurate model for thejet structure. In this work, we use the structured jet model predicted by 3D GRMHD simulations of a black hole, torussystem for the central engine of neutron star mergers. We apply this model to gravitational and electromagnetic waveobservations of GW170817, leading to a very tightly constrained observing angle for the system of θ v = 0 . ± . D L = 43 . ± H = 69 . ± − Mpc − .We also use Mote Carlo simulation to explore the prospects of applying this method to future detection of mergerevents. Assuming the current estimation of binary neutron star merger rate of ∼ − yr − , there may be afew neutron star merger events per year with electromagnetic counterpart observed in LIGO’s O4 stage (out of ∼ ∼
100 Mpc to 300 Mpc. In a fair fraction of these events, our line of sight will be atthe edge of the jet core, which is a different setup in comparison to the more inclined GRB 170817A. Such occurrencemay increase the difficulty of accurately constraining observing angle, but the constraints will still be much betterthan using gravitational waves only. With several such events, the Hubble constant measurement in this method mayachieve a higher precision and help to resolve the existing tension that appears among other methods.Our results imply that the inclination angle of the merger can be very well constrained, provided that we have areliable model for the jet. The inference of the observing angle is most sensitive to the geometry, especially the openingangle of the jet, which can rely on assumptions for the nature of the central engine. For the case of binary neutronstar merger, considering the current constraint of neutron star maximum mass M NS , max ≈ . M (cid:12) (Margalit & Metzger2017), it is reasonable to expect that the merger remnant of GW170817 (and of most binary neutron star mergers)collapses fast into a spinning black hole. Material with large amounts of specific angular momentum will surround theblack hole forming a geometrically thick accretion disk or torus. The accretion of the torus can result in magnetic fluxaccumulation onto the black hole. In a such system (i.e. a magnetized black hole with high spin), the Blanford-Znajekprocess seems to be responsible for launching a Poynting-flux dominated jet. The jet is then collimated by denseoutflows from the disc until it breaks out from the surrounding gas. Here, we adopted the Fern´andez et al. (2019)GRMHD setup in simulating such system. Although, clearly representing the state of the art, that work is limited toone model for the central engine. Clearly, a large number of realistic simulations with different initial conditions, toaccount, e.g., for the different progenitor masses and mass ratios are needed to investigate whether the findings on thejet structure are robust. Independent work (Beniamini et al. 2019), however, indicates a relatively narrow range ofshort GRB opening angle, and its independence to the luminosity function, supporting the promise of our approach.ACKNOWLEDGMENTSWe acknowledge support from the NASA ATP NNX17AG21G, the NSF AST-1910451 and the NSF AST-1816136grants.7REFERENCES Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a,PhRvL, 119, 161101,doi: 10.1103/PhysRevLett.119.161101—. 2017b, Nature, 551, 85, doi: 10.1038/nature24471—. 2018a, PhRvL, 121, 161101,doi: 10.1103/PhysRevLett.121.161101—. 2018b, Living Reviews in Relativity, 21, 3,doi: 10.1007/s41114-018-0012-9—. 2019a, Physical Review X, 9, 011001,doi: 10.1103/PhysRevX.9.011001—. 2019b, Physical Review X, 9, 031040,doi: 10.1103/PhysRevX.9.031040—. 2020, ApJL, 892, L3, doi: 10.3847/2041-8213/ab75f5Alexander, K. D., Margutti, R., Blanchard, P. K., et al.2018, ApJL, 863, L18, doi: 10.3847/2041-8213/aad637Ashton, G., H¨ubner, M., Lasky, P. D., et al. 2019, ApJS,241, 27, doi: 10.3847/1538-4365/ab06fcBeniamini, P., Granot, J., & Gill, R. 2020, MNRAS, 493,3521, doi: 10.1093/mnras/staa538Beniamini, P., Nava, L., Duran, R. B., & Piran, T. 2015,MNRAS, 454, 1073, doi: 10.1093/mnras/stv2033Beniamini, P., Petropoulou, M., Barniol Duran, R., &Giannios, D. 2019, MNRAS, 483, 840,doi: 10.1093/mnras/sty3093Beniamini, P., & van der Horst, A. J. 2017, MNRAS, 472,3161, doi: 10.1093/mnras/stx2203Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433,doi: 10.1093/mnras/179.3.433Buchner, J., Georgakakis, A., Nandra, K., et al. 2014,A&A, 564, A125, doi: 10.1051/0004-6361/201322971Crook, A. C., Huchra, J. P., Martimbeau, N., et al. 2007,ApJ, 655, 790, doi: 10.1086/510201Dietrich, T., Bernuzzi, S., & Tichy, W. 2017, PhRvD, 96,121501, doi: 10.1103/PhysRevD.96.121501Dietrich, T., Khan, S., Dudi, R., et al. 2019, PhRvD, 99,024029, doi: 10.1103/PhysRevD.99.024029Fern´andez, R., Tchekhovskoy, A., Quataert, E., Foucart,F., & Kasen, D. 2019, MNRAS, 482, 3373,doi: 10.1093/mnras/sty2932Ghirlanda, G., Salafia, O. S., Paragi, Z., et al. 2019,Science, 363, 968, doi: 10.1126/science.aau8815Gill, R., & Granot, J. 2018, MNRAS, 478, 4128,doi: 10.1093/mnras/sty1214Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJL, 848,L14, doi: 10.3847/2041-8213/aa8f41Gottlieb, O., Nakar, E., & Piran, T. 2019, MNRAS, 488,2405, doi: 10.1093/mnras/stz1906Hajela, A., Margutti, R., Alexander, K. D., et al. 2019,ApJL, 886, L17, doi: 10.3847/2041-8213/ab5226 Hallinan, G., Corsi, A., Mooley, K. P., et al. 2017, Science,358, 1579, doi: 10.1126/science.aap9855Hannam, M., Schmidt, P., Boh´e, A., et al. 2014, PhRvL,113, 151101, doi: 10.1103/PhysRevLett.113.151101Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, NatureAstronomy, 3, 940, doi: 10.1038/s41550-019-0820-1Kasen, D., Metzger, B., Barnes, J., Quataert, E., &Ramirez-Ruiz, E. 2017, Nature, 551, 80,doi: 10.1038/nature24453Kathirgamaraju, A., Barniol Duran, R., & Giannios, D.2018, MNRAS, 473, L121, doi: 10.1093/mnrasl/slx175Kathirgamaraju, A., Tchekhovskoy, A., Giannios, D., &Barniol Duran, R. 2019, MNRAS, 484, L98,doi: 10.1093/mnrasl/slz012Kumar, P., & Zhang, B. 2015, PhR, 561, 1,doi: 10.1016/j.physrep.2014.09.008Margalit, B., & Metzger, B. D. 2017, ApJL, 850, L19,doi: 10.3847/2041-8213/aa991cMargutti, R., Alexander, K. D., Xie, X., et al. 2018, ApJL,856, L18, doi: 10.3847/2041-8213/aab2adMetzger, B. D. 2019, Living Reviews in Relativity, 23, 1,doi: 10.1007/s41114-019-0024-0Mooley, K. P., Nakar, E., Hotokezaka, K., et al. 2018a,Nature, 554, 207, doi: 10.1038/nature25452Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018b,Nature, 561, 355, doi: 10.1038/s41586-018-0486-3Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., &Sievers, J. L. 2010, ApJ, 725, 496,doi: 10.1088/0004-637X/725/1/496Planck Collaboration, Aghanim, N., Akrami, Y., et al.2018, arXiv e-prints, arXiv:1807.06209.https://arxiv.org/abs/1807.06209Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., &Scolnic, D. 2019, ApJ, 876, 85,doi: 10.3847/1538-4357/ab1422Ryan, G., van Eerten, H., Piro, L., & Troja, E. 2019, arXive-prints, arXiv:1909.11691.https://arxiv.org/abs/1909.11691Sari, R., Piran, T., & Narayan, R. 1998, ApJL, 497, L17,doi: 10.1086/311269Tanvir, N. R., Levan, A. J., Gonz´alez-Fern´andez, C., et al.2017, ApJL, 848, L27, doi: 10.3847/2041-8213/aa90b6Troja, E., Piro, L., Ryan, G., et al. 2018, MNRAS, 478,L18, doi: 10.1093/mnrasl/sly061Troja, E., van Eerten, H., Ryan, G., et al. 2019, MNRAS,489, 1919, doi: 10.1093/mnras/stz2248Uhm, Z. L. 2011, ApJ, 733, 86,doi: 10.1088/0004-637X/733/2/868