Candy

Chemistry ANd DYnamics in protoplanetary disks

• Candy

  • About

  • Installation and Requirements

  • Quickstart

    • Static runs

    • Diffusive runs

    • Growth runs

    • Continuing a run

  • Input parameters

    • Model

    • Phys

    • Abundances

  • Output files

    • outfile

    • Pebble composition

  • Astrochem Changes

    • Reaction types

  • chemdiff

  • Bugs

  • Changelog

About

CANDY (Chemistry ANd DYnamics) is a protoplanetary disk model for determining the time-dependent chemistry and a 1D vertical slice of a dynamic disk. CANDY includes chemistry, vertical diffusion, pebble growth, and ice sequestration. Details of the model are given in Van Clepper et al. 2022.

CANDY contains two submodules, a modified astrochem (Maret & Bergin 2015) directory and the chemdiff directory, where the python wrapper is contained. CANDY works by iteratively switching between chemistry and dynamics for a vertical Column at a given location in a protoplanetary disk. The column is split into adjacent cells, within which the chemical and physical properties are assumed to be constant. Each is allowed to evolve chemically independently for some period of time before gas and ice species diffuse vertically through the disk between adjacent cells. Meanwhile, small grains grow into larger pebbles, removing ice species from the active chemistry.

Installation and Requirements

Astrochem dependencies are given in the astrochem wiki, while chemdiff dependencies include openmpi and mpi4py. Recommended installation is using conda to create a virtual environment with all dependencies. A new environment can be created with all needed dependencies with

conda create --name newcandy -c conda-forge sundials=5.7.0 python cython numpy matplotlib h5py autoconf automake libtool mpi4py openmpi

Active the enviornment with

conda activate newcandy

or

source activate newcandy

To compile Astrochem run the following commands from the Astrochem subdirectory

./bootstrap
./configure CPPFLAGS="-I$CONDA_PREFIX/include" LDFLAGS="-Wl,-rpath,$CONDA_PREFIX/lib -L/$CONDA_PREFIX/lib"  --prefix=$CONDA_PREFIX
make
make install

Because candy.py is run from the parent directory, there is no need to adjust your python path!

Quickstart

To run candy, simply run

python candy.py examples/growth_example/growth_example.ini

from the command line, where growth_example.ini is the input file you wish to use.

Static runs

In a static candy run, there is no diffusion nor is there any grain growth.

To set up a static run of candy, set

difftime = 0

in the model parameters of the input file. Additionally, the chemtime and tf parameters must be set equal to one another. This is because in a static run, chemistry is called only once per cell. To avoid ambiguity in how long to run chemistry for, these parameters should be the same.

Diffusive runs

In a diffusive run, there is diffusion of material throughout the column, but grains do not grow into pebbles. This can be done by setting the difftime and chemtime model parameters, but the physical parameter growth_height should be set to

growth_height = 0

This sets the height below which pebbles will grow to 0, shutting off any pebble growth.

Note: for diffusive and growth runs, diffusion is the time limiting step. Therefore we must have

difftime <= chemtime <= tf

Growth runs

A growth run contains diffusion and pebble growth (and is the default setup for CANDY). The growth_height parameter sets the height (in disk scaleheight) below which pebbles will grow. Any [1] growth_height > 0 will result in pebble growth. The resulting pebble composition is saved in the pebfile (default: pebcomp.out) in the outputdir directory.

Pebble compositions are stored normalized to the column density of hydrogen, with units of cm\(^{-2}\).

Continuing a run

(>= v1.2.0) If a run of candy is stopped for any reason, the run can be continued using the -C or --Continue flag followed by the number of the output to continue from. For example, to continue a run from the 8th output, call

python candy.py -C 8 examples/growth_example/growth_example.ini

from the command line. This will start the run using abundances from outputs ending in _t8 at the time touts[8]

Input parameters

Input files for CANDY are split into three sections: model inputs [model] which controls solver parameters, physical inputs [phys] which sets the relevent physical values for chemistry, and the molecular abundances [abundances], which sets the intial abundances of each chemical species in the cells.

Any parameters that are not explicitly set in the input file will use the default values.

Note: Any relative paths to files are always from the parent directory (where candy.py is located). For example, if you had an input file in an input directory (say newcandy/inputs/candyinput.ini), the path to the chmfile will be relative to the newcandy directory. Output files will automatically be placed into the output directory outputdir.

The chmfile path should always be relative to the newcandy home directory. Optionally, you can set chmfile = None to setup a run with no active chemistry.

Model

• chmfile

  • Chemical network file to be used by astrochem. Default: astrochem/networks/umist12_x.chm

• pebfile

  • name of file where pebble compositions will be stored. Default: pebcomp.out

• ti

  • Start time of candy integration in years. Default: 0

• tf

  • End time of the candy integration in years. Default: 1.0e6

• ncells

  • Number of cells to split the column into. Default: 50

• chemtime

  • Maximum time in years to run chemistry in each cell before diffusion. Default: 500

• difftime

  • Maximum timestep for diffusion between cells and pebble growth. Default: 100

