Skip to main content


Code Development

Even while we use our existing codes to do great science, our group is continually developing the next generation of numerical methods for evolving hydrodynamical systems.

Our group's numerical focus is on the idea of a "moving mesh". This is a numerical method first developed by Volker Springel and used to make the AREPO code, which is responsible for the famous "Illustris" simulation describing the formation of structure in the universe. Our group typically builds moving-mesh codes that exploit particular symmetries in the flow, in order to more accurately or more efficiently capture the fluid motion.

Check out our "Public Codes" page for codes that have already been developed by our group. Here are a few examples of numerical methods presently in development:

Sprout -- Hydrodynamics on a Homologously Expanding Grid

(details coming soon!)

Nuclear Hydrodynamics on a Moving Mesh

(details coming soon!)

Cube plot

Rayleigh-Taylor Instability in Supernovae

The James Webb Space Telescope has ushered in a new era of astrophysics, and we anticipate many breathtaking images of supernova remnants that other space telescopes have viewed before, but now with far greater resolution. This is important because we think there is a great deal of information living at these smaller angular scales.

Additionally, new transient searches like the Zwicky Transient Factory and LSST will observe a large number of core collapse supernovae (CCSNe). Light curves and measured ejecta velocities are often calculated using 1D models, for reasons of expediency; millions of 1D calculations can be performed in the months it can take to run a single end-to-end 3D calculation. On the other hand, some important effects necessitate a multidimensional treatment. One particular example our group has studied is the effects of 3D mixing by Rayleigh-Taylor instability (RTI) on Type IIp SN light curves.

Mixing is a crucial component to calculating SN light curves, and in absence of a good mixing model one must run a full 3D calculation, a single one of which can take months, in order to model a single light curve from first principles. After modeling these systems in 2D and 3D, we found a way to extract the essential physics of mixing found in 3D and to apply it to 1D models. This makes it possible to perform fast supernova calculations over large ranges of parameter space that would otherwise require expensive 3D calculations. The mixing model (now called "Duffell RTI") has now been incorporated into the hydrodynamics module of the MESA stellar evolution code so that the astrophysics community now has access to the physics extracted from our 3D models, and can apply this physics to their 1D models.

A broadly-defined project our group is pursuing is to characterize important 3D effects of turbulence from fluid instabilities in supernovae. It is known that the details of CCSNe depend on many 3D effects (not just RTI). By carefully studying these fluid instabilities in idealized scenarios and extracting fundamental scaling laws, it may become possible to massively simplify calculations that otherwise cost many millions of computing hours to execute. Our group is presently working on developing a novel characterization of supernova remnant anisotropy, in order to build a forensic tool for measuring the presence of different fluid instabilities and other 3D effects in different types of supernovae.

Circumbinary Disks

Stars are born when a molecular cloud suffers gravitational instability and collapses. What immediately forms is a protostar embedded in a protostellar disk. If a binary star system forms while the disk is still present, then the formative years of the binary occurred in the presence of an accreting circumbinary gas disk, whose presence could have a significant impact on the binary orbit and mass ratio.

Binary accretion is also fundamentally important for a different problem, on much larger scales: that of supermassive black holes accreting at the cores of galaxies. Each galaxy has a supermassive black hole at its core. When two galaxies merge, their supermassive black holes migrate towards the center of the merged galaxy, at which point one might expect a supermassive binary merger. Such mergers could be detectable in the gravitational wave background observed by pular timing arrrays (PTAs) such as NANOGrav and the international pulsar timing array (IPTA). However, all of this is contingent upon the true resolution of the "final parsec problem", or what mechanism is able to drive the binary toward merger within a Hubble time, once the binary has shrunk down to a separation of a few parsecs?

Gas accretion has been suggested as a likely mechanism (assuming that the gas can drive the binary inward), but if that is true then one expects gas accretion to leave an imprint on the observed spectrum of the gravitational wave background, as it would shorten the residence time of parsec-separation binaries, and therefore the amplitude of the signal at the corresponding binary frequencies.

Both of these very different systems (stars and supermassive black holes) exhibit the phenomenon of "binary accretion", which can be summarized as a fundamental physics problem: gas dynamics coupled to the time-dependent gravitational potential of an orbiting binary. This setup has been studied numerically over the past three decades, but the methodology, numerical resolution, and the duration of the calculations have improved greatly in the recent past.

