From Neutrino- to Photon-cooled in Three Years: Can Fallback Accretion Explain the X-Ray Excess in GW170817?

The neutron star merger GW170817 (Abbott et al. 2017a) was accompanied by radiation covering the electromagnetic spectrum, from radio to gamma-rays (Abbott et al. 2017b; Margutti & Chornock 2020). Thermal kilonova emission, powered by the radioactive decay of r-process elements (Li & Paczyński 1998; Metzger et al. 2010b; Barnes & Kasen 2013), was observed starting 11 hr after the merger (e.g., Arcavi et al. 2017; Coulter et al. 2017; Soares-Santos et al. 2017) and continued to be detected in the infrared up to 74 days (Villar et al. 2018; Kasliwal et al. 2019). The merger was also accompanied by a broadband nonthermal synchrotron afterglow (Haggard et al. 2017; Hallinan et al. 2017; Margutti et al. 2017), created by the shock interaction of a relativistic jet interacting with the circum-merger medium (e.g., Gottlieb et al. 2018; Granot et al. 2018; Lamb & Kobayashi 2018; Lazzati et al. 2018; Wu & MacFadyen 2019). Such relativistic jets are believed to give rise to short-duration gamma-ray bursts at cosmological distances, powered by the accretion of neutron star debris onto the newly formed black hole (BH; e.g., Narayan et al. 1992).

After peaking at t ≈ 150 days post-merger, the nonthermal radio and X-ray emission began to fade as a steep power law in time, which is consistent with theoretical predictions of the afterglow of a jet initially directed away from our line of site (Margutti et al. 2018; Nynka et al. 2018; Piro et al. 2019). However, a significant excess relative to the afterglow model was recently reported based on Chandra X-ray Observatory observations at t ≈ 3.4 yr (Balasubramanian et al. 2021; Hajela et al. 2021; Troja et al. 2021). A similar excess is not seen in the radio band, requiring a spectral flattening relative to the nearly fixed power law spectral shape observed at earlier epochs. Such a flattening is not expected in the standard afterglow model in which the spectral index of the electrons is fixed (e.g., Hajela et al. 2021). Instead, its presence suggests the onset of a new component of emission that predominantly contributes with an X-ray luminosity of L_X ∼ 5 × 10³⁸ erg s⁻¹ in the ∼keV range, albeit with significant uncertainty due to low photon statistics and sensitivity to the method of instrumental calibration.

Two main ideas have been put forward for generating an extra late component of emission from a neutron star merger. One possibility is that this component represents a different source of synchrotron emission, powered by the shock interaction of the kilonova ejecta with the interstellar medium (e.g., Nakar & Piran 2011; Kathirgamaraju et al. 2019; Nedora et al. 2021). Although the bulk of the merger ejecta expanded at velocities ≈0.1–0.2c (e.g., Villar et al. 2017), the dynamical ejecta can possess a small quantity of higher velocity material. The need to generate enough high-velocity matter to explain the X-ray luminosity can in principle constrain the neutron star equation of state (Nedora et al. 2021), to which the dynamical ejecta is sensitive (e.g., Radice et al. 2018). However, a drawback of this scenario is that it requires the shock-accelerated electrons to possess a flat energy distribution (p ≲ 2.1, where ${dN}/{dE}\propto {E}^{-p}$), in tension with the larger value p ≈ 2.2 inferred from the earlier afterglow (e.g., Margutti et al. 2018) and from other non-relativistic shocks such as those in supernovae (e.g., van Dyk et al. 1994).

Another potential source of X-ray emission occurs on much smaller scales: the accretion disk of the newly formed BH remnant ⁴ being fed by the fallback of debris from the merger (Hajela et al. 2021; Ishizaki et al. 2021a). As discussed by Hajela et al. (2021), an accretion-powered X-ray source has several appealing features. First, the X-ray luminosity is below the Eddington luminosity, L_Edd ≈ 8 × 10³⁸ erg s⁻¹ of the ≈ 2.6M_⊙ BH remnant (Abbott et al. 2017a), which is consistent with the generation of a radiatively efficient accretion flow down to the innermost stable circular orbit (ISCO) of the spinning stellar-mass BH. Second, the predicted temperature of the thermal disk emission at this luminosity is kT_eff ∼ 1 keV, which is close to that of the observed excess. Furthermore, the optical depth of X-ray photons through the kilonova ejecta can reach ≲1 on a timescale of years after the merger (e.g., Margalit et al. 2018), unshrouding the BH accretion funnel.

