Evidence for Dynamically Driven Formation of the GW170817 Neutron Star Binary in NGC 4993

The DECam images used in this work were taken as part of the DECam–GW follow-up program between the nights of 2017 August 17 and September 1, using ugrizY filters. We also use public ugrizY DECam data from 2015 June to avoid contamination in the transient region. In addition, we extract YJK data from the VISTA Hemisphere Survey (VHS; McMahon et al. 2013), covering the host galaxy. The images are coadded and registered to a common pixel scale (02636) using SWARP (Bertin et al. 2002) with 3.5σ clipping to remove cosmic-ray artifacts. An RGB coadded image of the galaxy is presented in Figure 1. We build a χ²-detection image from the r-, i-, and z-band data and run SExtractor (Bertin & Arnouts 1996) in dual mode on the coadded images without performing template subtraction.

Zoom In Zoom Out Reset image size

The photometry is corrected for galactic extinction. In order to compare the galaxy properties to a broader sample, we also use DES data from the first year of observations (Y1; Drlica-Wagner et al. 2017). We use MAG_AUTO magnitudes unless otherwise stated.

NGC 4993 was also observed during Hubble Space Telescope (HST) Cycle 24 (PropID 14840, PI: Bellini) using ACS in F606W. The data were publicly released in 2017 April and were accessed via the Hubble Source Catalog (HSC; Whitmore et al. 2016).

The 6dF Galaxy Survey (Jones et al. 2004) final release (Jones et al. 2009) includes an optical spectrum of the center of NGC 4993 with an estimated redshift (z = 0.009680 ± 0.000150).

Spectra of 14 galaxies with ${v}_{\mathrm{helio}}\sim 3000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and within a 1° radius of NGC 4993 were obtained in one target of opportunity exposure of the AAOmega spectrograph at the Anglo-Australian Telescope (AAT) on 2017 August 27. Of those, 10 spectral fits passed quality cuts. All the spectra used here are centered on their galaxy nucleus with a 2'' aperture.

We begin our study of NGC4993 with an analysis of its morphological properties, employing the CAS non-parametric light quantification (Conselice 2003) and parametric Sérsic light profile fitting with GALFIT (Peng et al. 2010). Both methods utilize a mask to exclude other sources in the image and the location of the kilonova event. The CAS system is able to pick out the salient features of galaxy morphology, allowing galaxy types to be assigned and identifying objects that are likely to have undergone a recent major merger (see Conselice 2003 for details). Meanwhile, fitting the light profile additionally provides us with an alternative estimate of the total magnitude and can reveal more subtle aspects of galaxy morphology within the residuals of the model-subtracted image.

GALFIT is run on the DES and VHS images in two ways: band-by-band and simultaneously across all bands using a modified version, GALFIT-m (Vika et al. 2013). In the second case the Sérsic fitting parameters are allowed to vary with wavelength as a second-order polynomial. We extract the point-spread function (PSF) model required by GALFIT from the coadded images with PSFEx (Bertin 2011) and initialize the fitting parameters based on measurements of the galaxy from sextractor. All parameters are left free without constraints, except for the central position in the single-band fits. This is allowed to vary by only ±1 pixel as it is already well constrained by SExtractor.

In order to assess the stability of GALFIT and obtain an estimate of the uncertainties on the measurements, each single-band run is performed 10,560 times, varying the inputs around their nominal values. We take the median as our final measurement and the standard deviation as the uncertainty.

Following the definitions given in Conselice (2003), we find concentration C = 3.348 ± 0.035, asymmetry A = 0.04 ± 0.01, and clumpiness S = 0.05 ± 0.05. These values are typical for an early-type galaxy. In Figure 2, we compare these values to field galaxies of similar masses (within 0.2 dex of NGC 4993) and redshifts (z < 0.2) from the GAMA survey and to a sample of sGRB hosts (C. Conselice et al. 2017, in preparation) taken in F814W imaging from HST. NGC 4993 stands out as peculiar with respect to other GRB hosts: such objects tend to lie on the more highly asymmetric side of late-type galaxies.

Zoom In Zoom Out Reset image size

The results from the single-Sérsic fit across all bands are summarized in Table 1 (the band-by-band fits give broadly consistent results). We find an increase in Sérsic index toward redder bands and a rotation in the position angle. This rotation of bluer versus redder bands suggests there could be two superimposed stellar populations with differing orientations. This may have arisen during the course of the galaxy's secular evolution but could also be caused by a minor galaxy merger, as indicated by the presence of shells.

