Probing high-density binary neutron star mergers with afterglow counterparts

The up-to-now only binary neutron star merger gravitational wave event with detected electromagnetic counterparts, GRB170817A, occurred in a rarefied medium with a density smaller than $10^{-3}-10^{-2}~{\rm cm}^{-3}$. As neutron star binaries are imparted kicks 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 a class of faster-merging binaries, which would coalesce in denser environments. Nonetheless, astronomical data is 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 sharp sensitivity of these to the circum-merger medium density. This method is based on a yet-to-be-collected sample of merger afterglows. 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 3$^{\rm rd}$-generation gravitational wave detectors, this method's potential will be fully realized as it will allow to probe mergers having occurred soon after the peak of cosmic star formation.


Introduction
The distribution of delay-times between the formation and merger of binary neutron stars (BNS) remains an open question. A binary's merger time depends strongly on its post-formation separation and eccentricity, both of which are poorly constrained in the general picture or may be variable from a system to another. Natal kicks imparted to systems by the second supernova allow them to migrate from their formation sites within dense star-forming regions to their more rarefied merger loci. However, for the binaries with the shortest delay-times, we expect the circum-merger medium to be significantly denser, owing to their shorter migration times.
Modelling of the afterglows of short gamma-ray bursts (GRB) is one probe of the environments of BNS mergers. 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).
The up-to-now only BNS merger gravitational wave (GW) event with afterglow counterpart, GRB170817A, occurred in a medium typical of early-type galaxies. Indeed, the circum-merger density was n 10 −3 cm −3 (e.g., Troja et al. 2019, and see Sec. 3.1 below). Throughout this work, we will use 'dense' to refer to media with densities significantly larger than this, i.e., with n 1 cm −3 .
In recent work, Gottlieb et al. (2019) and Duque et al. (2019) have studied the afterglows expected as counterparts E-mail: duque@iap.fr to GW signals from BNS mergers in present and future observing runs. Starting from a population model motivated by short GRB knowledge, 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 nonthermal electron population accelerated at the jet's forward shock front. Thus, should there be mergers in high-density environments, these will be over-represented in the afterglow population with respect to their actual number. In other words, the radio afterglow acts as a magnifying glass for these higher-density mergers.
It follows from this effect that, with a statistical fluxlimited 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 and, comparing to population predictions in order to compensate for the density-selection effect, constrain the intrinsic fraction of mergers in high-density media. This is the principle of the method we suggest.
It is hypothesized that short delay-time mergers occur in high-density media. Conversely, the progenitor of a merger observed to occur in a high-density medium is likely to have had a short delay-time. However, the exact link between the distribution of delay-times and that of circum-merger densi- ties is not clear, in particular because of the aforementioned uncertainty on the natal kicks and initial binary separation and eccentricity. Nonetheless, the method we suggest here is a first step to revealing the delay-time distribution of BNSs from their merger afterglows, in particular for the fastest-merging binaries.
This article is organized as follows. In Sec. 2, we motivate this new method of constraining high-density mergers by recalling some related observational and theoretical knowledge. In Sec. 3, we show that, for future BNS merger afterglows, multi-messenger observations will allow to estimate the circum-merger density and, thus, to determine the apparent fraction of high-density mergers quite accurately. In Sec. 4, we describe how these observations provide sharp 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 Sec. 5 we discuss the limitations of this method, and conclude in Sec. 6.

