[Beyond] Code in Astro: how to simulate a Protoplanetary Disc with DustPy

Paper Title: DustPy: A Python Package for Dust Evolution in Protoplanetary Disks

Authors: S. M. Stammler, Ti. Birnstiel

First-author institution: University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679, Munich, Germany

Status: Published in The Astrophysical Journal (Open Access)

Have you ever wondered what it’s like to do astrophysics research? Study the environment in which the Solar System might have formed? Model a protoplanetary disc? Planet formation is a vastly complex astrophysics problem, requiring advanced methods in both observations and theory. But at the core of almost all modern astrophysics research is code and models. Today’s Bite will walk you through the idea behind the state-of-the-art code DustPy, commonly used to model the evolution of dust and gas in a protoplanetary disc. 

If you’re new to simulations, I strongly recommend reading this previous Bite that explains the basics and jargon of simulations!

What’s a Simulation?

Complex systems require a careful approach. If you want to model the gravitational interaction between the Sun and Earth, you can solve the equations of motion for that system – the classic two body problem. For some simple set-ups like this, you can often write an analytical solution (i.e. equations that tell you the exact motion of the bodies). When you add a third body, like Jupiter, however, things become more complicated and you can’t always write an equation describing the bodies’ motions. When systems get complicated like this, we turn to numerical, computer-based simulations.

Simulating a Protoplanetary Disc

Protoplanetary discs are complex systems – they’re made of billions to trillions of dust and gas particles of different sizes, masses, temperatures, positions, and velocities. Unsurprisingly, it’s impossible to find an analytic solution for all these particles, so numerical simulations are required. Even then, we still need some simplifying maths – there’s no way we can model so much information at once! A common approach is to model the dust and gas as a fluid (i.e. you can model it like a smooth-moving liquid). This way, you can simulate the bulk motion of the dust and gas instead of individual particles, saving on computational costs. This is the idea of hydrodynamical (fluid-based) codes like Athena++ and PLUTO, where your quantities of interest behave like a fluid; check out this Astrobite for a more detailed view on how different simulation types work (grid-based, particle-based, [smoothed-particle] hydrodynamical, etc.).

So, we have a fluid that represents the dust and the gas. For some models, this simplification is sufficient, although others like to add additional physics on top of the motion of the dust and gas, like the growth of dust. Dust in protoplanetary discs should grow from the size of sand grains (<1mm) all the way up to planets ($\sim10^{4-5}$ km) – a size range that can span up to 13 orders of magnitude. My head hurts just thinking about this…

We need to reduce the simulation’s computational cost further, which we can do by reducing the number of dimensions in which we simulate the protoplanetary disc. It’s quite common to approximate protoplanetary discs using just one dimension – the radial distance from the star. This assumes that the disc is azimuthally symmetric around $\phi$, and that the vertical ($z$) component is in hydrostatic equilibrium. So, going to 1D basically takes a single radial slice of the disc whilst describing all of the vertical layers at once.

In 1D, we measure the dust and gas density in “surface density” $\Sigma$, measured in g cm$^{-2}$ (not g cm$^{-1}$, since we’re taking a radial slice!). Think of it like looking down on the surface of a protoplanetary disc: how much mass do you see in a column beneath a square centimetre? I’ve illustrated this in Figure 1.

At this stage, we have to pick how we want to model the dust growth. Dust growth is an extremely complex topic (people have written entire PhD theses on it!) and a lot of studies investigate dust growth explicitly. But, we can cheat directly modelling the growth of dust particles using sub-grid physics.

Sub-Grid Physics

This is where we model physics that is below the resolution of our numerical grid with other, sometimes analytical, models. The smallest spatial cell in our simulation might be 0.05 astronomical units (AU), but this is orders of magnitude larger than the biggest of dust particles, even as they grow. We therefore need another approximation: we will model the population of dust instead of simulating every particle per cubic centimetre. A population describes how the total number of particles of a given size per cubic centimetre evolves, rather than individual particles. This is commonly done with the Smoluchowski equation (see Figure 2).

This equation tells us how populations of dust evolve from growth and also fragmentation – it’s possible that two colliding dust particles won’t stick together, but instead fragment into smaller pieces! The equation is pretty monstrous, and requires clever tricks and computational power to solve.

That was only the assumptions?! This is getting complicated…

Don’t worry, we’re basically there! We have all our simplifying assumptions at hand: a 1D grid, modelling dust and gas as a fluid, and sub-grid physics to handle our dust growth. We’re now ready to simulate our protoplanetary disc, we just have to make a model.

Instead of leaving this as an exercise to the reader, we’ll pick a model that already exists: DustPy. This simulator does everything we’ve described above: gas and dust evolution, including dust growth and fragmentation, all in 1D. If you want the gory details on how it works, you can find them in the original paper describing the model in detail, and the newer DustPy paper.

DustPy tells us how the population of dust evolves over time; you can see an example in Figure 3 How many big and small particles are there? How does the dust move? How quickly does it grow and fragment? It provides an extremely detailed study of dust in protoplanetary discs, which form the building blocks of planets. Simulations like DustPy are crucial in studying how protoplanetary discs can be engines for planet formation.

If you would like to try DustPy, there is a thorough tutorial on the documentation page you can follow. You set up and run the model in Python, so it’s easy to work with. Note you need a Fortran compiler on your computer to run it. Simulations can also take a while if you run them for more than 500,000 years in simulation time! The basic tutorial can also take a little while to run (at least 10 minutes, depending on your computer). There are additional tutorials to extend the model, like including planets.

Okay, so how do we make a planet?

Congratulations, you now know how to simulate a protoplanetary disc! This is all quite complicated, and people spend whole careers working on simulations of protoplanetary discs. Hopefully this Astrobite has given you a taste of what a state-of-the-art model does.

Finally, it’s important to note that there’s a limit to how big dust particles can get. At some point, they’ll start to fragment into smaller particles instead of growing, or they might be transported away by the gas before they can get bigger. This size is often a few centimetres, which is nowhere near planet-sized! This problem is often called the “metre-sized barrier”. People dedicate whole papers, PhD thesis, and careers to study how to overcome this problem. Maybe you could study protoplanetary discs and help work out the story of planet formation!

Astrobite edited by Drew Lapeer.

Featured image credit: Andrews et al. 2018 

Further Resources:
“Intro to Astrophysical Simulation Jargon” Astrobite
“What is the meter size barrier?” Astrobite 
Dust growth and evolution in protoplanetary discs, Birnstiel 2024
DustPy documentation page
The Astrophysics of Planet Formation by Philip Armitage