3-D Stellar Hydrodynamics
Overview. This project is a joint effort between the teams of Paul Woodward, at the LCSE at the University of Minnesota, and Falk Herwig, at the University of Victoria. Herwig brings to this collaboration his deep knowledge of stellar evolution and the nuclear reactions that are involved. He has worked with the MESA stellar evolution code extensively and is an expert in stellar nucleosynthesis. His interest in investigating convective boundary mixing in 3D in the special cases he could identify as critically important from his study of 1-D stellar evolution was an initial motivation to begin this collaboration. Woodward brings to this project his many years of simulating convection in 3D, as well as his extensive work on high performance parallel computing in astrophysical fluid dynamics. The joint work between the two research teams began in 2006, with a combination of support from the DoE, the NSF, and, after Herwig took his present professorship at the University of Victoria, from NSERC in Canada. Computing resources for Herwig now are provided principally on Niagara, in Toronto, and for Woodward are now provided principally on Frontera, at TACC in Texas. Both Woodward and Herwig participate in the program of the JINA-CEE physics frontier center, supported by NSF, where Herwig is a member of the executive committee.
More recently, it has been confirmed impressively with the most modern asteroseismology data from the Kepler satellite [1,2] that the best agreement between stellar models and observations can be found with the exponentially decaying diffusive model for CBM , that has been first introduced into stellar evolution models of AGB stars [4,5]. The asteroseismology observations determined the best-fitting e-folding distance for CBM at the top of the H-burning core convection in main-sequence stars to lie in the range f = 0.017 and 0.030, while recent studies of double-lined eclipsing binaries find f ~ 0.014 … 0.021 for main-sequence stars more massive than 2 Msun . This parametric CBM model has since been deployed by many groups in low-mass stellar evolution [7-13], and calibrated against numerous observational phenomena. For example, models with this exponential CBM model agree better with observations of O in planetary nebulae [14,15], can explain the C and O enhancements observed in novae  and can reproduce the abundances observed in H-deficient post-Asymptotic Giant Branch (AGB) stars [17,18]. CBM plays an important role as a key mixing process for the s-process in AGB stars . In the most massive AGB stars (super-AGB stars), CBM has been demonstrated to lead to convective-reactive mixing episodes that may trigger the exotic intermediate neutron-capture process (i-process) nucleosynthesis and possibly eruptive mass-ejections and transients triggered by convective-reactive H/He mixing of H-ingestion. At the bottom of the C-burning shell in super-AGB stars on the other hand, CBM was found to possibly quench the C-burning shell, which may lead to the formation of hybrid C-O-Ne white dwarfs [20,21], and the implications, caveats and possible consequences of this scenario are presently being investigated [22-26].
In massive stars CBM in the exponential, diffusive model has been adopted beyond core He burning without detailed analysis of its effects [27,28] while the possible role of CBM in the later phases of massive star evolution has been speculated about . Considering the long history of studying CBM in post-He core burning stellar evolution in low-mass stars, this macrophysics effect has only recently been gathering attention in massive stars. Our team has just completed a systematic study of the effect of CBM in post-He core burning shells in a 25 Msun stellar evolution model . When comparing a benchmark case without (or with very small) CBM with three cases of CBM efficiencies covering the range of f = 0.012 to 0.032, nonlinear and non-monotonic variations of the core masses by 80% and of the compactness parameter of a factor ~2 are found. Differences in yields of elements such as C, Ne, Mg and Si can reach 1dex. Such large variations mean that whether or not this particular model will explode as a SN or collapse into a black hole will depend on the adopted CBM efficiency. Considering its widespread adoption, including, for example, in the widely used public stellar evolution code Module for Experiments in Stellar Evolution (MESA) [31-33], it is troubling that we have only very limited information on the efficiency of CBM in this model directly from hydrodynamic simulations. However, such 3-D hydrodynamic simulations are the only first-principles approach to investigate CBM in its various stellar regimes. One-dimensional simulations cannot capture the turbulent 3-D nature of the hydrodynamic instabilities that cause CBM. Such 1-D simulations can also not capture the often global modes of convective, large-scale flow that drive flow topology at the convective boundary, as became clear from the global nature of convection in thick shells . Furthermore, 1-D simulations would be incapable of reproducing global instabilities, such as the Global Oscillation of Shell H ingestion (GOSH) that we discovered in 3-D simulations of convective-reactive flow in H-ingestion simulations in He-shell flash convection in post-AGB stars  and low-Z AGB stars .
For the reasons given above, we believe it is of great importance to understand CBM efficiencies across stellar evolution better, and ideally to develop a predictive capability of CBM from hydrodynamic simulations that would be based on first principles and validated against the numerous observational constraints. This would address and possibly remove one of the key gaps in our understanding of the evolution and final fate of stars, and in turn allow stellar physics to serve down-stream research applications, such as galactic chemical evolution, stellar populations, supernova research, and even gravitational wave astronomy, much better. While there are numerous hydrodynamic simulations that describe the phenomenon of CBM, overshooting, and penetration [3,37,38] or Kelvin-Helmholtz mixing , there are very few studies that attempt to determine the parameters of this CBM model that is implemented in many stellar evolution codes. For He-shell flash convection in AGB stars the f value has been estimated as fbot = 0.01 for the bottom and ftop = 0.10 at the top, based on 2-D plane-parallel stellar hydrodynamic simulations of marginal quality . A significantly improved determination of f has been recently presented by our team for the top of the O-burning shell convection in a 25 Msun model star through 3-D hydrodynamic simulations. These were performed by applying the methods and code technologies that we used in the past to study He-shell flash convection [34-36], and that we will further develop and apply to new problems as this project continues. Based on these high-fidelity simulations of the entire shell in 4 pi on a 15363-cell Cartesian grid, we performed a diffusion coefficient analysis of the evolution of the spherically averaged abundance profiles and found that f = 0.03 at the top of the shell. This value is interestingly in the same range as those observationally determined from asteroseismology for main-sequence stars (see above), indicating some degree of convergence between the 3-D hydrodynamic approach and constraints from observations. However, the correspondence at this point is qualitative at best. We cannot use main-sequence constraints for CBM in the shells of post-He-core-burning massive stars.
Returning one more time to the investigation in  of the effect of CBM on post-He-core-burning evolution in massive stars, it should be mentioned that shell mergers are an important mechanism of amplifying the effect of CBM in non-linear and non-monotonic ways, depending on the efficiency of CBM shell mergers. As an example, we show the merger between the C- and Ne-burning shells in a simulation with f = 0.032 (Figure 1). The comparison sequence with lower CBM efficiency does not show this particular shell merger (but may show another one). The importance of shell mergers in the progenitor evolution of core-collapse supernovae has recently been emphasized in several studies [29,41], and mergers may induce asymmetric large-scale fluctuations of the convective velocity spectrum that could help successful explosions [42,43]. However, shell mergers are not only important for the final structure of the collapsing star, which will determine its ultimate fate as a black hole or supernova. We have recently shown in work that combines our 3-D hydrodynamic simulation capability with our multi-zone nucleosynthesis simulations that convective-reactive nuclear processing in a C- and O-shell convection zone merger would lead to an exclusive overproduction of odd-Z elements P, Cl, K and Sc . Interestingly, it is just these odd-Z elements that are under-produced in present galactic chemical evolution simulations. We show that if we assume that such C-O shell mergers occur in about 50% of all massive stars, which is a rate much higher compared to 1-D stellar evolution predictions, then this conundrum of nuclear astrophysics can be naturally resolved and observations of galactic chemical evolution be accounted for.
The above result for C-O shell mergers is reminiscent of the convective-reactive event that is found when H is entrained into He-shell flash convection in a post-AGB star  or in a low-Z AGB star . This situation gives rise to the i-process, or intermediate neutron capture process , and has more recently been associated with highly anomalous abundance features in C-enhanced metal poor (CEMP) stars that formed in the early universe. Specifically, the CEMP-r/s stars, that were previously thought to originate from a superposition of slow and rapid neutron-capture process agree very well with predictions of simplified i-process simulations [48,49]. However, none of these exploratory, and overly simplified simulations of i-process nucleosynthesis takes into account the recent finding that the H-ingestion flash has a highly non-radial and violent characteristic in 3-D simulations, where a Global Oscillation of Shell H-ingestion (GOSH) has been discovered . Likewise, the initial simulations of nucleosynthesis of mergers of O- and C-shells were based on stellar evolution simulations. However, our most recent work  combines nucleosynthesis simulations in 1D with convective mixing profiles that were derived from 3-D hydrodynamic simulations for conditions similar to those expected in full-blown mergers . These simulations were performed on Blue Waters and investigate the energetic impact of C ingestion into an O-shell convection zone, using two different nuclear network approaches to account for the complex nuclear burning in combined O and C burning.
The project we are describing here has 3 main objectives:
Objective 1, validating a useful 1-D model of CBM.
In , we focused on ingestion of stably stratified gas from above the convection zone driven by an O-burning shell in a 25 Msun star. This case is similar to one that had been studied earlier by Meakin and Arnett . In this case of the massive star, the resistance of the stably stratified gas to being pulled into the convection zone was less, and we were able to establish that two simulations, on grids of 7683 and 15363 cells, agreed closely on the value of the entrainment rate. In this study, we did not include the complicating positive feedback from burning of the ingested gas. However, further scrutiny of the determination of the f parameter for the CBM model (Fig. 22 in ) has motivated further investigation on how well the CBM process is determined. In 2017, we have therefore extended this study of entrainment from above the O-shell convection zone by performing new 3-D runs at varying luminosities and high grid resolution, and analyzing the combined series of previous mostly lower-resolution runs and the new runs. Two of these runs can be seen in Fig. 3, with the volume fraction of ingested gas in a slice through the center of the star, volume rendered so that the ingested concentration decreases as the colors go from red, to yellow, white, aqua, and finally to dark blue. The color ranges are set to emphasize the similar character of the entrainment despite the two simulations having driving luminosities differing by a factor of 25.
In Fig. 4, fits of 1-D diffusion coefficients of the 3-D simulation data from the two runs shown in Fig. 3 are given. The method to determine the diffusion coefficient is the same as in . Spherically averaged radial profiles of the mixing fraction of entrained gas are shown in each panel in Fig. 4 at two times. From these two profiles, we determine the effective diffusion coefficient, D, that could reproduce this change in a 1-D computation. We constructed special runs, with sinusoidally varying passively advected tracers, in order to measure accurately the low values of the diffusion coefficients in these cases. Not only do these two runs show similar character of the entrainment in Fig. 3, but their two effective diffusion profiles in Fig. 4 are also similar. This encourages us to expect that a single parameterized form of the effective diffusion coefficient can be used in a wide variety of circumstances that could occur in our studies of convective boundary mixing and shell mergers in massive stars. In , we fit results for the case at the left in Fig. 4 to a specific form for the diffusion coefficient that was introduced by Freytag et al. in  and has been widely used in 1-D stellar evolution studies. The results of this fitting procedure are shown in Fig. 5. We can think of the 1-D model as having two especially important parameters. The first is an overall scale factor that determines the size of the diffusion coefficient well inside the convection zone. The second is f, the fraction of a local pressure scale height over which the exponentially decaying diffusion coefficient D diminishes by one factor of e.
In Fig. 6, fits to the results of the complete series of 3-D simulations for this entrainment into the O-shell convection zone are shown. The parameter f of the 1-D model is plotted against the log of the driving luminosity. The value of f relates to the thickness of the mixing region in the spherically averaged radial profile. The results in Fig. 6 show an expected increase of f with luminosity as the entrainment grows more vigorous. There is also some indication that values of f agree more closely between runs at our two grid resolutions as the luminosity increases. However, at the lower luminosities in the figure, it appears that, although the entrainment rates agree to about 10%, the values of f are further apart. It could be the case that f is hitting a floor value at low luminosities that is related to the grid resolution. The two runs in the diagram for the high-resolution grid and at the far left and far right in the plot are the two runs shown in Fig. 3, with effective diffusion coefficient profiles shown in Fig. 4. As can be seen most clearly in Fig. 5, the thickness of the radial profile of the mixing fraction in the run at the left in Figs. 3 and 4, which corresponds to the blue point at the bottom left in Fig. 6 is about 0.15 Mm. The grid in this problem runs out to a radius of 9.8 Mm, so that on the 15363 grid, 0.15 Mm is 11.76 cell widths. Given the sub-cell resolution of our PPB advection scheme describing the mixing fraction distribution (see discussion in ), this average thickness of roughly 12 cells is likely to be sufficient, or at least very nearly so, in our higher-resolution runs for this problem. An average thickness of 6 cells in the half-resolution runs is more likely to be marginal.
Objective 2, assessing the conditions under which two convective shells in close proximity can merge.
In our work on this objective, we will accumulate detailed data from multiple 3-D simulations of CBM which we will make available to the research community, as described in a later section. This data will enable the formulation and testing of more elaborate 1-D models by ourselves or by others. For example, it could be useful in computing energy release and positive feedback on CBM from burning of ingested fuel to build into a 1-D model some measure of the variance of the concentration of this fuel at each radius. This would reveal limited but important information beyond the simple average concentration. Such information could easily be mined from the data our study will make available.
As we have found through a detailed stellar evolution based investigation, the nature of shell mergers and CBM will also be key in determining the fate of the first generations of massive stars. One of the goals of our stellar hydrodynamics simulation program remains to address the convective-reactive phases in massive Pop III stars. While we are just beginning to make progress in actual 3-D hydrodynamic simulations, collaborator Herwig and his team have studied the initial conditions and the stellar evolution environment of this situation . We found that the abundance pattern observed in metal-poor stars involving high observed abundance ratios of Na/Mg or Ca can be obtained in i-process conditions. Such conditions may be realized in the ultra-metal poor star models where a He-shell convection zone and a H-burning shell convection just above interact. The new analysis shows that the initial situation in massive Pop III stars is much like two interacting convection zone boundaries in close proximity, as in the O-C shell case, rather than the ingestion of H from a stable layer above a He-burning shell convection zone encountered in low-Z AGB stars or post-AGB He-shell flashes with H ingestion. Again, the initial question of when H-He convection zone interactions and possibly mergers happen would start with a general investigation of when and how close convective shells interact.
We will approach this objective by a thorough analysis of existing massive star stellar evolution models, such as the grid of massive star simulations as a function of mass and metallicity by the NuGrid collaboration , those from [30, 53] and possibly . In this large data base of stellar models we will characterize the conditions for nearly-missed shell mergers, as well as marginal and actual shell mergers. We will consider H/He mergers, He/C, C/O and C/O/Ne mergers. We will determine typical entropy profiles. Cases where the entropy barrier between the two shells is smaller are expected to be more likely to merge, while those with either steep entropy gradients between the shells, large entropy differences, or larger spatial separation are cases in which mergers are less likely. Based on this analysis we will categorize entropy profiles into different classes according to some yet to be developed metric of “mergeability” of convective shells. We will then perform numerical experiments that are close to conditions encountered in stellar evolution models, yet general and representative enough to not be tied to one particular situation. This approach ensures that we are investigating general processes and not a very peculiar situation that only occurs in a particular case.
We plan to settle on a relatively small, but representative number of cases, and then conduct simulations for some parameter space that allows us to reveal systematics that can be generalized. We plan to include both spherical, full 4 pi simulations, as well as zoomed-in plane-parallel simulations that allow us to focus a lot of grid resolution at the interface between the two zones at the expense of less realistic and parameterized driving of the flow at the boundary. The result will be a suite of simulations that explore dynamic, convective mixing in the face of a positive entropy gradient, a novel approach in stellar hydrodynamics that will allow us to expand the understanding of this important macroscopic stellar physics process. The ultimate goal is to arrive at a model that encapsulates conditions for mergers, and allows an improved treatment of mergers in global stellar evolution simulations.
Objective 3, simulating the consequences of a shell merger.
In order to address this objective, we will select suitable setups and cases identified in step 2 and will perform sufficiently high-resolution 3-D hydrodynamic simulations. Most of the work related to this objective will be done in the second half of the project, because it can greatly benefit from using the AMR capability that is under development for PPMstar. We will extract important information about how mixing and the resulting energy generation in shell mergers will feed back into the fluid flow, and possibly induce even more violent flows. We have seen evidence for this when we studied the H-ingestion in low-mass stars, where the GOSH amplifies the mass entrainment rate by a factor of ~100, thereby sustaining the instability. We will use this information from 3-D hydro simulations to inform nucleosynthesis simulations. We have adopted this approach previously [44,71] by (1) combining stellar evolution simulations to identify the initial stellar conditions, then (2) carrying out 3-D hydrodynamic simulations to learn about the character of mixing in the brief, dynamical time-scale, convective-reactive episode, and then (3) calculating the detailed nucleosynthesis using the state-of-the-art NuGrid 1-D multi-zone nucleosynthesis code, with mixing assumptions that are either directly derived from (see, for example, ), or at least informed by (see, for example,), the 3-D hydrodynamical simulation.
In this project we plan to take this integration one step further. One way of integrating 3-D hydro simulations with stellar evolution is the autocorrecting modeling approach described below. Another approach is geared to take a realistic nucleosynthesis network to the 3-D simulations. This would be especially important in very non-radial and violent cases of shell mergers or the H-ingestion simulations with GOSH in low-mass stars. We have developed a concept of how to integrate networks of sufficient size (128, 256 or 512 species) into the 3-D simulations, either via a co-processing or a post-processing approach. The key idea of our approach is the realization that we are dealing in convective-reactive flows with distributed burning that only requires a reduced spatial resolution to still give a good approximation of the spatial energy distribution, and that we only have to follow the nucleosynthesis on the time scale of the fluid flow, but not of the sound speed. Essentially, this is equivalent to an implicit approach for the nucleosynthesis at reduced spatial resolution, the reduced cost of which balances completely the added cost of including several hundred isotopes. We have started with the implementation of this concept. For example, we have completed an initial draft of a nucleosynthesis module to be used in this new capability. We have also planned out in detail the advection part and the data management aspects of this approach, and we plan to continue the development of this capability in the new code framework, with the goal to perform 3-D nucleosynthesis simulations on top of our convective-reactive hydro shell mergers toward the end of the project.
An important outcome of investigations addressing this objective will be increased understanding of the post-shell merger state of the star. If these shell mergers are anything like the violent and eruptive GOSH H-ingestion events we found in low-mass stars, then the 1-D spherically averaged state after such events is unknowable from 1-D stellar evolution. The goal, within the framework of our autocorrecting modeling approach, would be to determine a spherically averaged post-merger state in the best possible way, and then to continue with the stellar evolution simulation. If the first such shell merger is that of a H- and He-burning shell in a Pop III or ultra-metal poor massive star, the 1-D simulations may later approach again a situation in which shell mergers could happen, such as a C-O shell merger. We would then repeat the process, studying that next shell merger and determining its final 1-D spherically averaged state. Along the way, we would determine to our best ability the CBM efficiency using the approach adopted in objective 1. Overall, the outcome of this process will be a complete massive star simulation all the way to collapse with the advanced guidance on dynamic mixing processes in shell mergers and at convective boundaries from high-fidelity 3-D hydro simulations. We have established collaborations (NuGrid, JINA-CEE) through which we can ensure that the supernova explosion properties of such progenitor models can be investigated (e.g. Chris Fryer [LANL], Sean Couch [MSU]). This goal of following a massive star all the way to collapse will be enabled through the autocorrecting modeling approach.
References cited in the above text are listed in the PDF document at www.lcse.umn.edu/content-research-refs.pdf.