Perhaps the biggest open question regarding the accretion scenario is whether sufficient mass is returning to the BH years after the merger. Even assuming an initial binary with an unequal mass ratio, the quantity of dynamical ejecta falls short by orders of magnitude (e.g., Rosswog 2007). On the other hand, the large mass of the kilonova ejecta from GW170817 (M_ej ≈ 0.03–0.08M_⊙; Arcavi et al. 2017; Cowperthwaite et al. 2017; Villar et al. 2017) excludes a dynamical origin, but is consistent with arising in outflows from the post-merger accretion disk (e.g., Metzger et al. 2008, 2009; Lee et al. 2009; Fernández & Metzger 2013; Just et al. 2015; Siegel & Metzger 2017; Fujibayashi et al. 2018; Fernández et al. 2019). If a mass comparable to M_ej were to remain gravitationally bound to the BH and return at a sufficiently gradual rate, the observed X-ray excess could be produced (Hajela et al. 2021; Ishizaki et al. 2021a). On the other hand, if the BH accretion rate is too large on intermediate timescales (days to months, when the ejecta is still opaque at X-ray wavelengths), this could over-produce the kilonova luminosity via X-ray reprocessing (e.g., Kisaka et al. 2016; Matsumoto et al. 2018).

In this Letter we address these issues by means of hydrodynamical simulations of the post-merger BH accretion flow. Adopting disk models consistent with parameters inferred from GW170817 and numerical relativity simulations, and which are consistent with the evolution found in full magnetohydrodynamics (MHD; e.g., Fernández et al. 2019), we run the simulations for an unprecedentedly long time ≈35 s (≳10⁶ dynamical times at the BH horizon), as necessary to follow the evolution of the marginally bound debris into the full radiatively inefficient stage. Extrapolating our findings to timescales ∼ yr, we demonstrate that fallback from the disk outflow may in principle be sufficient to account for the X-ray excess of GW170817, while simultaneously evading an overproduction of the kilonova emission.

We simulate the long-term evolution of the accretion disk from GW170817 assuming that a BH formed promptly. The time-dependent evolution is carried out in axisymmetry with the hydrodynamic code FLASH version 3.2 (Fryxell et al. 2000; Dubey et al. 2009). The public version has been modified to include angular momentum transport with an imposed shear stress, following the prescription of Shakura & Sunyaev (1973) and Stone et al. (1999), neutrino emission and absorption using a leakage scheme for emission and an annular lightbulb for absorption (Fernández & Metzger 2013; Metzger & Fernández 2014), and a pseudo-Newtonian potential to model the gravity of the black hole (Artemova et al. 1996; Fernández et al. 2015). The equation of state is that of Timmes & Swesty (2000), modified such that neutrons, protons, and alpha particles are in nuclear statistical equilibrium for T > 5 × 10⁹ K, and accounting for changes in nuclear binding energy (nuclear dissociation and recombination).

The initial condition is a BH of mass M_bh = 2.65M_⊙ and dimensionless spin 0.8, as well as an equilibrium torus with initial mass 0.1M_⊙, which is consistent with values inferred for GW170817 (Abbott et al. 2017c; Shibata et al. 2017). The torus has constant initial electron fraction 0.1, entropy 8 k_B per baryon, and constant specific angular momentum. The domain outside the torus is filled with a low-density ambient medium. The boundary conditions are outflow in radius and reflecting in polar angle at the symmetry axis. The radial range extends from a radius halfway between the ISCO and the horizon, out to a radius 10⁵ times larger. The grid is discretized with 640 cells logarithmically spaced in radius and with 112 cells uniformly spaced in $\cos \theta$.