Elements of indirect evidence on high-density BNS mergers
Over the years, the question of high-density BNS mergers has been approached by various methods. However, as we show here, current data is inconclusive regarding the possibility 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 about time-scales as short as 20 Myr (Perna & Belczynski 2002;Belczynski et al. 2006). These correspond to tight binaries which 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, e.g., r-process element abundances in the Milky Way  or the redshift distribution of short GRBs (D'Avanzo et al. 2014). The two latter studies suggest a delay-time distribution with a slope −1, favoring a population of fast-mergers. However, it has been pointed out that population synthesis studies' conclusions 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 angle is the study of the delay-times and natal kicks of Galactic BNSs. This approach is limited by statistics and by the uncertainty in estimating the latter 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 the star formation may still be large, 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,  have shown that the majority of the observed BNSs have received relatively weak kicks at birth (v kick 30 km/s, see also Tauris et al. 2017).
Another approach is from the nature of short GRB host galaxies. On the one hand, these are found to be starforming 2 to 3 times more often than they are found to be elliptical galaxies (Berger 2014). This suggests higherdensity 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 (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 5 half-light radii and ∼ 20% within 1 half-light radius (Fong & Berger 2013;Berger 2014). This favors higher-density environments for the most centered ∼ 20% of systems. However, the 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 on 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 have 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 . This conclusion, however, is limited by short GRB afterglow-completeness. 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 modelling. These caveats may disallow a reliable estimation of the circum-burst density and explain this apparent contradiction. In recent years, with the detection of long-lived emission from GRBs with the Fermi-LAT (Ajello et al. 2018), the synchrotron self-Compton cooling channel has been realized to be an important ingredient of the physical picture (e.g., Beniamini et al. 2015) and could particularly affect circum-burst density estimates.
Finally, an independent approach to short-merger binaries (and thus large densities) comes from r-process abundance studies. The arguments in favor of short merger times have recently been summarized in some detail in Hotokezaka et al. . 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 r-process in general and on BNS mergers as its astrophysical site).
Furthermore, there is some theoretical support for the existence of fast-merging BNSs. Indeed, these can result from mechanisms such as (i) an efficient common envelope phase, reducing initial separation (e.g., Dominik et al. 2012) or (ii) a favorable kick at the second supernova, giving high eccentricity and thus rapid merger (e.g., Kalogera 1996). This last mechanism seems unlikely in the case of the Milky Way, where there appears to be no correlation between inferred delay-times and eccentricities in the observed systems (Beniamini & Piran 2019). 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), VLBI measurements (blue). Green triangles show the upper limit on n deduced from the yet-undetected kilonova afterglow. The preferred region for θv and n is highlighted in purple. See text for details and references.

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 has a sufficient density-completeness 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 insures 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 high-density events f HD from the apparent fraction f obs HD , i.e., that observed in the sample.
In this section, we will 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 to place individual events quite accurately in the θ v − n plane, as is done in Fig. 1 for the case of GRB170817A. First, we present in Fig. 1 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 red in Fig. 1, and are representative of three-interferometer constraints which 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 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 , in the case where the radio band lies in the [ν m , ν c ] portion of the synchrotron spectrum (Nakar et al. 2002). We then insert the equation relating the afterglow peak 'shape factor' η = ∆t/t 2 to the jet opening and viewing angles (Mooley et al. 2018b). 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. 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 parameters, 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 values on the left-hand side of Eq. 1 are valid for p = 2.2. We provide these relations in both the expanding and non-expanding 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 expanding jet case, the θ v −n relation no longer depends on the turnover time, which may reveal 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 to obtain a thin uncertainty region in the θ v − n plane. This constraint is shown in black in Fig. 1, where we have taken the values for GRB170817A from Mooley et al. (2018b) and the width of the uncertainty region is the sum of the 1-σ uncertainties on t p , F p , D and a typical uncertainty of 1 on log e,B .
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. (2018a) 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 coming 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 in the case of GRB170817A (Kathirgamaraju et al. 2019). The absence of rebrightening of the afterglow, interpreted as the non-detection of the emergence of this component two years after the merger, constrains the density to being n 10 −3 cm −3 (Kathirgamaraju et al. 2019, Fig. 3), as indicated in green in Fig. 1. Note, however, that this constraint is given for e = 0.1 in the corresponding shock, and allowing this parameter to assume values suggested for mildly-relativistic shocks ( e 10 −2 , Crumley et al. 2019) loosens the bound on n. Note also that in the case of a detection of the kilonova afterglow, the constraint takes the form of a band in the θ v − n plane, rather than an upper limit.
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 by some improvement of the GW constraint as more interferometers come online, though it should be modest (Veitch et al. 2012;Ghosh et al. 2016).
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 . 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 full-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 by the VLA (with a limiting sensitivity of 15 µJy) for the O3 run, according to the fiducial model of Duque et al. (2019), and placed in media with unique high (n = 1 cm −3 ) or low (n = 10 −3 cm −3 ) densities. Note that in this article, afterglow populations and derived quantities were obtained with full synchrotron simulations, including synchrotron self-absorption (Panaitescu & Kumar 2000). Synchrotron self-Compton effects can be 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 to estimate f obs HD for the sample.

