This site uses cookies. By continuing to use this site you agree to our use of cookies. To find out more, see our Privacy and Cookies policy.

The following article is Open access

Masses of White Dwarf Binary Companions to Type Ia Supernovae Measured from Runaway Velocities

, , , and

Published 2021 December 28 © 2021. The Author(s). Published by the American Astronomical Society.
, , Citation Evan B. Bauer et al 2021 ApJL 923 L34 DOI 10.3847/2041-8213/ac432d

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

2041-8205/923/2/L34

Abstract

The recently proposed "dynamically driven double-degenerate double-detonation" (D6) scenario posits that Type Ia supernovae (SNe) may occur during dynamically unstable mass transfer between two white dwarfs (WDs) in a binary. This scenario predicts that the donor WD may then survive the explosion and be released as a hypervelocity runaway, opening up the exciting possibility of identifying remnant stars from D6 SNe and using them to study the physics of detonations that produce Type Ia SNe. Three candidate D6 runaway objects have been identified in Gaia data. The observable runaway velocity of these remnant objects represents their orbital speed at the time of SN detonation. The orbital dynamics and Roche lobe geometry required in the D6 scenario place specific constraints on the radius and mass of the donor WD that becomes the hypervelocity runaway. In this Letter, we calculate the radii required for D6 donor WDs as a function of the runaway velocity. Using mass–radius relations for WDs, we then constrain the masses of the donor stars as well. With measured velocities for each of the three D6 candidate objects based on Gaia EDR3, this work provides a new probe of the masses and mass ratios in WD binary systems that produce SN detonations and hypervelocity runaways.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

The "dynamically driven double-degenerate double-detonation" (D6) scenario for producing Type Ia supernovae describes a double white dwarf (WD+WD) binary system in which dynamically unstable mass transfer ignites a supernova (SN) detonation that fully disrupts the accreting WD, possibly releasing the donor WD from the binary as a hypervelocity runaway star at approximately its former orbital velocity (Shen et al. 2018a, 2018b). Simulations have identified this scenario as a promising mechanism for igniting detonations in carbon–oxygen (C/O) WDs below the Chandrasekhar mass (Guillochon et al. 2010; Dan et al. 2011, 2012; Pakmor et al. 2012, 2013; Raskin et al. 2012; Shen & Moore 2014; Shen et al. 2018a; Tanikawa et al. 2018, 2019; Gronow et al. 2020, 2021). Additionally, the thin He envelope that sustains the initial He detonation can lead to an explosion that produces the spectral features of normal Type Ia SNe (Polin et al. 2019; Townsley et al. 2019; Boos et al. 2021; Shen et al. 2021a, 2021b).

Using astrometry from Gaia DR2, Shen et al. (2018b) identified three candidate D6 remnant objects, which we will refer to as D6-1, D6-2 (also known as LP 398-9), and D6-3. These hypervelocity stellar remnants all have Galactocentric velocities in excess of 1000 km s−1. The longer baseline of Gaia EDR3 provides data confirming the status of these three objects as hypervelocity runaway stars with greater precision. The velocities and trajectories of these three objects are difficult to explain with any other mechanism than liberation from a WD+WD binary system that is compact enough (Porb < 10 minutes) for the orbital velocities to exceed 1000 km s−1. It remains an open question what fraction of Type Ia SNe (or perhaps even peculiar thermonuclear SNe) the D6 channel may contribute to (e.g., see recent review by Soker 2019), but the extreme velocities of the observed D6 runaways clearly point to at least some fraction of thermonuclear SNe occurring in very compact WD+WD binaries. Our focus in this Letter is therefore to characterize the constraints on the progenitor WD+WD binaries that likely produced the three runaway objects observed by Shen et al. (2018b).

The D6 scenario specifically predicts an orbital configuration in which the donor WD fills its Roche lobe at the time of SN detonation when it is liberated, so the orbital dynamics and geometry of such a system place a strong constraint on the donor WD radius. Since the Roche-filling donor must be a WD, the WD mass–radius relation then constrains the possible mass of the donor WD based on the measured velocities of the D6 objects.