We are presently collaborating with Daniel D'Orazio at the Niels Bohr Institute in Copenhagen, performing calculations of accretion onto eccentric binary systems. Our results will enable predictions for demographics of observed binary star systems, as well as predictions for the accretion variability. Additionally we are participating in the Santa Barbara Binary-Disk Code Comparison, in collaboration with over a dozen other researchers on this topic. We have developed the "minimal binary-accretion" setup and are running the same identical setup using eleven different numerical prescriptions, in order to determine the necessary resolution for all codee to agree on the answer to this problem.

Gamma Ray Bursts and Neutron Star Mergers

Some relativistic transients, such as gamma ray bursts (GRBs), have been observed for many years and excellent data from prompt and afterglow emission are now publicly available due to excellent space-based observations including Fermi, Swift, Chandra, and XMM-Newton. Very recently, there has been a resurgence in interest in GRBs coinciding with the discovery of coincident gravitational waves, gamma rays, and kilonova emission from gravitational wave source GW170817, caused by a merging neutron star binary system. These observations signal a rebirth of the GRB community, as the multi-messenger detection has made it possible to investigate the event in new ways.

GRBs are typically highly relativistic, with Lorentz factors in the hundreds or sometimes thousands (velocities close to the speed of light, within one part in 106). Numerically modeling any sufficiently relativistic gas dynamics system is difficult; most codes fail when the Lorentz factor is in the tens. For such relativistic outflows, we have developed the JET code. One of the great advances of the moving mesh is that it allows the evolution of ultra-relativistic flows, with Lorentz factors in the thousands.

One recent example utilizing this code was in modeling the dynamics of the different relativistic components of the electromagnetic counterpart to GW178017. We have performed calculations of the relativistic jet colliding with the mildly relativistic ejecta responsible for kilonova emission. By performing a large parameter survey of jet and ejecta models, it was possible to extrapolate rules for whether the jet successfully breaks out, how long the jet takes to push its way through the ejecta, and how much thermal energy it imparts in the process. This made it possible to rule out several theoretical models regarding the interaction, and we were able to conclude based on this theory that a successful jet had made it through the ejecta, a statement which has since been confirmed observationally. We also showed that the observed 1.7 second delay between gravitational waves and gamma rays could be explained by the time the jet takes to push through the ejecta, which predicts a relationship between the ejecta mass, the GRB energy, and the delay time, a relationship that will be tested in future multi-messenger events over the coming decade.

GW170817 also included a radio afterglow signal, which is still being observed today. The rise and fall of the afterglow light curve probes the fact that the jet is "structured", i.e. has a nontrivial distribution of energy with angle. Recently, the afterglow light curve was well-fit using an analytical model our group had developed for a structured jet called the "boosted fireball". Our model was taken as initial conditions for a numerical hydrodynamics and synchrotron radiation calculation, in order to determine parameters the outflow structure. This has allowed for a very accurate measurement of the outflow structure, and has made it possible to rule out the idea that the light curve might have been produced by a spherical outflow.

Protoplanetary Disks

Planets are born in a gaseous disk (PDF). It is important to understand the interaction between the planet and the disk for many reasons, both theoretical and observational. Recent observational advances have resulted in new resolved and detailed images of protoplanetary disks, revealing interesting features which could be the result of a planet interacting with the disk. Of particular interest are the so-called "transition disks", or disks which exhibit observational signatures of gaps or cavities. There have been many new, spatially resolved observations of such disks in the past decade. Not only does this provide a new means of discovering planets, but it is a way of discovering planets in the formation process. This strongly motivates a research program to investigate the physical relationship between these disk structures and the properties of the perturbing planet.

The first application of the DISCO code investigated the gap opening process for low-mass planets. We demonstrated that low-mass planets can open gaps, contrary to previously existing theory. It used to be thought that gap opening necessitated a gas giant, but using DISCO we have shown that super-Earths or even Earth-mass planets can open gaps in sufficiently low-viscosity disks. This resulted in a new gap-opening criterion that overturned previous models, and a formula for the gap depth as a function of planet mass, disk viscosity, and scale height which has since been confirmed in several numerical and analytical studies. Our gap model has been used to analyze three low-density rings seen in direct imaging of the disk TW Hydrae. Using our model it was possible for observers to conclude that the embedded planets were at most 10 M, and given such limits the gaps were likely too wide to be carved out by planets.

We have also performed several studies on planetary migration. We have shown that gap-opening planets do not migrate on the disk's viscous timescale, contrary to established theory. This has dramatic implications for planet formation theory as it might explain why planets do not migrate into their host stars after opening a gap. Population synthesis models of planetary systems have shown that reducing the migration rate of massive planets can help to explain the present demographics of planetary systems.