We run two simulations using values of the viscosity parameter α = {0.03, 0.1}, which bracket the late-time (t > 1 s) behavior of the accretion rate in general relativistic magnetohydrodynamic (GRMHD) simulations (Fernández et al. 2019). Both models are evolved for 10,000 orbits at the initial density maximum of the torus (50 km), corresponding to ∼35 s. By this time, neutrino emission has decreased to negligible levels and and nuclear dissociation/recombination has ceased, because the maximum temperature in the disk is less than 5 × 10⁹ K. As summarized in Table 1, the masses ejected with positive energy ⁵ at a distance 10⁹ cm from the BH by the α = {0.03, 0.1} simulations are {1.7 × 10⁻², 3 × 10⁻²}M_⊙, respectively, on the lower end of the range of r-process production inferred from the kilonova of GW170817 (e.g., Villar et al. 2017).

Table 1. Simulation Parameters and Results

Model	Mbh	Md	α	tmax	Mej	${M}_{\mathrm{bnd}}({t}_{\max })$	${\left.\tfrac{{dM}}{{dE}}\right|}_{E=0,{t}_{\max }}$	${\dot{M}}_{\mathrm{bh},10}$	βbh	${\dot{M}}_{\mathrm{fb},10}$	βfb	LX
	(M⊙)	(M⊙)		(s)	(M⊙)	(M⊙)	(g2 erg−1)	(M⊙ s−1)		(M⊙ s−1)		(1038 erg s−1)
v03	2.65	0.10	0.03	35	1.7E-2	2.1E-3	1.1E+14	2.1E-5	2.41	4.0E-4	1.43	2–70
v10	2.65	0.10	0.10	35	3.0E-2	4.9E-4	1.7E+13	1.0E-5	2.10	1.0E-4	1.72	0.1–0.5

Note. Columns from left to right show: model name, BH mass, initial disk mass, viscosity parameter, maximum simulation time, ejected mass with positive net energy, mass in the computational domain with negative net energy at t_max, dM/dE around E = 0 and t_max (Figure 3), normalization and slope of power-law fit to accretion rate at the ISCO radius (Equation (1)), normalization and slope of power-law fit to the rate of change of bound mass in the domain (Equation (2)), and X-ray luminosity at t = 3.4 yr (t ≃ 10⁸ s), calculated from Equation (4) adopting the normalization ${\dot{M}}_{\mathrm{fb},10}$ at t = 10 s, and for power-law decay indices ranging from β = β_fb to β = 5/3.

Download table as:  ASCIITypeset image

Figures 1–3 summarize our results. We focus on the α = 0.03 model, insofar as it more closely matches the late-time accretion rate into the BH in GRMHD simulations (Fernández et al. 2019) and is the most promising in terms of generating a larger quantity of late-time fallback accretion (we return to a discussion of the α = 0.1 model at the end of the section). In Figure 1, a black solid line shows the accretion rate at the ISCO ${\dot{M}}_{\mathrm{bh}}$, while the red solid line shows the rate of decrease of gravitationally bound matter in the computational domain, ${\dot{M}}_{\mathrm{fb}}$. The difference, ${\dot{M}}_{\mathrm{out}}={\dot{M}}_{\mathrm{fb}}-{\dot{M}}_{\mathrm{bh}}$, represents mass unbound in outflows from the disk. Dotted lines in Figure 1 show power-law extrapolations of ${\dot{M}}_{\mathrm{bh}}$ and ${\dot{M}}_{\mathrm{fb}}$ to later times than followed by our simulations, of the form (see Table 1)

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{bh}}={\dot{M}}_{\mathrm{bh},10}{\left(\displaystyle \frac{t}{10\,{\rm{s}}}\right)}^{-{\beta }_{\mathrm{bh}}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{fb}}={\dot{M}}_{\mathrm{fb},10}{\left(\displaystyle \frac{t}{10\,{\rm{s}}}\right)}^{-{\beta }_{\mathrm{fb}}}\end{eqnarray} \tag{ 2 }$$

