Chemistry Along Accretion Streams in Protoplanetary Disks
Ellen M. Price, Ilse Cleeves, and Karin Öberg
18 June 2020
Hi! My name is Ellen Price, and today I'll be giving an overview of some of the work I've done as part of my doctoral thesis at Harvard; this particular project was published earlier this year. I would like to begin by thanking my coauthors, Ilse Cleeves and Karin Öberg, for their guidance and mentorship, and the organizing committee for giving me the opportunity to speak today.
Disk anatomy and processes
(inspired by Henning & Semenov 2013)

Though I don't need to remind anyone at this conference, it still bears repeating that protoplanetary disks are very complex and dynamic systems. Disks have strong vertical and radial temperature gradients, which I've shown here on the right side of the figure. On the left, I show several dynamical processes that can take place. For example, the grain population can evolve with time by growing, fragmenting, and settling to the midplane. What I really want to draw your attention to, however, is the process of accretion, where disk material spirals in towards the star and, eventually, may fall onto the star. For the rest of this talk, I'll be assuming that accretion dominates and that other transport processes are negligible. Furthermore, I will restrict my focus to the disk midplane.
Why include accretion?
Accretion serves to change the local conditions of any given gas parcel
Planets form from the material available to them in the protoplanetary disk (gas and solids)
So a planet forming in an accreting disk is different from an analagous planet forming in a disk without accretion
It is worth thinking about the expected outcome of including accretion and its implications for planet formation. Intuitively, any given parcel of material, as it moves, is experiencing a range of physical conditions, such as changing temperature and changing local surface density, which influences the level of irradiation from cosmic rays. A planet forming from this material, then, will have a different composition than a planet forming at the same location in a disk without accretion.
Why this method?
Want something simple enough to be tractable, but complex enough to tell us something interesting
Other models may try to solve everything globally, which is unnecessary under our assumptions
The method I will present is local and fast!
You could also reasonably wonder why we wrote our own disk model instead of using one of the amazing existing ones. We wanted to develop a method that is simple enough to be computationally tractable, but complex enough to tell us something interesting. In contrast to a full, global solve of disk chemistry, the method I present here is largely local and relatively fast. In our method, surface density is solved first, then chemistry along accretion streams or tracks. What do I mean by local and global in reference to solvers? A global solve must take into account the behavior at all radii to find the next timestep's behavior. A local solve, on the other hand, need only consider the local properties of the disk to get the next timestep.
Methods: Surface density solve
$$ \frac{\partial \Sigma}{\partial t} - \frac{3}{R} \frac{\partial}{\partial R} \left[R^{1/2} \frac{\partial}{\partial R} \left(\nu \Sigma R^{1/2}\right)\right] = 0 $$

Nonlinear diffusion equation for $\Sigma = \int \rho~\mathrm{d} z$, from Lynden-Bell & Pringle
Need two boundary conditions!
Solved with Crank-Nicolson timestepping, finite difference derivatives
We start with the classical nonlinear diffusion equation for the surface density of an alpha-disk from Lynden-Bell and Pringle. This is the only part of the method that involves a global solve, and it is still very fast because we solve the equation in one dimension. I would note that boundary conditions become very important at this stage, as they would for any nonlinear diffusion equation. We use a Crank-Nicolson implicit timestepping scheme and simple finite difference spatial derivatives; that sounds very simple, but it is sufficient for what we are doing here. Before I move on from this slide, I want to point out the function nu, which encodes the disk viscosity, because...
Temperature is a circular problem!
Viscosity depends on temperature! We encounter a chicken-and-egg problem in trying to solve for the temperature because of the couplings between variables that I'm showing here. The surface density sets the optical depth, which determines the dust temperature. In the midplane, thermal coupling between the dust and gas is a reasonable assumption. Yet the gas temperature appears in the expression for viscosity (through the sound speed term), which again determines surface density. We need a way to find a self-consistent solution!
Methods: Temperature solve
$$T = T_0 \left(e^{-\psi \tau} + \omega\right) e^{\beta_0 \log x + \beta_1 \log^2 x}$$

