THE THREE-DIMENSIONAL EVOLUTION TO CORE COLLAPSE OF A MASSIVE STAR

Real stars are not truly spherically symmetric. This is especially true for the interiors of massive stars at the end of their lives. As massive stars approach core collapse, the equation of state (EOS) becomes softer, cooling by neutrino emission drives nuclear burning ever more vigorously, and convective velocities increase. Nuclear burning couples to turbulent convection so that fuel is consumed in chaotic bursts. Core burning and thick shell burning are dominated by large scale modes of flow, which are of such low order that they do not cancel to a smooth spherical behavior. In 1D stellar evolution codes, convective mixing and energy transport is modeled using mixing-length theory (MLT), which is tuned to reproduce the solar photosphere (Asplund et al. 2009). Three-dimensional (3D) simulations of turbulent convection, however, demonstrate that MLT gives flawed representations of stellar convection, especially during the late burning stages (Arnett et al. 2015), and two-dimensional (2D) simulations of the late stages of massive stellar evolution have shown that the Si burning in the shell surrounding the iron core is strong and violent, generating large-scale fluctuations in the core of the star that will be present at the point of core collapse (Bazan & Arnett 1994; Arnett & Meakin 2011).

The state-of-the art in core-collapse supernova (CCSN) progenitor models (e.g., Woosley et al. 2002; Woosley & Heger 2007) is still 1D, and so simulations of the CCSN mechanism have all but exclusively used 1D initial conditions (ICs). Recent exploration of the impact of multidimensional progenitor structure on the CCSN mechanism (Couch & Ott 2013, 2015; Fernández et al. 2014; Müller & Janka 2014) has shown that physically motivated, yet highly parameterized, non-spherical structure in otherwise 1D progenitor models can have a qualitative and generally favorable impact on the CCSN mechanism. Clearly, there is a need for realistic 3D progenitor models that self-consistently capture the turbulent nuclear burning.

In this Letter, we present the first 3D simulation of the late stages of iron core growth via strong Si shell burning through the moment of core collapse. We concentrate on the final three minutes of Si burning in a 15 M_☉ star evolved initially in 1D. This timescale captures several (∼8) eddy turnover times of the convection in the Si shell and the synthesis of 0.2 M_☉ of iron. Gravitational core instability is attained when the iron core reaches its effective Chandrasekhar mass, which depends in general on the core's electron fraction and entropy (Baron & Cooperstein 1990). We find the convection driven by Si shell burning to be strong with vigorous fluctuations. The turbulent speeds reach several 100 $\mathrm{km}\;{{\rm{s}}}^{-1}$, with the largest eddies being roughly the full width of the Si shell. In order to assess the impact of realistic, physically self-consistent 3D progenitor structure on the CCSN mechanism, we follow the collapse of the 3D progenitor through core bounce and shock revival using approximate neutrino transport methods. Compared to spherically symmetric ICs, the presence of realistic 3D structure results in more favorable conditions for CCSN shock revival and robust explosion. We find that this is due principally to stronger post-shock turbulence, which provides a greater effective turbulent pressure that aids shock expansion (Murphy et al. 2013; Couch & Ott 2015). Turbulence in the context of the CCSN mechanism is generated by neutrino-driven buoyant convection (Murphy & Meakin 2011; Abdikamalov et al. 2014) and aspherical shock motion (Endeve et al. 2012). Progenitor asphericity introduces finite amplitude perturbations that enhance the growth rate of the instabilities that drive turbulence (e.g., Foglizzo et al. 2006). This first result shows that the final minutes of massive stellar evolution can, and should, be simulated in 3D.