Table 1.  Outputs from GALFIT Parametric Sérsic Fits Performed on the ugrizY DECam Coadded Images and YJKs VHS Data

Filter	MAG_AUTO	Mag	re	n		θ
u	14.24	14.15	61.8	3.2	0.15	−13.9
g	12.95	12.80	62.5	3.4	0.15	−12.8
r	12.08	11.90	63.5	3.7	0.16	−11.2
i	11.65	11.45	64.4	4.0	0.16	−9.9
z	11.34	11.13	65.3	4.3	0.16	−8.4
Y	11.13	10.96	65.7	4.4	0.16	−7.7
YVHS	11.27	11.00	65.9	4.5	0.16	−7.5
J	11.00	10.77	67.3	5.0	0.17	−5.2
Ks	11.08	10.68	72.9	6.7	0.19	+3.5
		±5 × 10−4	±0.07	±3 × 10−3	±4 × 10−5	±5 × 10−3

Note. The fit was joint across bands, allowing the effective radius, r_e (in pixels), Sérsic index, n, ellipticity,  = 1 − b/a, and position angle, θ, to vary with wavelength. One pixel corresponds to 02636. The final row lists indicative errors based on the single-band analysis.

Download table as:  ASCIITypeset image

The middle panel of Figure 1 shows DECam r-band residuals from GALFIT and the position of the transient. At least four shell structures are clearly visible. The surface brightness radial profile from the residual image shows an excess at the shell positions of $\sim 25\,\mathrm{mag}\,{\mathrm{arcsec}}^{-2}$. Closer inspection with HST data (right panel of Figure 1) reveals a possible further broad inner shell, on which the transient seems to lie, and obvious dust lanes (visible also as a negative residual in the DECam version). In summary, these results provide compelling evidence for a recent minor galaxy merger in NCG 4993, and the location of the kilonova event with respect to the shells leads us naturally to ask whether there is a causal connection between the two, for instance via dynamical interaction.

The r-band absolute magnitude from a 4 arcsec² region around the transient location in the galaxy-subtracted template image is −10.65. This luminosity implies a rather high stellar density in the locale of the BNS coalescence, implying that dynamical interactions between stars are more probable in this region compared to typical galaxy stellar densities.

From Figure 2 we see that clear major galaxy mergers are unusual among sGRB hosts. Furthermore, the other sGRBs are at cosmological distances and thus are mostly undergoing extensive galaxy formation through star formation or merging. If the hosts have to be related by some common features, this is an indication that NGC 4993 has undergone some merging activity, but a minor merger such that the bulk morphology is still elliptical. We thus explore the possibility that the kilonova was a result of a recent galaxy merger in NGC 4993.

We use pPXF (Cappellari & Emsellem 2004; Cappellari 2017) for the spectral fitting. It enables extraction of the stellar kinematics and stellar population from absorption line spectra of galaxies, using a maximum penalized likelihood approach. We use the Miles stellar libraries and fit over the wavelengths 4000–7409 Å, excluding the range 5500–5600 Å of the 6dF spectrum, where a strong sky line contaminates the flux.

We use LePhare (Arnouts et al. 1999; Ilbert et al. 2006) for the broadband spectral energy distribution (SED) fitting. We add a 0.05 systematic uncertainty in quadrature to the magnitudes. The simple stellar population (SSP) templates used are Bruzual & Charlot (2003), with two metallicities (${Z}_{\odot }$ and $2.5\,{Z}_{\odot }$), a Chabrier (2003) initial mass function (IMF), and a Milky Way (Allen 1976) extinction law. The SFH chosen is lognormal:

$$\begin{eqnarray}&&{\rm{\Psi }}(t,{t}_{0},\tau )=\displaystyle \frac{1}{t\sqrt{2\pi {\tau }^{2}}}{e}^{-\tfrac{{(\mathrm{ln}t-\mathrm{ln}{t}_{0})}^{2}}{2{\tau }^{2}}},\end{eqnarray} \tag{ 1 }$$

as it is the most representative family of models with only two parameters (Gladders et al. 2013). Here, t₀ and τ are the half-mass time and width.

Motivated by our morphological analysis, we allow for an additional burst of recent SF. This is modeled as a Gaussian centered at t_burst with a width of 10 Myr and peaking at a fraction 0.4–0.1 of the peak of the lognormal SFH (as no evidence for strong late SF is found).