• touts

  • list of output times in years for the simulation (if diffusive or growth run). If not supplied, default output times are generated semi-logarithmically based on ti and tf (touts = [1e3, 2e3, 3e3,…, 1e4, 2e4, 3e4,…, 1e5, 1.5e5, 2e5,…, tf]) Note: touts must be strictly greater than 0 (see Known Bugs)

• outputdir

  • Output directory where output files will be stored (trailing / is not needed, this is always interpreted as a directory). Default: r00

• outfile

  • Name of npz file where all outputs will be saved. Default: output_abuns.npz

Phys

• r

  • radial location of integration. Default: 30

• r_units

  • units of radius, options are ‘au’ or ‘cm’. Default: au

• alpha

  • Alpha viscosity parameter of the disk. Default: 1.0e-3

• chi

  • external UV radiation field strength in Draine Units. Default: 50

• cosmic

  • Cosmic ionization rate in \(\mathrm{s}^{-1}\). Default: 1.3e-17

• grain_size

  • size of solid grains in microns. Default: 0.1

• dg0

  • Initial dust-to-gas mass ratio. Default: 0.01

• opacity

  • The UV opacity of the dust in \(\mathrm{cm}^2\mathrm{g}^{-1}\). Default: 1.0e5

• growth_timescale_factor

  • Sets the timescale of growth for small grains into large pebbles. The growth timescale is defined as \(\tau_\mathrm{grow} = a/\Omega\epsilon\), where a is the growth_timescale_factor. Default: 1

• growth_height

  • Scaleheights below which pebbles will form. Default: 1

• growth_delay_time

  • Time in years to wait before beginning pebble growth. Grains and ice will grow into pebbles when time >= growth_delay_time. Default: 0

Abundances

The initial abundances of each molecular species. Abundances are assumed to be constant throughout the column intially, default abundances are the same as listed in table 1 of Van Clepper et al. 2022, reproduced below:

Molecule	Abundance
\(\mathrm{H}_2\)	5.00(-1)
\(\mathrm{He}\)	9.75(-2)
\(\mathrm{NH}_3\)	1.45(-6)
\(\mathrm{H}_2\mathrm{O}\)	1.18(-4)
\(\mathrm{CO}\)	6.00(-5)
\(\mathrm{N}_2\)	2.00(-5)
\(\mathrm{CH}_4\)	2.00(-6)
\(\mathrm{CH}_3\mathrm{OH}\)	1.00(-6)
\(\mathrm{H}_2\mathrm{S}\)	1.91(-8)
\(\mathrm{CO}_2\)	5.00(-5)
\(\mathrm{HCN}\)	3.50(-7)
\(\mathrm{grain}\)	2.20(-12)

Output files

As candy runs, it will create subdirectories named z00/, z01/, etc… for each cell in the Column. Within each of these folders will be astrochem_output HDF5 files for each output time from the simulation. The default behavior is for these many files to be deleted at the end of the CANDY run, and all of the outputs are gathered into one .npz file. CANDY also creates two plain text files: the pebble composition (default pebcomp.out) and output times, times.out.

outfile

The main output file, outfile in the model parameters, is saved as a zipped, uncompressed .npz file. This file can be opened in Python, by importing the numpy package and using the np.load() function (see documentation here). For the default outfile name of ‘output_abuns’, the file could be loaded as follows:

import numpy as np
cdout = np.load('output_abuns.npz')

This returns an NpzFile object, containing the saved numpy arrays that are output from CANDY. The available array names can be found by accessing the files property

cdout.files
# out: ['times', 'species', 'abundances', 'heights', 'density', 'extinctions']

These arrays can be indexed similar to a python dictionary, using square brackets. For example, ‘times’ contains a list of output times.

times = cdout['times']
times
# out: array([  100.,   200.,   500.,  1000.,  2000.,  5000., 10000.])

The different arrays have different shapes dependant on which dimension they vary. For an output with nt times, nz cells, and nspec species, the outputs are as follows:

• times: 1D array, shape=(nt,)

  • Output times in years.

• species: 1D array, shape=(nspec,)

  • List of all species considered in the network.

• abundances: 3D array, shape=(nt, nz, nspec)

  • Abundances relative to hydrogen of all species, at all heights and all times.

• heights: 1D array, shape=(nz,)

  • The elevation above the midplane at the center of each cell in cm.

• density: 1D array, shape=(nz,)

  • The number density of hydrogen in a given cell in \(\mathrm{cm}^{-3}\).

• extinctions: 2D array, shape=(nt, nz)

  • The visual extinction, \(A_v\) of a given cell at a given time. Convert to the UV optical depth using the relation \(A_v = \tau_{UV}/3.02\)

See below for a (quick) example of how to read abundances, this is also demonstrated in an example in the plot_results.py plotting script.

import numpy as np
# load data
cdout = np.load('examples/growth_example/output_abuns.npz')
times = cdout['times']
specs = cdout['species']
abuns = cdout['abundances']
elevs = cdout['heights']
# Say we want to find the midplane abundance of H2O ice after 1000 years of evolution
# find the correct indices using np.where() function
t = np.where(times == 1000)[0][0]       # time = 1000 years
j = 0                                   # midplane, cell index = 0
k = np.where(specs == 'grainH2O')[0][0] # species index of water ice
abuns[t,j,k]
# out: 1.6693156020606896e-05