Constraining high-density mergers with f obs 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 Sec. 3.
For the sake of illustration, suppose mergers occur in two different types of media: a high density n 2 and a low density n 1 ≤ n 2 . We are interested in inferring from multimessenger 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, i.e., the fraction of mergers occurring at density n to produce a detectable afterglow in the B band. It 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/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). They were determined from populations synthesized for the O3 run as in Fig. 2, but placed in media with densities unique within a population but varying from a population to another, and assuming two different distributions for the jet kinetic energies: one deduced from the short GRB luminosity function of Wanderman & Piran (2015) (hereafter, WP15), the other by that of Ghirlanda et al. (2016) (hereafter, G16), in both cases assuming typical short GRB durations and γ-efficiencies of 0.2 s and 20%. Also, we give the multi-wavelength afterglow recovery fraction r Mλ (n), which accounts for events detectable in at least one of the three bands.
As mentioned in Sec. 3.1, synchrotron self-absorption tends to decrease r 3GHz (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 Sec. 1, because of the strong dependence of the afterglow peak flux to the circum-merger den-sity (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 makes up for a method to strongly constrain the latter from having observed 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 were found to occur in a high-density medium is that of a binomial process with success probability p HD and N tries: Finally, as, 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 indeed clear in Fig. 4, where we have chosen n 1 = 10 −3 cm −3 , n 2 = 1 cm −3 , and we show the constraints which could be obtained from 10 events (as expected after 3 years of a O3-type run, Duque et al. 2019) among which 1, 3 or 5 in a high-density 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 high-energy 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 , i.e., before the onset of the self-absorption suppression. This also implies that the rate A&A proofs: manuscript no. main Posterior probability density function on fHD obtained after having observed fraction f obs HD of high-density mergers among 10 events, for varying f obs HD . In dashed line, the current constraint, obtained after the single low-density event GRB170817A. In dotted blue line, the constraint obtained with f obs HD = 5/10, but ignoring the selection effect, i.e., with r(n1) = r(n2) = 1. Left: assuming the population's jet kinetic energy distribution follows the short GRB luminosity function of G16. Right: same, for that of WP15. 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 Tab. 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 5 low-density events (e.g., observing exclusively lowdensity events during 18 months of a O3-type run, Duque et al. 2019) suffices to constrain f HD , at the 95%-confidence level, to being smaller than 18.5% (resp. 9.4%), assuming the short GRB luminosity function of G16 (resp. WP15).

Method limitations and applicability
A first limitation of the method presented here is the requirement of density-completeness of the sample above a certain afterglow flux. In other words, it requires the insurance that all detectable afterglows with fluxes above a limit were effectively detected. Only in this case can the modeldetermined 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 will not recall them. We only mention here that the difficulty to follow-up GW events linked to the size of the localization sky-maps should be met by large-field 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 highdensity mergers for reasons unrelated to the afterglow flux. For instance, afterglows of mergers occurring in denser media will peak at earlier times, favoring their detection during follow-up, regardless of their flux level. Consequently, the flux-related selection bias we consider here, and which is encapsulated in r(n), in fact, underestimates the bias towards 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 will 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, of course, 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 Sec. 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 which 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.

From high-density mergers to fast-merging binaries
In Sec. 1, we present this method of determining f HD as a first step towards constraining the population of fastmerging binaries required to explain various astrophysical data, as summed-up in Sec. 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 delaytimes from that of the merger densities is non trivial, because the actual medium hosting the merger depends on the locus of the second supernova, 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 commented in Sec. 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 Sec. 2 must have formed and enriched their hosts by shortly after the peak of cosmic star formation, i.e., at z ∼ 2. Thus, this method will remain ineffective with regards to this particular population, as long as we rely on 2 nd -generation GW instruments. However, with the perspective of detecting inspiral signals from systems at z 1 with 3 rd -generation interferometers such as the 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 the 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.

Conclusion
We have described a method of directly probing the binary neutron stars merging in high-density 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, as high-density mergers are naturally associated with fast-merging binaries, this method is a first step towards a new independent approach to the delay-time distribution.