The same templates are used to perform spatially resolved SED fitting across DES+VHS coadded images within 10 × 10 pixels, including the galaxy dust extinction. The other sources in the field are masked out using the segmentation map output by sextractor.

Figure 3 shows the best-fit model of the 6dF spectrum, which results in a reduced χ² of 1.22. An analysis of the mass fraction in age shows that part of the core galaxy stellar population has a supersolar metallicity, but the weighted mean value $\langle [M/H]\rangle =-0.012\pm 0.010$ is marginally consistent with solar metallicity. The mean age is 11.298 ± 0.054 Gyr, and the mass-to-light ratio is 5.23 ± 0.15 in the r band.

Zoom In Zoom Out Reset image size

The stellar model fit reveals the existence of weak ionized gas emission lines. However, the line ratios from the fit suggests they are produced by a harder ionizing source than star formation, formally lying in the AGN region of the Baldwin, Phillips, and Telervich (BPT; Baldwin et al. 1981) diagram. Blanchard et al. (2017) argue that there is a weak AGN present in the core of the galaxy on the basis of radio and X-ray emission, and so we conclude that there is no evidence of recent star formation from the 6dF spectrum, irrespective of the highly uncertain O iii/H_β ratio. A comparison of galaxies in the group using AAT spectra and classification by Kewley et al. (2006) is shown in Figure 3.

Given the evidence of dust presence in the HST study, we estimate the dust content using the Balmer decrement (Berman 1936) observed from the spectrum. The reddening is E(B − V) = 0.12 ± 0.50 in the case of I(H_α)/I(H_β) = 3.1, which is expected in the case of AGN activity. We therefore restrict our dust models to have reddening values 0.1, 0.2, 0.3, 0.4, and 0.5 in the photometric fits.

The photometric best-fitting template has a solar metallicity, a quickly declining lognormal SFH with t₀ = 3 Gyr, and τ = 0.1. A low reddening E(B − V) = 0.1 is preferred, and the stellar mass is (2.95 ± 0.65) × 10¹⁰ M_⊙. The inclusion of a late SFH burst is disfavored by the fitting apart from intermediate apertures.

Previous work found that the presence of dust lanes may bias the galaxy stellar mass from unresolved galaxy SED fits to lower values (Sorba & Sawicki 2015). The total stellar mass from fitting over the SExtractor segmentation map of NGC 4993 is (3.8 ± 0.20) × 10¹⁰ M_⊙, more than 1σ higher than the unresolved SED fitting. The specific SFR (sSFR) map from our pixel SED fits is shown in Figure 4, where the shell structure is clearly visible, suggesting that the sSFR is slightly more accentuated in the stellar halo compared to the inner parts. Younger ages (by ∼2 Gyr) are also preferred in the outer regions, though we still do not find evidence for a star formation burst at late times and explain our results by the stripping of stellar populations from the lower-mass galaxy in a minor dry merger. A dust model with E(B − V) = 0.1 is preferred in the inner few kiloparsecs, while E(B − V) = 0 is found outside. Despite the presence of dust lanes, an analysis of the HST photometry and a comparison with extinction models suggests that the effect of dust is not extreme, with reddening values that are consistent with 0.1 in the core. We therefore believe that the dust obscuration does not play a significant role in our SFR estimates.

Zoom In Zoom Out Reset image size

In Figure 4, we show a color–magnitude diagram for all the pixels within the field of view of the DECam data near the galaxy. The image has been cleaned of stars and other contamination, thus all points come from the galaxy itself. The position of the GW source, 106 offset from the center, is the cyan colored pixels, while the center of the galaxy is shown by the red points. This galaxy is well represented by a pixel "main sequence" that is bluer at fainter levels, which is typical of early-type galaxy color gradients (e.g., Lanyon-Foster et al. 2007). We conclude that there is no significant difference between the transient position and other outer light, although it is bluer than the core region. This further supports the scenario in which the BNS formation is not related to some particular recent star formation event in this region.

In the most accepted shell formation scenarios, the shells are stellar debris coming from the less massive, stripped galaxy, and the arcs form at the apocenter of the orbits of the infalling material (Quinn 1984).