The study of the CCSN mechanism with high-fidelity 3D simulations is still in its infancy. Early results, however, indicate that progenitors that explode successfully in 2D (e.g., Müller et al. 2012; Bruenn et al. 2014) may not in 3D (e.g., Tamborra et al. 2014). There are multiple physical and numerical reasons why we should expect 2D simulations to be artificially prone to explosion (Hanke et al. 2012; Couch 2013; Couch & O'Connor 2014; Takiwaki et al. 2014), though the most important is likely the inverse turbulent energy cascade found in 2D simulations (Couch & Ott 2015 and references therein). Observations indicate that massive stars explode successfully as CCSNe (e.g., Smartt 2009). The difficulty in obtaining explosions for massive stars in 3D simulations may indicate that something crucial is missing in our theory of the CCSN mechanism. More realistic, 3D progenitor models will undoubtedly have an important impact on the CCSN mechanism and may be crucial to obtaining robust explosions.

We proceed with a discussion of our simulation approach in Section 2. We describe the 3D evolution during the final three minutes prior to core collapse of a massive star in Section 3. In Section 4, we examine the impact of 3D progenitor structure on the CCSN mechanism. Finally, we discuss the implications of our results and conclude in Section 5.

Our simulation of stellar core collapse proceeds in two steps. First, we evolve a non-magnetic, non-rotating $15\;{M}_{\odot }$ star in 1D with the Modules for Experiments in Stellar Astrophysics (MESA; version 6794; Paxton et al. 2010, 2013) to the point of iron core collapse, when the infall velocity of the outer core reaches ∼1000 km s⁻¹. At a point prior to collapse, during quasi-hydrostatic Si shell burning when the iron core mass is around 1.3 ${M}_{\odot }$, we map the 1D MESA model into 3D using the FLASH simulation framework (Fryxell et al. 2000; Dubey et al. 2009).

For the MESA models,¹¹ we use the "Helmholtz" (EOS; Timmes & Swesty 2000). The nuclear reaction network was automatically extended during the evolution starting from a basic 8-isotope network and reaching a 21-isotope network ("approx21") by the end of the calculation. Standard mass-loss prescriptions were adopted (Vink et al. 2001; Glebbeek et al. 2009). We use the Ledoux criterion for convection including semi-convection, thermohaline mixing, and overshoot. We use a resolution parameter for the MESA calculations of 0.6, corresponding to 3646 zones in 1D at the time of mapping to 3D.

For the 3D FLASH simulation, we also employ the Helmholtz EOS and the same 21-isotope network, newly implemented in FLASH from the standalone public network.¹² The approx21 network includes a very approximate treatment of heavy element neutronization important for the near-collapse isotopic evolution of the iron core. Ideally, one would use a network with sufficient number of isotopes (∼60–100) to treat the core neutronization and URCA cooling directly. For simplicity, speed, and in order to maintain direct network equality between the 1D MESA models and the 3D FLASH simulations we use the reduced network.

We include the complete iron core in the 3D domain. The tools of stellar evolution model convective burning in essentially an average sense, neglecting multidimensionality and highly dynamic fluctuations. This gives rise to a practical problem for multidimensional simulation of stars. Convective progenitor models do not have consistent 3D turbulent flow, and when they are mapped onto a 3D grid as an initial state for hydrodynamic simulations, they always pass through a transient state during which a turbulent flow develops. This in not a problem of mathematics (the mapping) but one of physics (the turbulence). The 1D stellar models cannot provide the 3D information for the turbulent fluctuations (both amplitudes and phases). There are no known successful attempts to fake it; transients happen. Ideally the mapping from 1D to 3D should occur earlier in the evolution such that any initial transients are negligibly small.

In order to quell such initial transients due to the mapping from 1D to 3D, which typically manifest as strong radial waves as the star settles onto the new domain, we employ the hydrostatic initialization approach of Zingale et al. (2002) rather than initial damping (e.g., Arnett 1994). In this approach, the density profile is adjusted slightly while keeping the pressure profile fixed such that the equation of hydrostatic equilibrium is satisfied exactly throughout the model. This procedure is closed using the EOS. Nevertheless, some initial transient motion persists during the first ∼20 s of the 3D evolution. This erroneous expansion can effectively "de-evolve" the star, pushing the stellar structure into a state resembling earlier phases of the evolution. This results in the requirement for a prohibitively long period of simulation time to reach collapse. In order to overcome this lack of realism in the initial model, we enhance the parameterized rate of the neutronization reaction that converts ⁵⁶Fe to ⁵⁶Cr, thus enhancing the cooling of the inner part of the iron core to compensate for the inefficient rate of outward energy transport. The rate of this reaction was enhanced above the fiducial value by a factor of 50. This approach, though admittedly non-ideal, effectively damps the initial core expansion allowing us to reach iron core collapse in 3D in a reasonable amount of simulation time.

For the 3D FLASH simulation, we use Cartesian coordinates coupled with adaptive mesh refinement (AMR) in one octant of the full 3D sphere. The finest grid spacing, using 8 levels of refinement, is 16 km. The refinement level is reduced as a function of radius, with the first reduction occurring around a spherical radius of 2500 km, beyond the Si-burning shell. The entire domain is 100,000 km on a side. Along the octant symmetry planes we use reflecting boundary conditions while at the outer extents of the domain we use zero-gradient boundary conditions. For solving the hydrodynamics, we use directionally unsplit PPM (without contact steepening) and the HLLC Riemann solver. The unsplit PPM solver in FLASH does not utilize the "consistant multifluid advection" scheme of Plewa & Müller (1999). We include self-gravity assuming a spherically symmetric (monopole) gravitational potential. Due to the operator splitting between the hydrodynamics and the nuclear network, we find that, in order to maintain adequate coupling between the burning and the hydrodynamics, it is critical to limit the size of the time step so that the internal energy in any one zone changes by no more than 1% during the course of a single step.

We follow the violent Si shell burning and build up of the iron core mass in 3D for ∼160 s. This single simulation required approximately 350,000 core-hours on Stampede at TACC on 1024 cores. Visualizations of the 3D progenitor burning simulation are shown in Figure 1, where we show slice plots of key quantities along with volume renderings of the iron core and the radial component of the velocity in the Si-burning shell.¹³ In Figure 2, we show spherically averaged radial profiles of the density, ρ, electron fraction, Y_e, specific entropy, s, and convective velocity from the 3D FLASH simulation at three times: the transition from the 1D MESA model to 3D; 5 s prior to collapse; and the point of collapse. Angle averages, $\langle ...\rangle$, are taken over spherical shells and v_r is the radial velocity component. Evident from Figure 2 is that the convective speeds near collapse are typically >100 km s⁻¹. This is slightly larger than the comparable speeds found in the O-burning shell (Viallet et al. 2013). Because the nuclear burning is balanced on average by turbulent dissipation, the energy generation rate is related to the average velocity and the depth of the convection zone by $\epsilon \sim {v}^{3}/{\ell }$ (Arnett et al. 2009). This is, however, the average convective speed, and fluctuations increase the peak speeds. We see from Figure 1 that the peak speeds in the Si-burning shell can be several hundred km s⁻¹, reaching speeds near collapse of ∼500 km s⁻¹. This is not negligible when compared to nominal infall speeds for core collapse initial models (∼1000 km s⁻¹). The speed of the convection increases as collapse approaches and the core contracts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 2 we also show final 1D MESA models at the point of collapse considering two different scenarios: one in which the neutronization reaction rate is enhanced by the same amount as in the 3D FLASH simulation (blue dashed lines) and the other in which we do not enhance this reaction rate above the fiducial value found in MESA's approx21 network (cyan dashed lines). Stellar collapse is highly dynamic and the model profiles change rapidly once gravitational instability sets in. Thus, for the sake of fair comparison, we consider all models at the point when the central densities have reached the same value as that of the final 3D FLASH simulation ($7\times {10}^{9}$ g cm⁻³). Varying the rate of neutronization has very little effect on the 1D MESA profiles at the point of collapse. It does, however, dramatically change the time it takes to reach collapse. For the fiducial neutronization rate, the MESA model takes 1000 s to reach collapse from the model we used for the 3D ICs (i.e., the black lines in Figure 2). With the neutronization rate enhanced as we have done for the 3D FLASH simulation, the MESA model requires just 40 s to reach collapse. This is faster even than the FLASH simulation because the MESA model, obviously, does not experience the initial transient pulsation that the 3D simulation must go through until vigorous convection is established. The final iron core masses of all models are similar, though correlated with the time it takes to reach collapse: 1.46 M_☉ for the enhanced-rate MESA model, 1.51 M_☉ for the fiducial-rate MESA model, and 1.50 M_☉ for the 3D FLASH simulation.

Figure 3 shows the time evolution of the central density, central electron fraction, the mass of the iron core, the total kinetic energy, and the net nuclear energy release rate in the Si-burning shell for the 3D FLASH simulation. We define the Si-burning shell as the region of the star where both the iron and silicon mass fractions exceed 0.001. During the simulated time, the iron core mass grows from 1.3 M_☉ to 1.5 M_☉, while the central Y_e decreases from 0.432 to ∼0.428 at collapse. The initial transient wave during the first 20 s is evident as an excess of radial kinetic energy in the bottom panel of Figure 3. After the brief transient, the kinetic energy is dominated by convective motion in the Si-burning shell until the collapse begins, around 160 s. The growth in the kinetic energy corresponds directly with a positive increase in the net nuclear energy release rate, which is negative because it is dominated by strong cooling from the central core. We expect that the convective kinetic energy is not unduly influenced by initial transient behavior because it reaches a near steady-state long before the end of the 3D simulation, and because we have simulated several convective eddy turnover times. Assuming an average convective speed of 100 km s⁻¹, and a width of the Si-burning convective region of ∼1000 km (see Figure 2), the turnover time is approximately ${\rm{\Delta }}{t}_{\mathrm{eddy}}\sim 2{\rm{\Delta }}r/{v}_{\mathrm{con}}\sim 20$ s. Driven convection/turbulence reaches a quasi-steady-state on roughly only a few turnover times (e.g., Meakin & Arnett 2007; Radice et al. 2015).

Zoom In Zoom Out Reset image size

Figure 1 also shows the initiation of convection from the initially spherically symmetric model (right-most panel of top row). As in the 2D simulations of Arnett & Meakin (2011), the Si shell burning is inhomogeneous with alternating regions of net cooling by neutrinos and net heating by nuclear burning. The iron core shows overall contraction while the Si shell is depleted by the burning. The bottom row of Figure 1 shows that, near collapse, the surface of the iron core is significantly distorted from spherical symmetry. At this time, the average-weighted standard deviation of the radial coordinate of the surface of the iron core is 9.5%. The radial velocity also shows that the velocity fluctuations are large in both scale and amplitude. We quantify the spatial scale (as well as strength) of the turbulent convection in the Si shell by measuring the turbulent kinetic energy power spectrum in spherical harmonic basis (for details on how this is computed, see Couch & O'Connor 2014). For the 3D progenitor simulation near collapse, this is shown in the top panel of Figure 4. The turbulent energy spectra peak at quite small ${\ell }$ (large scale), around ${\ell }\approx 4$. This implies that the largest eddies, which carry the bulk of the convective/turbulent kinetic energy, have roughly the same radial extent as the convective shell, precisely as observed in Figure 1. We have also analyzed the convection using the vector spherical harmonic framework of Chatzopoulos et al. (2014) and find that the vector spectra also peak around ${\ell }\approx 5$.

Zoom In Zoom Out Reset image size

The final stages of nuclear burning in massive stars can excite waves that could have important implications for mass loss and angular momentum transport just prior to core collapse (Meakin & Arnett 2007; Quataert & Shiode 2012; Shiode & Quataert 2014; Fuller et al. 2015). We find that the Si shell burning in the minutes before collapse drives non-spherical gravity waves that propagate both outward away from the core and inward into the iron core. Previous multidimensional studies of late-stage stellar burning (e.g., Bazan & Arnett 1994; Meakin & Arnett 2007; Arnett & Meakin 2011) saw such phenomena but lacked fully dynamic inner cores and were not able to assess them accurately. Such waves persist throughout the Si burning phase and are present at the point of collapse, as evident in the spherically averaged convective velocity shown in Figure 2. The wave velocity amplitudes at the start of collapse are tens of km s⁻¹. The contraction, and ultimately unstable collapse, of the iron core amplifies these waves (Lai & Goldreich 2000) and they become quite substantial around the moment of core bounce.

Our primary interest in simulating the final stages of stellar evolution in 3D is to assess whether realistic multidimensional progenitor structure has a significant impact on the CCSN mechanism. In order to achieve the 3D simulation of the final minutes of a massive star's life up to the point of core collapse we had to make certain approximations: (1) use of a simplified nuclear reaction network and adjusted electron capture rates, (2) simulating only the final moments of evolution after O shell burning has ended (for this model), and (3) use of a 1D MLT initial model. In this section we present a brief examination of the impact of such idealized but fully 3D ICs on the CCSN mechanism.

We carry out two 3D Newtonian simulations of collapse, bounce, and shock-revival including approximate neutrino physics. For one simulation, we use the full 3D progenitor simulation continued from the moment of collapse $({t}_{3{\rm{D}}}\approx 160\;{\rm{s}})$. In the other, we angle-average the final 3D progenitor model, washing away all non-radial velocities and non-spherical structure, to produce 1D ICs which are used in an otherwise fully 3D simulation. Our simulation approach is essentially identical to that of Couch & O'Connor (2014) and Couch & Ott (2015), with the only major difference being that our 3D simulations here are carried out only in one octant, rather than the full star. We use Cartesian coordinates with AMR, yielding a finest grid spacing of ∼0.5 km and an effective angular resolution of ∼05. Neutrino effects are incorporated using a multispecies leakage scheme that includes charged current heating and pre-bounce deleptonization is approximated using the density-dependent approach of Liebendörfer (2005). We enhance the local charged current heating rates by 6% in both cases in order to yield explosions.

Figure 5 summarizes our investigation of the impact of 3D ICs on the CCSN mechanism. The top panel of Figure 5 shows the average shock radius of the full 3D ICs simulation (red lines) and the angle-averaged 1D ICs simulation (blue lines). The average shock radii between the two simulations is nearly identical until around 100 ms after bounce. This corresponds with the time at which the Si shell interface is accreted through the shock surface. At this time, the shock in the full 3D ICs model begins to expand more rapidly than for the 1D ICs model, behavior that continues as both models transition to explosion. The second panel of Figure 5 shows the neutrino heating efficiency, ${\eta }_{\mathrm{heat}}={Q}_{\mathrm{net}}{({L}_{{\nu }_{{\rm{e}}},\mathrm{gain}}+{L}_{{\bar{\nu }}_{{\rm{e}}},\mathrm{gain}})}^{-1}$, where ${Q}_{\mathrm{net}}$ is the net charged current heating rate in the gain region, which is divided by the sum of the electron-type neutrino and antineutrino luminosities at the base of the gain layer. We find essentially no difference in the neutrino heating efficiencies between the two simulations. The same can be said also for the neutrino luminosities and the total gain-region heating rates. Clearly, the cause for the divergence in the results between the 3D ICs and the 1D ICs is not due principally to differences in the neutrino heating or average gain region matter dwell times (e.g., Murphy & Burrows 2008; Murphy & Meakin 2011; Dolence et al. 2013). Our present results are consistent with the picture in which the differences are due primarily to different strengths of post-shock turbulence, and the attendant effective turbulent pressures that aid shock expansion (Couch & Ott 2015). The third panel of Figure 5 shows the total turbulent kinetic energies in the gain region for both simulations. The 3D ICs simulation has greater turbulent energy after about 100 ms, following accretion of the Si interface. There is also a larger "burst" of turbulent energy immediately post-bounce caused by the presence of aspherical, convection-generated waves in the inner part of the core right at bounce.

Zoom In Zoom Out Reset image size

We also find differences in the turbulent kinetic energy spectra between the 3D ICs simulation and the 1D ICs simulation. The bottom two panels of Figure 4 show the turbulent energy spectra from the two collapse simulations at different radii: ∼300 km (ahead of the shock) and ∼150 km (in the gain region). The spectrum for the 1D ICs simulation for the region ahead of the shock, which should be essentially a spherically symmetric accretion flow, quantifies the magnitude of the perturbations excited by our use of a Cartesian grid. The turbulent energy in this region for the full 3D ICs is an order of magnitude or more greater than this at all scales (except ${\ell }=4$). This translates into greater turbulent energy in the gain region, as seen in the bottom panel of Figure 4. The 3D ICs result in greater turbulent energy on large scales, between about ${\ell }=6-10$, precisely where it is most effective at aiding shock expansion (Hanke et al. 2012; Couch & O'Connor 2014; Couch & Ott 2015). The excess of power in the quadrupole ${\ell }=2$ mode is likely due to our use of an octant domain.

The greater strength of turbulence excited by realistic 3D ICs also results in a greater diagnostic explosion energy, as seen in the bottom panel of Figure 5. Here, the diagnostic explosion energy, ${E}_{\mathrm{exp}}$, is the total energy of all gravitationally unbound material with net positive radial velocity (see definitions in, e.g., Müller et al. 2012; Bruenn et al. 2014). We note, however, that at the end of our simulations the explosion energies are not near their asymptotic final values, so caution should be used when interpreting the differences in the diagnostic explosion energies.

We have carried out the first 3D simulation of the final growth of the iron core in a massive star up to the point of gravitational instability and collapse. The violent Si burning in the shell surrounding the iron core drives large-scale, strong deviations from spherical symmetry that is not captured by 1D models. We show that this has a positive impact on the favorability for explosion via the delayed neutrino heating mechanism for CCSNe. We were forced to make a number of approximations, e.g., the use of a much reduced nuclear network, modifying electron capture rates, using MLT 1D initial models, simulating only one octant of the full star, and focusing on only one initial progenitor model. Thus, crucial aspects of the final iron core, such as the electron fraction and entropy, may be affected. Additionally, the presence of a Si-burning shell at the moment of core collapse may be progenitor model-dependent (see, e.g., Woosley et al. 2002).

Compared to 2D simulations of Si shell burning (Arnett & Meakin 2011), the maximum convective speeds we find are somewhat smaller, in agreement with 2D versus 3D results found for O burning shells (Meakin & Arnett 2007). The speeds we find are smaller, also, than those assumed in the parameterized velocity fluctuations employed by Couch & Ott (2013). As a result, the impact we find on the CCSN mechanism may be less dramatic than that found by Couch & Ott (2013), though it is still significant. We show that the biggest impact of 3D ICs is the enhancement of the strength of post-shock turbulence. Greater turbulence behind the shock leads to a greater effective turbulent pressure and, thus, more favorable conditions for shock expansion and explosion (Murphy et al. 2013; Couch & Ott 2015).

Realistic 3D simulations of the convective nuclear burning in massive stars, comparable to those used here, also suggest modifications to the MLT algorithms used in stellar evolution calculations (Arnett et al. 2015). Here, we have focused on the impact that 3D ICs have on the CCSN mechanism, showing that the breaking of spherical symmetry has a significant and positive impact on the likelihood for explosion. This is a critical finding since robust neutrino-driven explosions have been notoriously difficult to achieve across the broad range of (spherically symmetric) progenitors studied to date (see, e.g., Janka et al. 2012). The first 3D simulations including detailed neutrino physics indicate that explosion may be harder to attain than for 2D (Tamborra et al. 2014).

The CCSN "problem," i.e., the persistent difficulty in achieving robust explosions, has been with us for nearly a half century since the neutrino mechanism was introduced (Colgate & White 1966; Arnett 1968). Most theoretical investigations of the CCSN mechanism have, however, employed spherically symmetric, non-rotating, non-magnetic stars. It is well known that all stars rotate to some degree, are endowed with magnetic fields, and are spherically symmetric only in the average sense. We have shown here that the detailed 3D structure of massive stellar cores, that is a natural and unavoidable result of the final stages of nuclear burning, has a big and important impact on the CCSN mechanism. Additional investigation of the roles that rotation and magnetic fields play in both the pre-collapse evolution of massive stars and in the CCSN mechanism itself is needed, but taken together with realistic 3D progenitor structures, it is conceivable that attaining energetic explosions that reproduce the observable and statistical properties of CCSNe (Clausen et al. 2015) may not be so challenging. That is, there might not be a CCSN "problem" per se, we may have just been using progenitor models that do not occur in nature, and given a set of progenitor models more representative of real massive stars, the "problem" may go away.

This work is a first step toward addressing these issues. It serves to a large extent as a proof-of-principle indication that the final minutes of massive star evolution can and should be simulated in 3D. Doing so has a significant impact on the final structure of the core and on the CCSN mechanism. There is significant progress that must still be made in accurately simulating CCSN progenitors, including better models for mixing during the long, evolutionary phases of the star's life, and inclusion of realistic treatments for rotation and magnetism.

We thank J. Fuller, C. Graziani, C. Meakin, C. Ott, and M. Zingale for many valuable conversations. S.M.C. is supported by the National Science Foundation under award No. AST-1212170. E.C. would like to thank the Enrico Fermi Institute for its support. F.X.T. thanks NASA for partial support under TCAN award NNX14AB53G. This research was partially supported by the National Science Foundation under award No. AST-1107445 at the University of Arizona. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This research used computational resources at ALCF at ANL, which is supported by the Office of Science of the US Department of Energy under contract No. DE-AC02-06CH11357, and at TACC under NSF XSEDE allocation TG-PHY100033.