where the fit parameters are obtained from t = 10 s until the end of each simulation. Note that the outflow rate inferred from the decrease of the total bound mass in the computational domain computed here (${\dot{M}}_{\mathrm{fb}}$) differs from the unbound mass flux at a fixed sampling radius, as typically reported in disk studies, as the latter decays more steeply with time than the former as the bound disk expands viscously past the sampling radius (e.g., Figure 2) and outflows emitted at larger radii are missed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of the torus follows that described in previous works (e.g., Fernández & Metzger 2013). At early times, the accretion rate is higher than the characteristic value ${\dot{M}}_{\mathrm{ign}}\sim {10}^{-2}{M}_{\odot }$ s⁻¹ (brown dotted–dashed line in Figure 1) above which the disk is sufficiently hot and dense to cool by neutrinos (e.g., Narayan et al. 2001; Chen & Beloborodov 2007; De & Siegel 2020), resulting in a high accretion efficiency (${\dot{M}}_{\mathrm{fb}}\approx {\dot{M}}_{\mathrm{bh}};$ ${\dot{M}}_{\mathrm{out}}\ll {\dot{M}}_{\mathrm{bh}}$). However, at later times as the accretion rate drops ${\dot{M}}_{\mathrm{bh}}\lt {\dot{M}}_{\mathrm{ign}}$, weak interactions become slow and freeze out, causing a transition to a radiatively inefficient state (e.g., Metzger et al. 2008). Once the maximum temperature in the disk decreases below 5 × 10⁹ K, internal energy changes due to nuclear dissociation or recombination of alpha particles also stops. At this point, viscous heating is no longer balanced by any significant cooling process and the disk becomes susceptible to outflows that carry much of its mass in unbound outflows, i.e., ${\dot{M}}_{\mathrm{out}}\approx {\dot{M}}_{\mathrm{fb}}$, suppressing accretion onto the BH ${\dot{M}}_{\mathrm{bh}}\ll {\dot{M}}_{\mathrm{fb}}$ (e.g., Blandford & Begelman 1999; Quataert & Gruzinov 2000).

Figure 2 shows a snapshot of the net specific energy E at the end of the α = 0.03 simulation (t ≈ 35 s). In the equatorial plane, gravitationally bound matter (E < 0) released by outflows during the early phases of the torus evolution, extends out to radii r_out ∼ 10¹⁰ cm where E ≈ 0. The marginally bound matter at these large radii is not rotationally supported, but instead moves ballistically with a freefall time ${t}_{\mathrm{orb}}\sim 2\pi {({r}_{\mathrm{out}}^{3}/{GM})}^{1/2}$ comparable to the system age ∼ t. Fresh matter can become unbound from the system (E > 0) only as a result of viscous heating, which occurs primarily in the rotationally supported disk on smaller scales r_disk ≲ 3 × 10⁷ cm. Pockets of matter that reach positive energy at r ∼ r_disk expand outward along the top and bottom boundary between the disk and polar funnel, where it mixes and joins a quasi-spherical unbound outflow at r ≳ r_out (see, e.g., Figure 5 in Fernández et al. 2015).

With this picture in mind, our interpretation for the evolution of the unbound mass, ${\dot{M}}_{\mathrm{fb}}$, is as follows. Marginally bound debris, released during early phases in the disk evolution and extending out to radii r_out set by the condition t_orb ∼ t, is continuously falling back to join the inner rotationally supported disk. However, upon returning this material is promptly unbound by viscous heating due to the radiatively inefficient nature of the accretion flow and the low binding energy of the fallback matter at late times (e.g., Blandford & Begelman 1999; Li et al. 2013), as suggested by Rossi & Begelman (2009). Thus, matter becomes unbound at a rate that follows the rate it falls back to the inner disk.

To check this interpretation, Figure 3 shows a mass histogram of net specific energy E at three different times in the α = 0.03 simulation. Of particular interest is the value of dM/dE near E ≈ 0, as this quantity controls the rate of mass fallback to the disk of the marginally bound debris. If mass falls back to the disk ballistically on the orbital time $t={t}_{\mathrm{orb}}=2\pi {({a}^{3}/{{GM}}_{\mathrm{bh}})}^{1/2}$, where a is the semimajor axis and ∣E∣ = GM_bh/2a, then the fallback rate is given by (e.g., Rees 1988)

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}\approx 4\times {10}^{-4}{M}_{\odot }\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{dM}/{dE}{| }_{E\approx 0}}{{10}^{14}\,{{\rm{g}}}^{2}\,{\mathrm{erg}}^{-1}}\right){\left(\displaystyle \frac{t}{10\,{\rm{s}}}\right)}^{-5/3}.\end{eqnarray} \tag{ 3 }$$