Based on our results, we believe that NGC 4993 experienced a dry minor galaxy merger with still visible signs. The shells are expected to be washed out within a time that depends on the velocity dispersion at their position. We estimate the shell survival time in two ways: based on the velocity dispersion of the galaxy, as well as the velocity dispersion of the shell itself. From the 6dF spectrum the line-of-sight central velocity dispersion is ${\sigma }_{v}=(160.0\pm 9.1)\,\mathrm{km}\,{{\rm{s}}}^{-1}$. We estimate its value at the position of the transient. The velocity dispersion of early-type galaxies drops from its central peak value at larger radii, and observations show that the maximum drop to the outer parts of ellipticals near the effective radius is ∼40% of the central value (Emsellem et al. 2004). Based on the distance of the shell from the center, R ≈ 4 kpc, we estimate that the dynamical time at this radius is t_dyn ≡ R/σ_v ≈ 60 Myr (the line-of-sight velocity is relevant here, given the shell's geometry, but, e.g., if we assume a 3D isotropic velocity dispersion it would reduce the dynamical time by $\sqrt{3}$).

So far, we have no measurement of the shell's velocity dispersion, but estimates from the literature suggest for similar shells in other galaxies ${\sigma }_{v}\approx 20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Quinn 1984). This would give a dynamical timescale of t_dyn ≈ 192 Myr. Detailed simulations of shells in other galaxies suggest that survival time could be even larger that 1 Gyr, depending on the assumed scenario (Pop et al. 2017).

The survival time of the shell could be used as an upper limit for the time the minor merger took place, i.e., t_dyn ≥ t_merg, so we estimate t_merg ≲ 200 Myr.

If the BNS was formed as such in a shell, then we would have expected to see evidence for recent star formation, but we find no indication of this. In the absence of star formation, it is plausible that the BNS coalescence was triggered by a dynamical process, e.g., NS–NS capture or the destabilization of a pre-existing wide-separation binary. These processes will be quite sensitive to the stellar density, which, given the Sérsic index and the luminosity from the residual image found in Section 3.2, is high in the center of NGC 4993 and around the transient position. If this dynamical hypothesis is true, then the delay time Δt_NSM between the BNS formation and coalescence is ≲200 Myr. On the other hand, Blanchard et al. (2017) find a median delay time of ${11}_{-1.4}^{+0.7}\mathrm{Gyr}$ under the assumption that the binary was formed through secular SF.

If the binary formation is related to dynamical processes in galaxy merging, as we are investigating here, then this is most likely to happen in galaxy groups and low-mass clusters. According to the 2MASS catalog (Tully 2015), NGC 4993 resides in a group of which we analyze the remaining seven galaxies. A spectral analysis shows that NGC 4993 is not the only galaxy showing AGN activity (see Figure 3), but it is peculiar in terms of age, metallicity, and mass-to-light ratio. It shows an older stellar population (the mean age of the other 13 galaxies is Log (Age) = 9.56 ± 0.17), lower metallicity (mean: M/H = −0.31 ± 0.11), and higher M/L_r (mean: 2.41 ± 0.45) than the average. The group has a projected virial radius of R_vir = 0.36 Mpc and a line-of-sight velocity dispersion ${\sigma }_{v}=143\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Tully 2015). The crossing time is therefore ${t}_{\mathrm{cr}}\sim {R}_{v}/(\sqrt{2.5}{\sigma }_{v})\sim 1.6\,\mathrm{Gyr}$.

If galaxy mergers are correlated to BNS coalescence, future GW studies could possibly concentrate on galaxy groups (but note that these are crowded regions and therefore matching candidates to a host could be difficult). In order to have precise measurements of H₀, one needs to identify the host galaxy redshift clearly. When the match is clear, the properties of the type of host galaxy found could help future studies to select the right host galaxy or create galaxy catalogs of likely hosts for GW–EM follow-up and untriggered kilonova searches (Doctor et al. 2017). In fact, large photometric surveys such as DES, LSST, or WFIRST are expected to observe kilonova events at redshifts beyond the sensitivity of GW experiments, where the angular separation between galaxies decreases (Scolnic et al. 2017).

We derive a constraint on neutron stars merging rate at time t by using

$$\begin{eqnarray}&&{R}_{\mathrm{NSM}}(t)=\alpha {R}_{\mathrm{NS}}(t^{\prime} ),\end{eqnarray} \tag{ 2 }$$

where α is the fraction of neutron stars that are in binaries, t'= t − Δt_NSM, and the fraction of mass of formed stars that are NSs is