In this Letter, we present updated velocity distributions for the three candidate D6 donor remnants based on Gaia EDR3 in Section 2. In Sections 3 and 4, we then describe how these velocities can be used along with the WD mass–radius relation to constrain the masses of these objects as donor WDs with the geometry required by the D6 scenario. Surprisingly, none of the resulting mass inferences point to the three candidate remnant objects being descended from standard ≈0.6 M WDs. Instead, the relatively low velocity of D6-2 requires a large radius for its configuration as a Roche-filling donor star that implies it would have been ≲0.2 M as a WD donor unless it was significantly inflated by tides, while the extreme velocities of D6-1 and D6-3 suggest that they were most likely ≈0.8–1.1 M as donor WDs in the D6 scenario. In Section 5 we explore the implications of our calculations for the current measured rotation period of D6-2, and we discuss some further implications of our mass inferences in Section 6.

2. Velocities of the D6 Stars

We compute updated velocities for the D6 stars using data from Gaia EDR3 4 (Gaia Collaboration et al. 2021; Lindegren et al. 2021). We sample the parallax and proper-motion distributions using the ensemble Markov Chain Monte Carlo sampler emcee (Foreman-Mackey et al. 2019), assuming Gaussian uncertainties and taking into account the covariances reported by Gaia. We apply an exponentially decreasing space-density prior on the distance, with a scale length of 1350 pc (e.g., Astraatmadja & Bailer-Jones 2016), and apply an empirical correction 5 to the parallax zero-point (Lindegren et al. 2021). We utilize radial velocities from Shen et al. (2018b) for D6-1 and D6-3, and Chandra et al. (2021) for D6-2. We display the posterior probability distributions of the velocities in Figure 1, along with the previous measurements from Gaia DR2 for comparison. We estimate 1σ uncertainties using the 16th and 84th quantiles of the posterior samples. The final heliocentric velocities are ${2160}_{-160}^{+200}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for D6-1, ${1010}_{-50}^{+60}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for D6-2, and ${2270}_{-400}^{+600}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for D6-3.

Figure 1.

Figure 1. Posterior distributions of total heliocentric velocity for the D6 stars, using astrometry from Gaia EDR3 (solid lines) and Gaia DR2 (dashed lines).

Standard image High-resolution image

In the D6 scenario, the measured velocity of the WD donor remnant is produced almost entirely by its orbital velocity vorb at the time of companion explosion, which is ≳1000 km s−1 in WD+WD binaries with a high-mass primary. The donor WD is also expected to receive a kick from the ejecta with a magnitude of a few hundred km s−1 perpendicular to the direction of orbital velocity, but this can only change the total velocity by a few percent when added in quadrature with the orbital velocity (Wheeler et al. 1975; Taam & Fryxell 1984; Marietta et al. 2000; Hirai et al. 2018; Bauer et al. 2019).

The observed D6 remnants lie within a few kpc of the Sun (Shen et al. 2018b), and we expect that they were most likely ejected from progenitor binaries near the Sun with roughly co-moving Galactic orbits. Therefore, the measured velocities in the heliocentric frame are likely the best representation of the orbital velocities at which they were ejected from the progenitor binary. Accounting for alternative progenitor orbits relative to local Galactic rotation could introduce an additional velocity uncertainty of up to ±200 km s−1 depending on the orientation of observed runaway velocity. For example, when compared to the velocities reported above the heliocentric frame, converting to the Galactocentric frame using astropy (Astropy Collaboration et al. 2018) yields velocities that are ≈200 km s−1 slower for D6-1, ≈100 km s−1 faster for D6-2, and nearly unchanged for D6-3. However, we adopt the heliocentric velocities for the remainder of this work, as they are most likely to best represent the orbital ejection velocities relevant for characterizing the orbital dynamics and geometry of the WD+WD progenitor binary.

3. The Donor Radius

The D6 scenario predicts that the donor star is a Roche-filling WD when its companion explodes as a supernova and ejects it from the system as a hypervelocity runaway. This places a strong constraint on the donor WD radius as a function of its measured velocity. We define the mass ratio of the system as qMdon/Macc, where Mdon is the mass of the donor WD that will become the observed D6 runaway, and Macc is the mass of the accretor WD that explodes as a Type Ia SN. The orbital velocity vorb of the donor star about the center of mass of the binary system can then be expressed as

Equation (1)