By the final two snapshots, the energy distribution is not evolving significantly around E ≈ 0. Using the value dM/dE∣_E≈0 ≈ 1.2 × 10¹⁴ g² erg⁻¹ inferred from the final snapshot, we find dM/dt ≈ 5 × 10⁻⁴ M_⊙ at t = 10 s. This is in excellent agreement with the rate at which matter is being unbound from the disk on the same timescale (${\dot{M}}_{\mathrm{fb},10};$ Table 1). However, the best-fit power-law index for ${\dot{M}}_{\mathrm{fb}}$ of β_fb ≈ 1.43 (Equation (2)) is somewhat shallower than the value ≈5/3 predicted by Equation (3) for constant dM/dE.

Zoom In Zoom Out Reset image size

Naive extrapolation of the BH accretion rate (dotted black line in Figure 1) to the present epoch t ≈ 3.4 yr would fall short of that required to explain the observed X-ray excess by many orders of magnitude, if one assumes

$$\begin{eqnarray}&&{L}_{{\rm{X}}}\approx \eta {\dot{M}}_{\mathrm{bh}}{c}^{2},\end{eqnarray} \tag{ 4 }$$

where radiative efficiency of η ≈ 0.1 is the maximum value for a BH of spin 0.8 (Bardeen et al. 1972). However, by the epoch t ≈ 3.4 yr, the rate of fallback accretion has become sub-Eddington (brown dotted–dashed line in Figure 1), at which point the photons are no longer trapped in the flow and are able to diffuse out, radiatively cooling the disk (e.g., Begelman 1979). Thus, as with the early neutrino-cooled epoch, we expect that ${\dot{M}}_{\mathrm{bh}}$ will again rise to match ${\dot{M}}_{\mathrm{fb}}$ starting a few years after the merger, in which case ${L}_{{\rm{X}}}\approx \eta {\dot{M}}_{\mathrm{fb}}{c}^{2}$ for L_X ≲ L_Edd ≈ 8 × 10³⁸ erg s⁻¹. This transition is illustrated schematically with a dashed line in Figure 1.

Interpreted in this way, we see that the BH accretion rate at t ≈ 3.4 yr may indeed be sufficient to explain the X-ray excess from GW170817, which we show for comparison with a blue square in Figure 1 (Hajela et al. 2021; Troja et al. 2021). We include a larger error bar on the luminosity value than the formal one (by a factor of 2 in either direction), to account for uncertainties in (1) the bolometric correction of the disk emission into the Chandra bandpass; (2) residual absorption by the kilonova ejecta; (3) geometric beaming of the X-ray emission by the ∼ L_Edd accretion flow.

The lateness of the transition to a radiative efficient disk is also important for not overproducing the kilonova emission. The kilonova ejecta will remain opaque to X-rays from the inner accretion disk for at least the first year of the explosion, only becoming optically thin around the present epoch (Margalit et al. 2018; Hajela et al. 2021). Most of the accretion power that emerges during earlier epochs will be absorbed and reprocessed into optical/infrared emission (e.g., Kisaka et al. 2016; Matsumoto et al. 2018). It is thus of equal importance that the maximal accretion power $\sim 0.1{\dot{M}}_{\mathrm{bh}}{c}^{2}$ remains below the observed luminosity of the kilonova emission. The latter is shown in Figure 1 by blue squares where ultraviolet/optical/infrared (UVOIR) data is available (e.g., Cowperthwaite et al. 2017), and as blue triangles where 4.5 μm detections serve as lower limits (Villar et al. 2018; Kasliwal et al. 2019). Although the results are sensitive to how rapidly the flow becomes radiatively efficient, we cannot rule out that some or all of the infrared emission detected by Spitzer is powered by accretion instead of radioactivity.

Although our α = 0.03 model appears promising to explain the X-ray excess of GW170817, the α = 0.1 model is less so (see Table 1). The rate at which matter is unbound (∼falling back) in our α = 0.1 simulations decays more steeply, β_fb ≈ 1.7 compared to the α = 0.03 case, while its normalization (derived either directly from the rate at which fallback is unbound or from the energy distribution dM/dE around E ≈ 0), is a factor of ≈3–7 times lower. We surmise that the reason for this difference is that a higher value of α leads to more vigorous heating of the debris that returns to the rotationally supported disk, thereby launching more powerful outflows that act to unbind material that would otherwise fall back. This is consistent with the larger mass ejected by the model with α = 0.1 (Table 1). We also find that at late times the α = 0.1 model displays large amplitude (factor ∼10) oscillations in the accretion rate on timescales of ∼10 s, making our power-law fits to the temporal dependence less reliable than in the α = 0.03 case.