Assume this flexible, parametric form for the temperature
Use RADMC-3d to generate “true” temperatures everywhere
Fit the function, update the surface density, and repeat until convergence
To get around this problem, we use an iterative procedure that includes a Monte Carlo radiative transport code, RADMC-3d. We assume the temperature function shown here and fit it over time and radius. You might wonder what theory leads to a temperature function like this one, but, truth be told, this is just a very flexible form that seems to work in practice. We also note here that we need to assume a parametric form, rather than simply using the RADMC-3d output, because we need derivatives of temperature in our surface density solver. After the temperature function is fitted, we update the parameters to the surface density solver and iterate again, eventually reaching convergence.
Methods: Solving for tracks
Once the surface density profile is set, we numerically solve for the accretion tracks through the disk by integrating the velocity field with various starting radii. Here, I'm showing a few such tracks and how different fields change as a parcel moves along the track. The 5 au track turns out to be most interesting, and you can see that it experiences a drastic rise in temperature but also a falling cosmic ray rate as it spirals in. I also show the local densities and surface densities for comparison.
Methods: Evolving chemistry
Use a C++ implementation of a modified Fogel chemistry model
Approximately 600 species and 6000 reactions
Finally, once we have computed the conditions along a track, we can simply run the chemistry forward under the appropriate physical conditions. We use a chemical reaction network that includes about 600 species and 6000 reactions. Most importantly, the network we use includes both gas-phase processes and surface chemistry. That surface chemistry is restricted to small grains, since we do not include growth or fragmentation of solids in this model.
Cautionary tale: Reaction order matters!
This is a purely numerical artifact that occurs because of finite floating point precision
When adding numbers of disparate orders of magnitude, the order in which they are added really matters
Experiment: What is the sum of $1$ and $10^{-9}$ added $10^9$ times in single precision?
I would like to mention at this point something important that I learned while working on this, for anyone working on similar models: The order in which you add together reaction rates can really matter here and can slightly shift your results in comparison to other models. You can understand this better by performing a simple experiment: Add ten to the minus nine to the number 1, ten to the nine times. You would obviously expect that the sum would be 2, but, if you do this calculation in single precision floating point, you will find that the answer is 1! Because rates can be very small or very large, the problem of adding them together in assembling the ODE is analagous. If anyone has questions about that, I'm happy to answer online or offline.
Methods: How should we measure change?
Every track has a starting radius, which we set, and an ending radius, which we clip to 1 au or 1 Myr
We compare the composition of a moving parcel to its corresponding static counterparts at the initial and final radii of the track
Imagine a parcel spiraling in towards the star. We can measure its composition at all times, but what do we compare those measurements to? We chose to evolve parcels of the same initial composition at the endpoints of the track, keeping everything else the same but not allowing the parcel to move. We call these models the static models. Then, at 1 Myr, we can compare the compositions of the moving, dynamic parcel with the two static ones. We decided that comparing the 5 au parcel to the static model at 1 au was the best comparison to make.
Results: Changing fields
Here, I am showing several panels that communicate how various fields — the density, surface density, and midplane temperature — change over time. You can see that all fields drop over time, which is expected as material accretes onto the star. You might wonder how we get a vertical structure since we only solved the surface density equation. We assume a vertical Gaussian gas distribution, which uniquely specifies the vertical structure of the disk.
Results: Accretion is important!
Here, I am showing one panel of a figure from our paper to highlight how accretion can really change the chemical outcome of a disk, particularly close to the star. Along the track beginning at 5 au, we see orders of magnitude enhancement and depletion of different species, particularly hydrocarbons, which have been highlighted with colored rings; the other colors indicate the heaviest atom in the species plotted. The most abundant species, like water and carbon monoxide, see virtually no change. The dashed, horizontal, gray line shows where a species experiences no change, so the ratio between the abundances in the dynamic and static model is unity. Above that line, the species is enhanced and, below that line, it is depleted. We show analagous figures in the paper for a few different cases.
Why does this happen?
Chemistry is (usually) fastest at high temperatures and high densities
Cosmic ray flux is highest at low surface densities, so CR-driven chemistry can happen far out in the disk and then the products travel inwards
Cosmic rays play a huge role in this phenomenon. As we all know, chemistry is typically fastest at high temperatures and high densities. Cosmic-ray driven chemistry, however, can happen rapidly in the outer disk because of the low surface density, despite its lower temperatures and lower densities than the inner disk. The material then migrates inward and is eventually found in the inner disk, where higher temperatures may encourage chemistry that re-processes it. In the end, much of these results can be attributed to cosmic rays, so measuring the cosmic ray rate better could have important implications for a disk model like this one!
Takeaways
Accretion changes the compositions along streams of material in the disk, potentially changing the compositions of planets that form there
Signs of accretion (like enhanced hydrocarbons) might be observable with JWST if vertical mixing is strong and lofts midplane material into the upper disk layers
Stay tuned for Part 2!
Finally, I'll go over a few conclusions from this talk. Accretion really does change the chemical outcome of the protoplanetary disk, and much of that change can be attributed to the changing cosmic ray rate along a track of material. Signs of accretion, like enhanced hydrocarbons in the inner disk, might be observable with a telescope like JWST if vertical mixing lofts midplane material into the upper layers of the disk. I am also working on a follow-up project for this one, so stay tuned for Part 2!
Questions?
Thanks for your attention and, with that, I'll take any questions you may have.