where a is the orbital separation. The Roche lobe radius RRL of the donor WD can be expressed using the Eggleton (1983) approximation:

Equation (2)

The fact that the WD donor fills its Roche lobe (Rdon = RRL) at the time of SN detonation allows us to combine the above equations to express the donor radius as

Equation (3)

This expression for the donor radius depends only on the orbital velocity of the donor and the masses of the two stars in the binary. In the D6 scenario for producing a Type Ia SN, the accretor is expected to be a C/O WD with mass in the range Macc ∈ (0.85 M, 1.15 M). We select this accretor mass range to cover the full range of plausible D6 scenarios for producing normal Type Ia SNe (Blondin et al. 2017; Shen et al. 2021a, 2021b).

4. Mass Constraints

The fact that the donor star in the D6 scenario must be a WD implies that the donor radius must also fall along the WD mass–radius relation, which is independent of the radius estimate provided in Equation (3). The measured velocities of the D6 objects therefore provide a constraint on their possible masses as donor stars as shown in Figure 2, where we plot the range of possible radii according to Equation (3) with Macc ∈ (0.85 M, 1.15 M) and vorb from the Gaia EDR3 measurements described in Section 2, along with WD mass–radius relations representing the full range of potential WD temperatures and envelope configurations.

Figure 2.

Figure 2. Donor radii according to Equation (3) compared to WD mass–radius relations. The solid portions of the curves represent the mass–radius relation given by the Bédard et al. (2020) C/O models. The sections of the mass–radius relations covered by models with He cores are indicated by the dashed portions of the curves, and regions where we are extrapolating are indicated by the dotted portions. Models with H envelopes are represented by black curves (H-env), and models with He envelopes are represented by purple curves (He-env). For each D6 object, the inner shaded band represents possible radius solutions using the median posterior velocity with a range of possible accretor masses assuming Mdon < Macc and Macc ∈ (0.85 M, 1.15 M). The outer shaded bands represent the additional radius uncertainty introduced by the ≈1σ velocity uncertainty. D6-3 looks very similar to D6-1 on this plot except that it has a much larger velocity uncertainty, so we only display D6-1 and D6-2 here for clarity.

Standard image High-resolution image

4.1. WD Mass–Radius Relations

The donor WD in the D6 scenario may have a thin outer H envelope when it first fills its Roche lobe and begins to transfer mass. Depending on the mass ratio and compactness of this envelope, dynamically unstable mass transfer may or may not begin right away. If mass transfer immediately proceeds to a dynamically unstable runaway, then the WD will respond to mass loss adiabatically and expand to fill its Roche lobe even as the H envelope is quickly stripped away, so that the radius including this envelope determines the orbital velocity at the time of instability and SN detonation. However, the H envelope may transfer stably for a time. If the envelope can maintain thermal equilibrium, the orbit may continue to shrink toward the more compact radius of the underlying He envelope before dynamical instability sets in. We therefore present mass–radius relations for WD models both with and without H envelopes, which represent an upper and lower bound on the relevant radius of a WD at a given mass and temperature. In reality, some WD+WD systems may fall somewhere between these two limiting cases once mass transfer sets in before ultimately experiencing instability (Shen 2015). Figure 2 shows that the presence of an H envelope has only a small impact on the overall radius of massive WDs, so it will not have a large impact on mass inferences in that regime, but for lower-mass WDs it can substantially influence the result.

WDs in compact binaries may also experience tidal heating that injects energy into the WD interior or envelope and causes it to be inflated relative to the radius of a cooler WD of the same mass (Fuller & Lai 2012; Piro 2019). We use Teff as a proxy for this inflation in Figure 2, since hotter WD models also have more inflated radii, and Teff is easily accessible in the grids of WD models that we use to construct mass–radius relations (Althaus et al. 2013; Bédard et al. 2020). In reality, the interior structure of a tidally heated WD may be more complicated than the structure of a WD cooling model that might produce a similar radius and effective temperature, but Teff provides a loose representation of the possible amount of inflation that an extremely tidally heated WD might achieve.