$$\begin{eqnarray}&&{R}_{\mathrm{NS}}(t^{\prime} )=\int {{dM}}_{\star }{\rm{\Phi }}({M}_{\star }){\rm{\Psi }}({t}_{\star }){{\rm{\Theta }}}_{\mathrm{NS}}({M}_{\star }),\end{eqnarray} \tag{ 3 }$$

with ${\rm{\Phi }}({M}_{\star })$ being the IMF, Ψ(t_⋆) is our best-fit SFH, and Θ_NS(M_⋆) is 1 for star mass ranges of 8 M_⊙ < M < 20 M_⊙, and zero otherwise. We drop the metallicity dependence in Θ_NS because we only consider a solar metallicity for the galaxy, as a result of our spectroscopic fit. t_⋆ is the time when the progenitor of the NS was formed, therefore satisfying t' = t_⋆ + t_life, with t_life being the lifetime of the progenitor before becoming an NS. We assume a t_life = 0.02 Gyr, but our calculation is insensitive to this choice as the typical lifetime of these massive stars (∼0.01–0.03 Gyr) is much shorter than the timescale over which the SFH found for NGC 4993 is changing at late times. We assume a Chabrier IMF, but this choice is not relevant, as we are only exploring the high-mass end of the IMF. Assuming α = 0.002 and the distribution of Δt_NSM from Vangioni et al. (2016, their Figure 3 for solar metallicity), and our best-fit SFH from Equation (1) with t₀ = 3 Gyr and τ = 0.3, we get an NS formation rate of ${R}_{\mathrm{NS}}^{\mathrm{gal}}={3.6}_{-3.6}^{+28}\times {10}^{-5}\,{\mathrm{yr}}^{-1}$ and a BNS merger rate of ${R}_{\mathrm{NSM}}^{\mathrm{gal}}={5.7}_{-3.3}^{+0.57}\times {10}^{-6}\,{\mathrm{yr}}^{-1}$ for the whole galaxy. Errors reflect the uncertainty on the SFH, which dominates our errors: they represent the two central quartiles of the rates distribution computed with the SFHs of the pixel SED fitting over the galaxy.

Given the sensitivity of the BNS merger event rate to the recent SFR of a galaxy, it is somewhat surprising that GW170817 occurred in an old, early-type galaxy. We therefore ask what is the probability of observing such an event in any early-type galaxy within the LIGO-detectable volume. To make this estimate, we integrate the stellar mass function of early-type galaxies from Weigel et al. (2016) and scale the per-solar-mass rate from Equation (2) to the mass contained within the LIGO-detectable volume (radius 80 Mpc). We find ${R}_{\mathrm{NSM}}^{\mathrm{early}}={23}_{-14}^{+2}\,{\mathrm{yr}}^{-1}{\mathrm{Gpc}}^{-3}$ resulting in ${0.038}_{-0.022}^{+0.004}$ expected events. This calculation assumes that the SFH of NGC 4993 is representative of local early-type galaxies. In fact, much of the mass will be contained in more massive, and on average older and less star-forming, galaxies. We contrast this with a similar calculation for all galaxy types, using the cosmic SFR density from Gladders et al. (2013), finding ${R}_{\mathrm{NSM}}^{\mathrm{all}}\approx 270\,{\mathrm{yr}}^{-1}{\mathrm{Gpc}}^{-3}$ and ∼0.5 expected events.

This result shows that it is unlikely that we observed one such BNS merger with LIGO over the combined nine months of operations in an early-type galaxy. The assumptions in the calculation include the fraction of NS that form in binaries (α = 0.002) and the delay time distribution, both coming from binary star models (where the progenitors of the BNS were already a bound system) and satisfying Milky Way constraints. If the BNS formation mechanism is via dynamical interaction, our result could point to a higher value of α or a shorter Δt_NSM for systems that recently underwent a galaxy merger, more so for those that have high stellar density (such as early-type galaxies). It is therefore of interest to know the fraction of galaxies similar to NGC 4993 that show similar signs of a galaxy merger in the form of visible shells. We select galaxies from the first year of DES data with size, surface brightness, and Sérsic index within 10% of the best-fit values for NGC 4993. We find 1100 such galaxies, and we visually inspect them to identify shell galaxies. Only 15% of these objects display shells, and so NGC4993 is unusual among early-type galaxies. This is far from conclusive evidence for a merger origin of BNS events. However, the coincidence of evidence for a recent merger in a galaxy for which a BNS event was otherwise improbable is compelling.