An obvious caveat to our conclusions is that we have employed an α-viscosity prescription in lieu of a self-consistent treatment of angular momentum transport via the magneto-rotational instability (MRI). MHD simulations of the post-merger disk evolution that resolve the MRI find larger mass ejection than the α-viscosity hydro simulations (e.g., Siegel & Metzger 2017; Christie et al. 2019; Fernández et al. 2019). Mass outflow rates are usually measured at a fixed extraction radius, but as pointed out earlier, over long timescales the bound portion of the disk expands outward beyond any extraction radius and outflow rates cut off steeply with time at that location. While the time dependence of outflow rates in MHD and hydrodynamic simulations is similar when measured in this way (Fernández et al. 2019), one cannot assume that the time dependence of the bound mass in the computational domain will be the same owing to the additional mass ejection provided by magnetic processes other than dissipation of MRI turbulence, which are not present in hydrodynamic simulations. On the other hand, accretion disk masses up to ≈ 0.2M_⊙, which are twice as large as those adopted in our simulations, are compatible with numerical relativity merger simulations and the quantity of disk wind ejecta in GW170817 (e.g., Shibata et al. 2017). Our simulations also neglect the impact of additional heating of the bound debris by the late-time radioactive decay of r-process nuclei, which can potentially impact the matter energy distribution and rate of late-time fallback (e.g., Metzger et al. 2010a; Desai et al. 2019; Ishizaki et al. 2021b).

The surprisingly close neutron star merger GW170817 offers a unique opportunity to witness a single accretion flow evolve from an neutrino-cooled phase in the first second after the merger to a photon-cooled state a few years later, in between which is a long period of radiatively inefficient accretion. Using long-term axisymmetric hydrodynamical simulations of the post-merger accretion disk system, we have taken preliminary steps to evaluate the feasibility of late-time accretion onto the newly formed BH as a source of the recently discovered excess X-ray luminosity observed at 3.4 yr after the merger. Taking the rate of decline of bound gaseous mass as a proxy for the rate at which matter is returning to the disk at late times (and being unbound in disk winds), we find that a power-law extrapolation of ${\dot{M}}_{\mathrm{fb}}$ in our α = 0.03 simulation, yields a late-time accretion rate onto the BH, which is in principle sufficient to explain the excess.

One prediction of the accretion-powered origin is that the X-ray excess should, possibly after a brief brightening period (if the ejecta is still marginally opaque to X-rays at the present epoch), begin to decay, roughly as L_X ∝ t^−β where β ≈ 5/3 (Hajela et al. 2021). This contrasts with alternative kilonova afterglow model, for which L_X is predicted to decay slower, or even rise, in time (e.g., Kathirgamaraju et al. 2019; Nedora et al. 2021), and should eventually be accompanied by radio emission in excess of the ordinary afterglow contribution.

B.D.M. acknowledges support from NSF (grant AST-2002577) and NASA (grants NNX17AK43G, 80NSSC20K0909). R.F. acknowledges support from the Natural Sciences and Engineering Research Council (NSERC) of Canada through Discovery Grant RGPIN-2017-04286, and by the Faculty of Science at the University of Alberta. The software used in this work was in part developed by the U.S. Department of Energy (DOE) NNSA-ASC OASCR Flash Center at the University of Chicago. Data visualization was done in part using VisIt (Childs 2012), which is supported by DOE with funding from the Advanced Simulation and Computing Program and the Scientific Discovery through Advanced Computing Program. This research was enabled in part by support provided by WestGrid (www.westgrid.ca), the Shared Hierarchical Academic Research Computing Network (SHARCNET, www.sharcnet.ca), Calcul Québec (www.calculquebec.ca), and Compute Canada (www.computecanada.ca). This research also used resources of the U.S. National Energy Research Scientific Computing Center (NERSC), which is supported by the DOE Office of Science under Contract No. DE-AC02-05CH11231 (repository m2058).