For models that include H envelopes, we construct mass–radius relations for both cool (Teff = 10,000 K) and hot (100,000 K) donor WDs using the models of Althaus et al. (2013) and Bédard et al. (2020) as follows:

  • 1.  
    For Mdon ≥ 0.45 M, we interpolate from the C/O WD models of Bédard et al. (2020) using Sihao Cheng's WD_models package. 6 We adopt the thick hydrogen envelope models in this regime, so these represent an upper limit for the radius at a given mass, which will translate into a slightly slower velocity at Roche lobe overflow. It is also possible that donor WDs of mass Mdon ≳ 1.05 M would have O/Ne cores that are slightly more compact than the mass–radius relation shown on the high-mass end here. However, the more massive accreting WD in the D6 scenario should have a C/O core to explode as a normal Type Ia SN, and the donor WD needs to be less compact than the accretor, so we adopt the C/O mass–radius relation even for the most massive possible donors. In any case, O/Ne models would only be slightly more compact and therefore would not achieve qualitatively higher orbital velocities before mass transfer.
  • 2.  
    For Mdon < 0.45 M, we interpolate based on the helium-core WD model grid of Althaus et al. (2013), which reliably covers a mass range of about 0.2–0.43 M for Teff≤ 30,000 K. It is also possible that a significant fraction of WDs in compact binaries with Mdon ≈ 0.3–0.5 M could descend from hot subdwarf stars and would therefore have C/O cores rather than He cores (Zenati et al. 2019; Kupfer et al. 2020a, 2020b; Bauer & Kupfer 2021; Romero et al. 2021; Schwab & Bauer 2021). For simplicity, we do not include this possibility in our plots in this paper, but it should be noted that some WDs in this mass range could have somewhat more compact cores than the He models presented here. He cores below 0.3 M cannot sustain He burning, so the lowest mass WDs matching the radii that we infer for D6-2 in Figure 2 should reliably have He cores.
  • 3.  
    For the lower-mass regions where we do not have model coverage, we extrapolate from the lowest mass for which we have reliable models, as shown by the dotted curves in Figure 2. For Teff = 100,000 K, there is no coverage from the lower-mass Althaus et al. (2013) model grid because He-core WDs generally do not reach this high temperature, but the C/O models of Bédard et al. (2020) have radii that extend up to roughly the region of interest for Mdon ≈0.5 M, and a mild extrapolation of the C/O models to lower masses covers the entire relevant radius region. For the lower-temperature models, we extrapolate the He-core model relation when Mdon ≲ 0.2 M.

We follow a similar approach to construct the mass–radius relation for models with no H envelope, using the DB WD models of Bédard et al. (2020) for C/O WDs when Mdon ≥ 0.45 M. For the lower-mass regime where we expect WDs to have He cores, we construct our own model grid using MESA version 15140 (Paxton et al. 2011, 2013, 2015, 2018, 2019). 7 We construct these models starting with a 1.1 M star that ascends the RGB until its He core reaches a mass of 0.35 M. We then remove all of the outer hydrogen envelope to form a semidegenerate He ball of nearly uniform composition, representing a proto-He WD with no H envelope. To make a grid of He WD models, we use a relaxation process to rescale this model to masses in the range 0.1–0.43 M, and then let the models cool to 10,000 K to form a grid of He WDs suitable for a mass–radius relation.

Once again, for the hotter (100, 000 K) set of models, He WDs do not generally reach this temperature during isolated evolution, so we use the Bédard et al. (2020) DB WD models for Mdon ≥ 0.45 M, and we show a rough extrapolation of that mass–radius relation for lower WD masses at this temperature.

4.2. Donor Mass from Measured Velocity

The donor mass can be expressed as a function of donor velocity by setting Rdon from Equation (3) equal to the WD mass–radius relation. We show results based on WD mass–radius relations at both cool and hot effective temperatures in Figure 3 for models that include H envelopes, and in Figure 4 for models that have no H envelope. It is immediately clear that D6-1 and D6-3 require rather massive donor WDs to produce their high velocities, with their median velocities favoring Mdon ≈ 1.0 M depending on Teff and envelope composition.

Figure 3.

Figure 3. Solutions for WD donor mass as a function of velocity assuming a cool WD (Teff = 10, 000 K) mass–radius relation in the left panel and a hot WD (Teff = 100, 000 K) relation in the right panel. These mass–radius relations include H envelopes in both cases (see Figure 4 for the He envelope case). The gray band gives the possible range of solutions for accretor masses in the range Macc ∈ (0.85 M, 1.15 M), subject to the constraint that Mdon < Macc. The blue, orange, and green shaded regions give the ±1σ velocity ranges for D6-1, D6-2, and D6-3, respectively, and the corresponding vertical lines mark the median posterior velocity for each object. The lower panels show the Gaia EDR3 velocity posterior distributions from Figure 1 (arbitrary scale).