Pebble composition

The pebble composition file is saved as plain text, default name pebcomp.out. This file keeps track of the ice species that are sequestered away from the active chemistry on the mantle of growing pebbles. The abundances of each ice species are normalized to the column density of hydrogen, rather than the number density of hydrogen. This is because the entire column will only have one reservoir of pebbles. The routine for removing ice is in the grow_grains() function in the chemdiff/diffusion.py module, and this is written to the pebcomp.out file in the save_pebbles() function in chemdiff/candyio.py.

The pebble composition file has three columns, the first is the time of the output in years, second the name of the ice species, and the third the normalized abundance.

Astrochem Changes

The astrochem folder is copied from the astrochem github page, complete with documentation and installation instructions. Changes have been made within the src folder, and Astrochem can be installed as usual. If you already have astrochem installed on your machine, you can simply copy the src folder from here into your astrochem directory (replacing the default astrochem/src folder), then reinstall as usual.

All of the chemical rate reactions are calculated in the rate() function in astrochem/src/rates.c (line 226). The function takes in the alpha, beta, and gamma parameters in addition to the reaction type from the chemical network file for each reaction. The rates also depend on the physical parameters, including column densities, temperature, extinction, and ice abundances.

Reaction types

A full list of chemical reactions and their associated types considered in the expanded astrochem is included below

description	reacno
Electron-grain recombination	-1
H2 formation on grains	0
HD formation on grains	100
Cosmic-ray ionization	1
General, two body reactions	2-12
Photo-ionizationa and dissociation	13
Shielded dissociation of H2	14
Shielded dissociation of HD	114
Shielded photo-ionization of C	15
Shielded dissociation of C16O	16
Shielded dissociation of C17O	17
Shielded dissociation of C18O	18
Freezeout onto grains	20
Thermal desorption	21
Cosmic-ray desorption	22
Photo-desorption	23
Hydrogenation	25
Cosmic-ray induced reactions	41
Cosmic-ray desorption of CO	42
Cosmic-ray reactions involving He	43
secondary xray ionization of … H	60
… H2	61
… other molecules	62
Photoelectron production from PAH	70
Charge exchange with PAH	71
Alternate charge exchange with PAH	72
creation of excited H2	90
de-excitation of H2	91
reactions with excited H2	92

Some reaction rate calculations are also adjusted to be more similar to the DALI calculations (Burderer et al 2012). A full description of the handling of the chemical calculations is in the appendix of Van Clepper et al. 2022.

chemdiff

The python wrapper for for running astroCHEM with DIFFusion. This is entirely written in python, and successively calls astrochem in parallel at different vertical locations ranging from disk midplane to 5 scale heights. Input parameters are readin from custom .ini files – examples have been provided for growth, static, and diffusion setups.

The chemdiff directory contains all of the python modules. Of particular interest: - chemistry.py - Calls astrochem in parallel for each of the cells. - diffusion.py - Diffuses species between cells and contains grain growth functions. - disk.py - Contains the disk parameters, including surface density and temperature profiles. - utils.py - Contains miscellanious, possibly useful, functions for analysis.

Bugs

If you find a bug please let me know via email or provide a pull request to fix the issue! Below are some of the known bugs that I am planning to fix, but haven’t yet.

• KeyError: ‘grain’

  • This error can arise if the species ‘grain’ is not included somewhere in the network file. If ‘grain’ is not in the network file, then the abundance of that species is not returned, and candy loses track of the proper grain abundance. This can be fixed by appending a reaction where ‘grain’ is present but does nothing. For example:

  grain -> grain    0.00e+00  0.00e+00  0.00e+00  2  6830

• tout = 0

  • touts listed must be strictly greater than 0, otherwise astrochem will attempt to run for zero years resulting in an error.

• oserror: unable to synchronously open file

  • This error typically means a .h5 file is corrupted, generally this means that there was an issue with running astrochem, and so the astrochem_output.h5 file cannot be read. To troubleshoot this, I recommend cd into outdir/z00 and manually calling astrochem on the input.ini file that is generated by candy.

Changelog

Below is a list of versions and any changes made between successive versions. It is always recommended to use the most recent version when possible. - 1.2.0 + Include option to continue previous run using -C or --Continue flags + 1.2.1 * fixed bug when checking for if grain is defined in init_abuns + 1.2.2 + fixed bug with double counting ice between astrochem output and pebble abundances

• 1.1.0

  • Grain abun <-> dust-to-gas ratio functions now take grain_size in microns as second argument

  • Include warnings for initial grain abundance and compare with initial dust-to-gas ratio

  • Print the current version at start of the run

• 1.0.0

  • Initial version, major changes from CANDY

[1]

The growth_height must be larger than the z-value of the lowest cell in the column. The column is 5 scaleheights tall, split between ncells cells. Cell z-values are stored for the center of each cell, so growth_height must be above \(0.5*5/\) ncells.