Probing binary neutron star mergers in dense environments using afterglow counterparts

The only binary neutron star merger gravitational wave event with detected electromagnetic counterparts recorded to date is GRB170817A. This merger occurred in a rareﬁed medium with a density smaller than 10 − 3 − 10 − 2 cm − 3 . Since kicks are imparted to neutron star binaries upon formation, and due to their long delay times before merger, such low-density circum-merger media are generally expected. However, there is some indirect evidence for fast-merging or low-kick binaries, which would coalesce in denser environments. Nonetheless, present astronomical data are largely inconclusive on the possibility of these high-density mergers. We describe a method to directly probe this hypothetical population of high-density mergers through multi-messenger observations of binary neutron star merger afterglows, exploiting the high sensitivity of these signals to the density of the merger environment. This method is based on a sample of merger afterglows that has yet to be collected. Its constraining power is large, even with a small sample of events. We discuss the method’s limitations and applicability. In the upcoming era of third-generation gravitational wave detectors, this method’s potential will be fully realized as it will allow us to probe mergers that occurred soon after the peak of cosmic star formation, provided the follow-up campaigns are able to locate the sources.


Introduction
Upon the second supernova leading to their formation, binary neutron stars (BNS) are kicked away from their dense star-forming birth regions (Blaauw 1961;Boersma 1961;Fryer & Kalogera 1997), allowing them to migrate to a different environment before merging (Portegies Zwart & Yungelson 1998).Because of the slow rate of orbital decay, this migration is generally expected to be long enough for the merger to occur in a rarefied medium (e.g.Bloom et al. 1999).This was the case for the up-to-now single BNS merger gravitational wave (GW) event with afterglow counterpart, GRB170817A, which occurred in a medium with density n 10 −2 cm −3 (Hallinan et al. 2017;Hajela et al. 2019, and see Sect. 3.1 below).
Supernova kicks are poorly constrained in the general picture (Podsiadlowski et al. 2005) and may be variable from a system to another (e.g., Podsiadlowski et al. 2004).Overall, they define the system's velocity during migration and affect its initial separation and eccentricity and therefore its merger time (Brandt & Podsiadlowski 1995;Kalogera 1996;Belczyński & Bulik 1999).For systems with the lowest kick velocities or shortest delay-times, we expect the distances covered during migration to be shorter than for the rest of the population, leading to the possibility of binaries merging in environments that are much denser than those encountered by systems with long migrations.We refer to these events, with densities n 1 cm −3 , as "high-density mergers".
Furthermore, there exist some theoretical stellar evolutionary pathways leading to short-delay or low-kick systems and therefore high-density mergers (e.g., Ivanova et al. 2003;Liu & Lai 2018;Secunda et al. 2019).Whether these mechanisms are realized is not certain and, as we show below, current data is inconclusive regarding the importance of this class of mergers.
In the electromagnetic domain, modelling of the afterglows of short gamma-ray bursts (GRB) is one probe of the environments of BNS mergers (e.g., Fong et al. 2015).Unfortunately, because of the poor localization of most short GRBs and of the relative faintness of their afterglows, the X-ray afterglow of only a small fraction have been found, and less than a handful have detected afterglows in the radio band (D'Avanzo 2015).
In recent works, Gottlieb et al. (2019) and Duque et al. (2019) have studied the afterglows expected as counterparts to GW signals from BNS mergers in present and future observing runs.Starting from a population model motivated by knowledge of short GRB, Duque et al. (2019) have found that the fraction of afterglows detectable in the radio band sharply increases with the density n of the medium hosting the mergers.This is due to the fact that (i) radio frequencies ν R are expected to fall between the injection and cooling frequencies ν i and ν c of the synchrotron slow-cooling regime for the bulk of the population, and (ii) in this regime, the afterglow peak flux scales as F p ∝ n p+1 4 (e.g., Nakar et al. 2002), where p ∼ 2−3 is the spectral index of the non-thermal electron population accelerated at the jet's forward shock front.Thus, should there be mergers in high-density environments, these would be over-represented in the afterglow population with respect to their actual number.In other words, the radio afterglow acts as an amplifier for these higher density mergers.
Given a statistical flux-limited sample of BNS merger afterglow counterparts endowed with sufficient completeness in circum-merger density estimates, one can determine the apparent fraction of high-density mergers.Starting from this number, by estimating the amplification factor related to the high-densityselection effect from population models, one can constrain the intrinsic fraction of mergers in high-density media.This is the principle of the new method we propose in order to study the class of high-density mergers.
As we develop later on, this method should allow us to constrain the number of high-density mergers, even after a small number of GW events with afterglow counterpart.The exact link between the rate of high-density events and the distribution of delay times and kick velocities is not clear, in particular because of the aforementioned uncertainty on the supernova kicks.Nonetheless, the method we suggest here is a first step toward studying the delay-time distribution of BNSs from their merger afterglows, at least for the fastest merging or low-kick binaries.
This article is organized as follows.In Sect.2, we explain the motivation for developing this new method of constraining highdensity mergers by recalling some related observational and theoretical knowledge.In Sect.3, we show that, for future BNS merger afterglows, multi-messenger observations will allow the circum-merger density to be estimated and, thus, will enable the apparent fraction of high-density mergers to be determined quite accurately.In Sect.4, we describe how these observations provide significant constraints on the population of high-density mergers even with a limited sample of afterglows, exploiting the sensitivity of the afterglow flux to the circum-merger density.Finally, in Sect. 5 we discuss the limitations of this method, and conclude in Sect.6.

Indirect evidence regarding BNS mergers in dense media
Theoretically, mechanisms exist that lead to fast-merging or lowkick systems.Among these are (i) an efficient common envelope phase, that reduces initial separation (e.g., Ivanova et al. 2003;Dominik et al. 2012) and merger time, (ii) a favorable supernova kick, that causes high eccentricity and thus rapid merger or a small migration velocity (e.g., Kalogera 1996), (iii) the formation of the BNS by dynamical capture in a migration trap within an active galactic nucleus disk (Secunda et al. 2019), or (iv) the interaction of the BNS with another compact object therein (Liu & Lai 2018;Fernández & Kobayashi 2019).The frequency with which these actually occur is still unclear.
Over the years, a body of indirect evidence on high-density mergers has emerged.However, as we show here, current data is inconclusive regarding the importance of this class of mergers.
First, some population synthesis studies suggest the existence of a "fast" channel for BNS mergers, and, thus, a delay-time distribution featuring a peak around time-scales as short as 20 Myr (Perna & Belczynski 2002;Ivanova et al. 2003;Belczynski et al. 2006).These correspond to tight binaries that undergo a third mass transfer episode, and merge while still within star-forming regions in dense environments.These conclusions are corroborated by population study predictions on, for example, r-process element abundances in the Milky Way (Côté et al. 2017) or the redshift distribution of short GRBs (D'Avanzo et al. 2014).The two latter studies suggest a delaytime distribution with a slope − 1, favoring a population of fast mergers, and therefore possibly mergers in dense external media.However, it has been pointed out that the conclusions of population synthesis studies are somewhat sensitive to the assumptions on the physics of the common envelope phase (Dominik et al. 2012) or the distribution of natal kicks (Safarzadeh & Côté 2017).
A second approach is the study of the delay times and kick velocities of Galactic systems.This approach is limited by statistics and by the uncertainty in estimating these from observations.However, finding short delay times or weak natal kicks can imply that a significant fraction of double neutron star mergers should occur in regions where star formation may still be significant, and in turn, the densities are large too.Recently, Beniamini & Piran (2019) have shown that at least 10−20% of Galactic systems are born with delay times of less than 100 Myr between formation and merger.Furthermore, Beniamini & Piran (2016) have shown that the majority of the observed BNSs received relatively weak kicks at birth (v kick 30 km s −1 , see also Tauris et al. 2017).
Another approach is to consider the nature of short GRB host galaxies.On the one hand, these are found to be star-forming two to three times more often than they are found to be elliptical galaxies (Berger 2014).This suggests higher density media for a significant fraction of mergers.This is particularly noteworthy since up to a redshift z 1, that is, where short GRB hosts can be seen, elliptical and star-forming galaxies share roughly equal fractions of the cosmic stellar mass (Bell et al. 2003).This suggests that short GRBs are preferentially found in lower mass galaxies, and thus experience larger external densities on average (Zheng & Ramirez-Ruiz 2007).Also, the observed host galaxy offset distribution has a median value of 1.5 half-light radii, with ∼20% of objects lying outside five half-light radii and ∼20% within one half-light radius (Fong & Berger 2013;Berger 2014).This favors higher density environments for the most centered ∼20% of systems.However, host-galaxy completeness of typical samples is small.Moreover, the offset distribution relies on a correct identification of the host galaxy, and may be grossly overestimating the true offset if, for example, the true host is a fainter, unobserved galaxy of lower mass or higher redshift (e.g., Behroozi et al. 2014).
Insight into short GRBs occurring in dense environments also comes from GRB afterglow observations.On the one hand, Nysewander et al. (2009) have shown that (i) short and long GRBs present a similar correlation between X-ray flux and gamma-ray fluence, (ii) above a gamma-ray fluence threshold of 10 −7 erg cm −2 , optical afterglows are detected in almost all short GRBs and (iii) short and long GRB afterglows have similar radio-to-X-ray flux ratios.These results prompted Nysewander et al. (2009) to suggest that short GRBs have similar or larger external densities to long GRBs, with typical values that may be as large as 1 cm −3 .For a selected sample of short GRB early afterglows, O'Connor et al. ( 2020) have found that less than 16% of events took place at densities smaller than 10 −4 cm −3 , suggesting that few short GRBs occur in very rarefied media.On the other hand, short GRB afterglow catalogs such as Fong et al. (2015) or Berger (2014) do not exhibit a population of high-density afterglows.Similarly, these studies are limited by poor afterglow sampling, parameter degeneracy in photometry fitting and, often, by a lack of the synchrotron self-Compton cooling component in the radiation modeling.In recent years, with the detection of long-lived emission from GRBs with the Large Area Telescope onboard Fermi (Ajello et al. 2018), the synchrotron self-Compton cooling channel has been realized to be an important ingredient of the physical picture.As the Compton parameter affects the position of the cooling frequency, using the cooling break in the X-ray band to estimate the density while disregarding the synchrotron self-Compton effect can particularly bias the result (Beniamini et al. 2015).These caveats may impede a reliable estimation of the circum-burst density and explain this apparent contradiction.Multi-messenger determination of the viewing angle θ v and circum-merger medium density n in the case of GRB170817A.We present 1-σ confidence regions (solid line: median; dashed line: 68% confidence limits) obtained from the GW data assuming the source localization (red), the radio afterglow's properties around its peak (black, see Eq. ( 1)) and very long baseline interferometry imaging measurements (blue).Green triangles show the upper limit on n deduced from the (as yet) undetected kilonova afterglow.The preferred region for θ v and n is highlighted in purple.The text gives details and references.
Finally, an independent approach to short merger binaries comes from r-process abundance studies.The arguments in favor of short merger times, and therefore possibly mergers in dense environments, have recently been summarized in some detail in Hotokezaka et al. (2018a) and Beniamini & Piran (2019).A prevalence of short merger times is implied by (i) observations of r-process enriched stars in ultra-faint dwarf galaxies (Beniamini et al. 2016), (ii) the large scatter of rprocess abundances in extremely metal-poor stars in the Milky Way halo (Argast et al. 2004;Tsujimoto & Shigeyama 2014;Wehmeyer et al. 2015;Vangioni et al. 2016;Beniamini et al. 2018), (iii) the declining rate of deposition of radioactive 244 Pu and 247 Cm on Earth (Hotokezaka et al. 2015;Wallner et al. 2015;Beniamini & Hotokezaka 2020) and (iv) the declining rate of [Eu/Fe] as a function of [Fe/H] observed in Milky Way stars for [Fe/H] −1 (Matteucci et al. 2014;Côté et al. 2016;Komiya & Shigeyama 2016;Hotokezaka et al. 2018a;Simonetti et al. 2019).However, these conclusions rely on knowledge of the rates and r-process yields of BNS mergers, core-collapse and thermonuclear supernovae, all of which are still a matter of debate (see Cowan et al. 2019 andHotokezaka et al. 2018a for reviews respectively on the rprocess in general and on BNS mergers as its astrophysical site).

Determining the apparent fraction of high-density mergers from afterglow observations
We now describe the method we suggest to directly probe the class of high-density mergers.Our method relies on a sample of afterglow counterparts to GW signals from BNS mergers, which would have a sufficient completeness in density above a certain limiting afterglow flux.Population models such as Gottlieb et al. (2019) or Duque et al. (2019) apply criteria based on afterglow flux levels, and thus provide predictions on detectable events.Therefore, applying a flux cut to a sample of detected afterglows ensures that the sample actually represents all the detectable events above the threshold.This in turn allows one to safely use the predictions from population models to compensate for the density-selection effect and infer the intrinsic fraction of highdensity events f HD from the apparent fraction f obs HD , that is, the one observed in the sample.
In this section, we describe how to estimate f obs HD for a sample of afterglow counterparts to BNS mergers.This can be done by inferring the densities of individual events from multi-messenger observations, or directly on the level of the entire sample.

Measuring the viewing angle and density for a single merger event
Combining the GW and electromagnetic (EM) information channels allows one to place individual events quite accurately in the θ v − n plane, as has been done in Fig. 1 for the case of GRB170817A.
First, in Fig. 1 we present the constraints on θ v obtained from the GW data using the information on the event localization from the EM counterpart, as was found by Finstad et al. (2018).These are marked in Fig. 1, and are representative of three-interferometer constraints that can be obtained in the favorable case where the source is pin-pointed thanks to the detection of the kilonova or early afterglow.
Second, we plot the constraint arising from the properties of the light curve of the radio afterglow around its peak.We start from the equation for the 3 GHz afterglow peak flux F p and peak time t p as a function of the jet parameters, in the case where the radio band lies in the [ν m , ν c ] portion of the synchrotron spectrum (Nakar et al. 2002).Combining these two equations in order to write the ratio F p /t 3 p , we eliminate the jet's kinetic energy from the calculation.We then insert the equation relating the afterglow peak "shape factor" η = ∆t/t 2 to the jet opening and viewing angles (Mooley et al. 2018a).Here, t 2 is the onset time of the afterglow's decreasing phase, and ∆t is the afterglow turnover time, counted between the end of the afterglow's increasing phase and the onset of its decreasing phase.This last operation eliminates the jet's opening angle from the calculation, and, finally, we obtain the following relation between observable quantities (left-hand side) and the jet parameters (right-hand side): where D is the luminosity distance to the event, e and B are the usual shock microphysics parameters1 , and α is such that the forward shock Lorentz factor is Γ ∝ t −α .For a jet plowing through a uniform medium, α equals 3/8 for a non-expanding jet, and 1/2 for a jet with sound-speed lateral expansion (Rhoads 1999).The numerical normalization values on the left-hand side of Eq. ( 1) are valid for p = 2.2.We provide these relations in both the expanding and nonexpanding jet hypotheses, which are extreme options regarding the jet lateral dynamics.The actual dynamics should lie somewhere in between, and the discrimination between both can be done on the basis of the post-peak afterglow temporal slope (e.g., Lamb et al. 2018).We note that, in the case of an expanding jet, the θ v − n relation no longer depends on the turnover time, which may prove difficult to measure in the poorly sampled afterglows of marginally detectable events.
Fortunately, the strongest dependencies here are in the measurable quantities t p , F p and D, rather than on the uncertain e and B , allowing us to obtain a thin uncertainty region in the θ v − n plane.This constraint, which requires only data on the afterglow around its peak, is shown in Fig. 1, where we have taken the values of afterglow observables for GRB170817A from Mooley et al. (2018a).Here, the width of the uncertainty region is obtained by propagating the 1-σ uncertainties on t p , F p , D and adding an uncertainty of 0.3 (resp.2) on log e (resp.log B ), deduced from the scatter of its value in GRB jet forward shocks (Beniamini & van der Horst 2017;Nava et al. 2014;Santana et al. 2014).Third, we include the viewing angle constraints from the very long baseline interferometry (VLBI) imagery of the radio remnant.By comparing high-resolution imagery of the remnant to synthetic images based on jet models, Mooley et al. (2018b) and Ghirlanda et al. (2019) were able to constrain the viewing angle to the region shown in blue in Fig. 1.
Finally, we add the constraint that comes from the nondetection of the so-called "kilonova afterglow".This is expected radiation from the forward shock formed by the mildly relativistic material responsible for the kilonova signal on the external medium (Hotokezaka et al. 2018b;Nakar et al. 2018;Kathirgamaraju et al. 2019).Due to the small Lorentz factor and smooth velocity structure of this ejecta, this afterglow component is expected to peak within a decade post-merger in the case of GRB170817A (Kathirgamaraju et al. 2019).The absence of rebrightening of the afterglow, interpreted as the nondetection of the emergence of this component two years after the merger, already constrains the density to n 10 −3 cm −3 (Kathirgamaraju et al. 2019, Fig. 3).
In addition, a detection of the kilonova afterglow would allow for an actual measurement of the density, and not only an upper limit.However, we note that, in both cases, the constraint depends on the assumed value for e in the corresponding shock, which is still uncertain for mildly relativistic shocks.Allowing this parameter to assume values suggested for such shocks ( e 10 −2 , Crumley et al. 2019) by particle-in-cell simulations and by observations of young supernova remnants (Morlino & Caprioli 2012) loosens the bound on n.Therefore, we advise prudence on the use of the kilonova afterglow for measurements of the density.More details on this last point may be found in Sect. 5.
As seen in Fig. 1, the combination of the constraints from the GW, the afterglow light curve and VLBI measurement and the kilonova afterglow leads to θ v ∈ [24, 28] • and log n/cm −3 ∈ [−5, −3] (all 1-σ confidence intervals) for GRB170817A.Disregarding the kilonova afterglow constraint because of the aforementioned uncertainty on e in the corresponding shock, the range of inferred densities becomes log n/cm −3 ∈ [−5, −2].
Such a combination of constraints is only obtained if all the possible multi-messenger observations are made.Using these after a number of events, an estimate of f obs HD can be obtained.It is clear from Fig. 1 that GW and VLBI data crucially narrow down the constraint on θ v .Unfortunately, VLBI remnant imagery will likely become impossible in most cases as the GW horizon increases and we expect its contribution to vanish for most events as of the start of the O3 run (Duque et al. 2019).In the future, this may be compensated for by some improvement in the GW constraint as more interferometers come online, though it will probably be modest (Veitch et  Fig. 2. Corner plot of luminosity distance, 3 GHz afterglow peak flux and time of peak of two populations of mergers: one in density of 10 −3 cm −3 (yellow), and another in 1 cm −3 (blue).Shown here are synthetic populations for radio-GW jointly detectable events as expected from the population model of Duque et al. (2019) for the current O3 run and taking the Very Large Array (VLA) as the limiting radio instrument, with a 3 GHz sensitivity of 15 µJy.
An advantage of this multi-messenger estimation of n is the use of Eq. ( 1), which requires the properties of the radio afterglow around its peak only and thus is applicable even for faint or poorly-sampled afterglows.Also, it can easily be adapted to other bands, such as the optical, provided they lie between ν m and ν c and the afterglow is not outshined by the kilonova.However, Eq. ( 1) is valid only for small densities, when the effects of synchrotron self-absorption in the forward shock are negligible.As illustrated later in Fig. 3, this is no longer the case as soon as n 10−100 cm −3 , depending on the distribution of jet kinetic energies of the population.Nonetheless, from Fig. 3, one expects that at these densities, the X-ray afterglow will be readily accessible and n can be estimated from fully-fledged afterglow fitting, containing more physics than Eq.(1).

Using n -θ v correlations in the sample of merger afterglows
If such follow-up observations are not done and the only available data are GW and afterglow photometry, f obs HD can still be retrieved at the level of the observed sample thanks to important density-dependent correlations in the afterglow peak properties.
In Fig. 2, we plot the distributions of the distance, 3 GHz afterglow peak flux and peak time for two populations of mergers, in high-or low-density media.These are the distributions for the mergers predicted to be detectable jointly in GW and in the radio band by the VLA (with a limiting sensitivity of 15 µJy) for the O3 run and supposing they are placed in media with unique high (n = 1 cm −3 ) or low (n = 10 −3 cm −3 ) densities.These distributions, as all the afterglow populations mentioned throughout this article, were generated as in Duque et al. (2019).That is, progenitor binaries were assumed uniform in space within the Fig. 3. Afterglow recovery fraction in X-ray, optical, and radio bands as function of circum-merger medium density, for a population with energy distribution function deduced from G16 (left) or WP15 (right).We note the effect of synchrotron self-absorption on the recovery fraction in the radio band as of n 10 cm −3 .
GW horizon and isotropic in jet direction, which we suppose is the direction of the system's angular momentum.For each binary, the jet's energy was sampled from an energy distribution function (in Fig. 2, this was deduced from Ghirlanda et al. 2016, see details in Sect.4) and the afterglow radiation was computed using the full synchrotron spectrum, including self-absorption, in the thin shell regime, supposing a top-hat jet with relativistic deceleration dynamics.In the sample, events were deemed GWand radio-detectable by applying thresholds on their GW signalto-noise ratio and afterglow peak flux, respectively.Synchrotron self-Compton effects were ignored in this analysis, as frequencies are always well below ν c .
In particular for t p and F p , the distributions are qualitatively different.The low-density mergers accumulate around the limiting flux, showing that the bulk of the population is undetectable, whereas the high-density mergers present a peak at the mJy level.The combination of these population-level correlations with an adequate statistical treatment of afterglow observations should allow one to estimate f obs HD for the sample.

HD
We now illustrate our method of constraining high-density mergers starting from their apparent fraction f obs HD obtained from multi-messenger follow-up campaigns, as shown in Sect.3.
For the sake of illustration, suppose mergers occur in two different types of media: high-density (n 2 ) and low-density (n 1 ≤ n 2 ).We are interested in inferring from multi-messenger BNS merger observations the intrinsic fractions f HD and f LD = 1− f HD of mergers occurring respectively in media of densities n 2 and n 1 .
For a certain EM band B, let r B (n) denote the "afterglow recovery fraction" at density n, meaning the fraction of mergers occurring at density n to produce a detectable afterglow in the B band.This is provided in Fig. 3 for the X-ray (1 keV), optical (r), and radio (3 GHz) bands, assuming detection limits respectively of 10 −15 erg s −1 cm −2 (50 ks exposure of Chandra in 0.5-8 keV band), magnitude 24 (space telescope routine observation) and 15 µJy (18 ks exposure of VLA in 2-4 GHz band).The plotted r B (n) were determined from populations synthesized for the O3 run as in Fig. 2, but placed in media with densities that are constant within a population but varying from one population to another.Furthermore, they assume two different distributions for the jet kinetic energies: one deduced from the short GRB luminosity function of Wanderman & Piran (2015, hereafter, WP15), the other from that of Ghirlanda et al. (2016, hereafter, G16).In both cases, we deduced the jet kinetic energy distribution from the short GRB luminosity function assuming typical short GRB durations and γ efficiencies of 0.2 s and 20%.Also, we give the multiwavelength afterglow recovery fraction r Mλ (n), which accounts for events detectable in at least one of the three bands.
The luminosity functions found in G16 and WP15 were deduced from distinct GRB data sets and using methods with distinct hypotheses.Among GRB luminosity function studies, they represent two extremes in terms of population steepness, with WP15's luminosity function being much more bottomheavy than that of G16.For our purposes, G16 can be understood as optimistic with regards to afterglow detectability, and WP15 pessimistic.
As mentioned in Sect.3.1, synchrotron self-absorption tends to decrease r 3 GHz (n) as of n 10−100 cm −3 , which appears clearly in Fig. 3.This leads us to consider other bands (and most prominently the X-ray) for the estimation of n in individual events.Therefore, we shall consider r Mλ as the relevant recovery fraction in what follows.
As explained in Sect. 1, because of the strong dependence of the afterglow peak flux to the circum-merger density (F p ∝ n p+1 4 ), we have r(n 1 ) r(n 2 ).Therefore, mergers in high-density media should be over-represented in the observed population with respect to their intrinsic fraction f HD .This establishes a method to effectively constrain the latter following the observation of only a few of these high-density events.
The probability of observing a high-density merger is Furthermore, after observing N afterglow counterparts to GW, the likelihood that a fraction f obs HD will be found to occur in a high-density medium is that of a binomial process with Posterior probability density function of f HD obtained after having observed fraction f obs HD of high-density mergers among ten events, for varying f obs HD .The dashed line shows the current constraint, obtained after the single low-density event GRB170817A.The dotted blue line shows the constraint obtained with f obs HD = 5/10, but ignoring the selection effect, i.e., with r(n 1 ) = r(n 2 ) = 1.Left: assuming the population's jet kinetic energy distribution follows the short GRB luminosity function of G16.Right: same, for that of WP15.success probability p HD and N tries 2 : Finally, since according to Bayes' theorem with no prior information on f HD we have p( f HD | f obs HD , N) ∝ p( f obs HD | f HD , N), a constraint on f HD follows.Given the high sensitivity of the fraction r(n) to the density, we expect these constraints to be tight even with a small number of events.This is clear in Fig. 4, where we have chosen n 1 = 10 −3 cm −3 , n 2 = 1 cm −3 , and we show the constraints that could be obtained from ten events (as expected after three years of an O3-type run, Duque et al. 2019) among which one, three or five are in a highdensity medium.We observe that the constraints do not center around f obs HD and are tighter than if the bias towards high-density events were ignored, as can be seen by comparing the solid blue curves with the dotted blue curves.This illustrates the 'magnifying effect' of the selection by the afterglow.
The slope of the jet energy function is steeper for WP15 than for G16.This implies that, overall, G16 predicts more highenergy events than WP15.This explains why r(n) is systematically larger for G16 than for WP15, at least in the regime where F p ∝ En p+1 4 , that is, before the onset of the self-absorption suppression.This also implies that the rate at which afterglows are recovered by increasing the density is greater for WP15 than for G16.In terms of recovery fraction, this is expressed by saying that the contrast r(n 2 )/r(n 1 ) is larger for WP15 than for G16, which naturally leads to tighter constraints, as is clear from Fig. 4.
In the case where no high-density events are observed, upper limits on the intrinsic fraction f HD can be deduced.This is done in Table 1, where we report the 95%-confidence level upper limits deduced from the observation of N events, all in low-density media.It appears that the observation of only five low-density events (e.g., observing exclusively low-density events during 18 months of an O3-type run, Duque et al. 2019) suffices to constrain f HD , at the 95%-confidence level, to being smaller than 2 Here we denote the binomial coefficient a Table 1.95%-confidence level upper limits on f HD obtained after observing no high-density events among N afterglows, in two short GRB energy function distribution hypotheses.

Discussion
We have presented a method of effectively constraining the class of BNS mergers that occur in high-density media.It is based on the observation of their afterglow counterparts.We will now discuss the limitations, conditions for application and possible extensions of this method.

Method limitations and applicability
A first limitation of the method presented here is the requirement that the sample be density-complete above a certain afterglow flux.In other words, it requires the certitude that all detectable afterglows with fluxes above a limit were effectively detected.
Only in this case can the model-determined recovery fraction r(n) be used to infer f HD from f obs HD .As the observational biases resulting in practical limitations to these detections are discussed in Duque et al. (2019), we do not repeat them here.We only mention that the difficulty in following-up GW events linked to the size of the localization sky-maps should be met by largefield facilities such as the Zwicky Transient Facility (Bellm et al. 2019), and by future high-cadence survey instruments such as the Large Synoptic Survey Telescope (Ivezic et al. 2008).In practice, density-completeness will be difficult to obtain, and an uncertainty on f obs HD must be taken into account in applying this method.
Furthermore, there is a selection bias towards high-density mergers for reasons unrelated to the afterglow flux.For instance, afterglows of mergers occurring in denser media should peak at earlier times, favoring their detection during follow-up, regardless of their flux level.Consequently, the flux-related selection bias we quantified here in r(n) actually underestimates the bias toward high-density events.
Similarly, there is a selection bias towards bright afterglows regardless of the events' circum-merger density.For instance, events closer or brighter in gamma-rays should be better localized by the GW or GRB data, easing their follow-up, regardless of the circum-merger density.These density-unrelated biases towards afterglow detection actually correlate positively with afterglow flux and thus, statistically, with density.Therefore, once again, the bias towards high-density events we consider here is underestimated.
This method is not applicable to the population of cosmological short GRBs for which densities have been estimated, for two main reasons.The first is that the densities claimed for this population are deduced from uncertain fits, as argued in Sect.2, and that only a small fraction of GRBs have a claimed density.Thus, the resulting f obs HD would be quite uncertain.The second is that, for these regular short GRBs, the afterglow detectability depends more on factors that are not density-related, such as (i) the availability of sufficiently rapid follow-up observations and other human factors, or (ii) the quality of the localization of the GRB, which is linked to its prompt properties and not to its afterglow.Also, for regular GRBs, the expected recovery fraction r(n) should be determined through a population model selecting events on joint GRB-afterglow detection, instead of on joint GW-afterglow.
In Sect.3, we mentioned the kilonova afterglow as an alternative means of measuring the merger environment density regardless of the viewing angle, as allowed by the quasi-isotropy of this signal.Nonetheless, we caution against the feasibility and robustness of such a measurement.First of all, as shown by particle-in-cell simulations and the observation of young supernova remnants (Crumley et al. 2019;Morlino & Caprioli 2012), mildly relativistic shocks are expected to be poor electron accelerators, with e up to orders of magnitude lower than in relativistic shocks.Therefore, seeing as the afterglow flux scales with e , the kilonova afterglow should be significantly fainter than the jet's afterglow and unlikely to be detectable in most cases.Furthermore, in the typical case of a low-density medium, this signal is expected to peak up to a decade post-merger, posing some challenge to its detection in follow-up campaigns.Finally, the kilonova afterglow light curve depends on the minimal velocity of the merger ejecta and on its radial velocity structure, both of which are still uncertain for lack of modelling and observation history.Therefore, although the kilonova afterglow signal's quasi-isotropy dismisses the degeneracy between the density and the viewing angle, its use introduces some uncertainty to the measurement, which is thus rendered not robust.

From mergers in dense environments to fast-merging and low-kick binaries
In Sect. 1, we presented this method of determining f HD as a first step towards constraining the population of fast-merging binaries required to explain various astrophysical data, as summed-up in Sect. 2. First of all, an astrophysically interesting constraint on the densities of circum-merger media should be given as a continuous distribution of densities within the population, and not only as a fraction of high-and low-density mergers as we have shown here for simplicity.A continuous (parametric) distribution of densities does not pose any mathematical problems and can be included in this method.
Second, constraining the distribution of merger delay times from that of the merger environment densities is non-trivial, because the medium hosting the merger effectively depends on the locus of the second supernova in the galaxy, on the kick it imparts to the binary system, on the galactic potential, and on the galactic density profile.All of these are uncertain or variable from system to system.As stated in Sect.2, efforts to tackle these effects on the level of population-synthesis models have been done, and are ongoing (O'Connor et al., in prep.).Nonetheless, untangling the effects of all these factors remains difficult.
Third, our method relies on observing the afterglow counterparts to GW inspiral signals, and thus can only inform us on the high-density mergers within the horizon of the GW instruments.However, the fast-merging binary population suggested by the r-process element observations mentioned in Sect. 2 must have formed and enriched their hosts shortly after the peak of cosmic star formation, that is, at z ∼ 2. Thus, this method will remain ineffective with regards to this particular population, as long as we rely on second-generation GW instruments.However, with the prospect of detecting inspiral signals from systems at z 1 with third-generation interferometers such as the Einstein Telescope (Punturo et al. 2014) or Cosmic Explorer (Reitze et al. 2019), the constraining power of this method becomes larger and extends to the redshifts where fast-merging binaries are a matter of debate.In this context, a complete description should require a redshift-varying fraction f HD (z), the addition of which is a straightforward extension of our method.At these redshifts, however, detection of the kilonova may reveal challenging and the localization of the source needed for multiwavelength follow-up should be ensured directly by detection of the afterglow by wide-field X-ray instruments such as Theseus (Amati et al. 2018) or radio survey facilities such as the Square Kilometer Array (Dewdney et al. 2009).

Conclusion
We have described a method of directly probing the binary neutron stars that merge in dense environments, based on the observation of binary neutron star merger afterglows and exploiting the high sensitivity of these to the circum-merger medium density.Its constraining power is large and, since high-density mergers are naturally associated with fast-merging or low-kick binaries, this method is a first step toward a new independent approach to the delay-time and kick velocity distributions.
Fig.1.Multi-messenger determination of the viewing angle θ v and circum-merger medium density n in the case of GRB170817A.We present 1-σ confidence regions (solid line: median; dashed line: 68% confidence limits) obtained from the GW data assuming the source localization (red), the radio afterglow's properties around its peak (black, see Eq. (1)) and very long baseline interferometry imaging measurements (blue).Green triangles show the upper limit on n deduced from the (as yet) undetected kilonova afterglow.The preferred region for θ v and n is highlighted in purple.The text gives details and references.
Fig.4.Posterior probability density function of f HD obtained after having observed fraction f obs HD of high-density mergers among ten events, for varying f obs HD .The dashed line shows the current constraint, obtained after the single low-density event GRB170817A.The dotted blue line shows the constraint obtained with f obs HD = 5/10, but ignoring the selection effect, i.e., with r(n 1 ) = r(n 2 ) = 1.Left: assuming the population's jet kinetic energy distribution follows the short GRB luminosity function of G16.Right: same, for that of WP15.