Standard image High-resolution image
Figure 4.

Figure 4. Same as Figure 3 but using the mass–radius relations for WDs with no H envelope, i.e., with He envelopes.

Standard image High-resolution image

WDs in binaries should generally have ample time to cool to Teff ≈ 10,000 K before coming into contact due to gravitational wave radiation and beginning mass transfer (e.g., Maoz et al. 2018; Cheng et al. 2020). The left panels of Figures 3 and 4 would then suggest that D6-2 was an extremely low-mass WD with a He core and Mdon ≲ 0.2 M to explain its low observed velocity. On the other hand, the right panels show that D6-2 is also marginally consistent with having been a very hot ≈0.5 M C/O WD that was inflated when mass transfer began due to its high temperature. This configuration could possibly be produced by tidal heating in the compact WD binary as it approached contact (Fuller & Lai 2012; Burdge et al. 2019; Piro 2019), though theoretical predictions for the total amount of tidal heating and corresponding inflation in the outer layers are complex and still uncertain in this regime (Burkart et al. 2013, 2014; Fuller & Lai 2013, 2014; Yu et al. 2020). According to Piro (2019), the total luminosity available from tidal heating can be on the order of ≈1 L, which could translate into a temperature of Teff ≈ 30– 40,000 K for a WD with the radius we calculate for D6-2. This could contribute somewhat toward the inflation of a WD in a binary system like D6-2. However, in the Porb = 6.9 minute WD+WD system ZTF J1539+5027 observed by Burdge et al. (2019), the lower-mass WD that will eventually become the donor is very cool and inconsistent with significant tidal heating, while the higher mass WD in that system is much hotter. Burdge et al. (2019) also argue that it is difficult for tidal heating to explain a temperature much hotter than 40,000 K for a WD+WD system at this orbital period.

Tidal heating has a much smaller impact on mass inferences for D6-1 and D6-3, where it only pushes to slightly higher masses above ≈1.0 M according to the right panels in Figures 3 and 4.

5. Rotation

Chandra et al. (2021) have recently identified a period for the D6-2 remnant that most likely corresponds to its current rotation period. The current radius of the D6-2 donor remnant is significantly inflated compared to its radius prior to the SN companion explosion, and Chandra et al. (2021) discuss the relationship between the current observed rotation rate and the potential tidally locked rotation rate of the D6-2 donor in the binary just prior to SN detonation. The required radius for D6-2 in Figure 2 further refines the possible inferences that can be made about rotation in the D6-2 system. Assuming tidal locking for D6-2 as a donor star, conserving angular momentum in its outer layers predicts that its current rotation period should be

Equation (4)

where Robs = 0.2 R is the current measured radius of the D6-2 remnant object according to Chandra et al. (2021). The system likely would have had an orbital period of Porb = 3–7 minutes at the moment of SN detonation, and therefore Equation (4) and Figure 2 predict that the current observed rotation period should be Prot ≈ 2–5 hr. However, Chandra et al. (2021) report a measured rotation period of Prot ≈ 15.4 hr. This may require that D6-2 experienced some angular momentum loss after interacting with SN ejecta to reach its current rotation period. Alternatively, because D6-2 was so much less compact than the other D6 objects, its more loosely bound surface layers could have experienced significant stripping when interacting with SN ejecta. The current surface could have originally been deeper material at a smaller radius that inflated more to reach the current observed radius, so that Equation (4) would predict a longer rotation period than the value given by using Rdon.

6. Discussion

We have calculated WD donor masses for D6 runaways using observed Gaia EDR3 velocities and WD mass–radius relations for several combinations of possible envelope configurations and effective temperatures. Our conclusions on the high donor mass (Mdon ≳ 0.8 M) for both D6-1 and D6-3 hold across all effective temperatures and nondegenerate envelope layers we sample. Since the donor must be less massive than the accretor, this implies that two of the three known D6 stars had progenitor binaries composed of two massive white dwarfs, both likely >0.8 M. No massive WDs are currently known to show large radial-velocity changes indicating another massive companion, though such searches are challenging due to their intrinsic faintness (Rebassa-Mansergas et al. 2019). Additionally, most radial-velocity surveys have limited numbers of massive white dwarfs; just 31 out of 643 targets (<5%) in the ESO SN Ia Progenitor surveY (SPY) have a mass >0.8 M (Napiwotzki et al. 2020). If WD+WD systems containing two massive WDs are indeed intrinsically rare, this may suggest that the D6 channel only contributes a fraction of explosions to the overall thermonuclear supernova rate.

The low donor mass of the D6-2 system in Figure 3 implies that it would have originated in a system with a low mass ratio. A system with such a low mass ratio would not necessarily experience dynamical instability at the onset of mass transfer (Marsh et al. 2004; Dan et al. 2012). However, this dynamical instability is a prerequisite for producing the detonation that would eventually lead to the liberation of the donor as a hypervelocity runaway. Therefore, the existence of the runaway D6-2 donor remnant lends further support to the arguments of Shen (2015) and Brown et al. (2016) that most low-mass-ratio double-WD binaries are disrupted at contact rather than experiencing stable mass transfer from an ELM WD to form an AM CVn system (e.g., Wong & Bildsten 2021).

Alternatively, the qualitatively slower velocity of D6-2 compared to the other D6 candidates may instead suggest that it had a different compact binary evolution origin. For example, the highest velocity tail of the runaway velocity distribution for He star donors to WD companions calculated by Neunteufel (2020) and Neunteufel et al. (2021) could be marginally compatible with the velocity of D6-2. This could require an unlikely random alignment of orbital ejection with a halo star orbit, as has been suggested for the runaway He star US 708 (Brown et al. 2015; Geier et al. 2015; Bauer et al. 2019; Liu et al. 2021). However, Hermes et al. (2021) have argued that relatively slow rotation (Prot ≳ 10 hr) is problematic for the runaway He star donor channel, and it is also somewhat difficult to explain for a runaway WD donor as discussed in Section 5. Instead, D6-2 could plausibly be related to the LP 40–365 class of runaway stars, which have been suggested to be bound remnants of Type Iax SNe (Jordan et al. 2012; Vennes et al. 2017; Raddi et al. 2018a, 2018b, 2019; Hermes et al. 2021). The three candidates identified as part of that class so far have velocities somewhat below 1000 km s−1, but the distributions of kick velocities in that scenario are uncertain, and predictions for final velocities of remnants are not constrained well enough to completely exclude a velocity as high as that of D6-2. Pakmor et al. (2021) have also proposed yet another possibility for producing runaways with velocities up to 1000–1500 km s−1 from binaries containing a massive He-C/O WD (e.g., Zenati et al. 2019; Pelisoli et al. 2021) donating material onto a C/O WD companion. In that case, the donor's thick He envelope may experience a detonation and release the accretor WD as a hypervelocity runaway. This speculative scenario may achieve velocities sufficient to explain D6-2, but not D6-1 or D6-3.

In conclusion, our calculations in this work have enabled new inferences about the masses of the observed D6 runaways that fall on either end of the WD mass spectrum, while objects near the 0.6 M peak of the WD mass distribution are absent so far. The relatively low velocity of D6-2 points to a low donor WD mass, introducing new puzzles for its interpretation as a D6 donor remnant in the context of growing observational data about its present-day properties, such as surface rotation and circumstellar material (Chandra et al. 2021). On the other end of the velocity scale, D6-1 and D6-3 both consistently point toward origins in massive WD+WD binaries regardless of modeling assumptions. This clearly motivates further observational followup work, especially detailed spectroscopic analysis to understand the complicated surface compositions and structures of D6 stars.

We thank Warren Brown and Charlie Conroy for comments and suggestions. V.C. is supported in part by the James Mills Peirce fellowship at Harvard University. Financial support for K.J.S. and J.J.H. was in part provided by NASA/ESA Hubble Space Telescope program #15871. K.J.S. also received support from NASA through the Astrophysics Theory Program (NNX17AG28G). The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.

Footnotes

Please wait… references are loading.
10.3847/2041-8213/ac432d