Time-dependent chemistry

Time-dependent model on a fixed structure

If a restart.fits.gz is available (i.e. from a steady-state model), one can simply use the time_dependent flag. In that case, the structure, temperature, and radiation field are kept fixed, and the disk chemistry is solved at every point in time-dependent mode from certain initial concentrations. These initial molecular concentrations are taken from a solution of the same reaction network for typical molecular cloud conditions, using a long integration time (For more details, see time_dependent_mc.f).

If you want to run such a time-dependent model, activate in Parameter.in ...

.true.   ! time_dependent           : solve time-dependent disk chemistry?
.true.   ! molecular_cloud          : initial concentrations from mol.cloud.?
6        ! N_age                    : no of successive output ages
1.0 1.E+3 1.E+4 1.E+5 1.E+6 1.E+7   : disk ages in years

It is recommended to use a converged model for this exercise and use the same parameter, grid resolution etc. All physical properties in the disk (gas and dust temperature, local radiation field, density structure) will stay the same. A number of output files for successive disk ages will be created: ProDiMo_0001.out, ProDiMo_0002.out, ... for the above-specified disk ages. The heating/cooling rates are recalculated with the resulting abundances, but the sum of heating and cooling rates will, in general, not be equal. The output files do not contain certain quantities, namely kappa_Ross, tau_Ross and the timescales tauchem, taucool, and taudiff.

The first step of time-dependent models is to evolve the chemistry from an initial "atomic" chemical state, under molecular cloud conditions, see Helling & Woitke (2014). In Parameter.in, you can set those molecular cloud conditions, and age of the cloud:

1.E+4       ! nH_mc
10.0        ! Tg_mc
10.0        ! Td_mc
10.0        ! Av_mc
0.1         ! a1_mc [micron]
1           ! N_age_mc
1.7E+5      ! [years]
1.36E-17    ! CRI_mc
.true.      ! Xray_workaround

These are the recommended values for TMC-1 in the UMIST2012 paper (McElroy+2012, A&A 550, A36) and also the default values in ProDiMo (so one does not have to set it). A very young age of TMC-1 is used to avoid too high O2O_2 concentrations, which would be in conflict with observations. Xray_workaround avoids some problems when doing the molecular cloud model with X-rays chemistry (not always required). To set the initial atomic state, modify Species.in as ...

100
H      1.E-20
H+     1.E-20
H-     1.E-20
H2     1.0
H2+    1.E-20
H3+    1.E-20
H2exc  1.E-20
He     1.0
He+    1.E-20
C      1.E-20
C+     1.0
C++    1.E-20
CH     1.E-20
CH+    1.E-20
...

it is recommended to set H2 and H as above (note that the positive notation means that the abundances are taken from Elements.in and a '1.0' indicates into which species that abundance is put initially), and put all other element abundances into He, C+, N, O, S+, Si+, Na+, Fe+, Ne, Ar.

Time-dependent chemistry with self-consistent gas temperature from thermal balance

The time-dependent mode time_chem_disk has more options but is also used for the time_dependent mode described above.

The advantage of the time_chem_disk mode is that it is more flexible. For example, one can also calculate the gas temperature or restart from another time-dependent model run. In the following, some typical use cases are described.

The following switch combination will ask ProDiMo to run first a default molecular cloud model, apply all the abundances to all the grid points, and use it as starting abundances for time-dependent chemistry for seven disk ages given in years (1, 1000, 1e4, 1e5, 5e5, 1e6, 2e6). At each time step, the chemistry is consistent with the thermal balance. The time for the model to reach thermal steady-state is low (a few 100 years) such that the age steps can safely increase. The outputs are (ProDiMo_01.out, FlineEstimates_01.out, restart_01.fits.gz, ...) where the 01 means that the outputs correspond to the first age of the disk (_02, _03, ...). The main switch is time_chem_disk.

.false.    ! restart
.true.     ! molecular_cloud
.true.     ! time_chem_disk
7          ! N_age
1.0 1e3 1e4 1e5 5e5 1e6 2e6

One can restart from any ProDiMo run:

.true.     ! restart
.false.    ! molecular_cloud
.true.     ! time_chem_disk
3          ! N_age
3e6 4e6 5e6

One needs to copy restart_0*.fits.gz to restart.fits.gz before launching ProDiMo and after renaming the other output files since the count will start again at 1. In the example, one starts with the abundances at 2e6 yrs and computes the chemistry and thermal balance for 3, 4, and 5 Myrs. The structure of the disk has not changed.

If the gas temperature is frozen

.true.  ! freeze_Tgas

then the chemistry will be run with the gas temperature from a previous run (restart.fits.gz).

The flag:

.true. ! Tgas_once

Tells the code to run the chemistry with the temperature in restart.fits.gz and recompute the thermal balance only once at the end of the chemistry.

One can, off-course, not freeze the gas temperature. In that case, the thermal balance will be computed at the end of each chemical step, and if the balance is not reached, the code will rerun the chemistry with a new guess temperature. On average, 6-10 guesses are needed to find the self-consistent Tgas and chemistry.

Setting the initial abundances

It can also be useful to set the initial species abundances for a time-dependent model by hand (in Species.in). This can be done like this:

.false.    ! molecular_cloud
.true.     ! time_chem_disk_ifspecies  : take initial abundances from Species.in

Please note in that case, it is recommended to use Species.in with the negative notation. Below is an example (e.g. CO has an absolute abundance of 10^-4.2868). This way, one can define the abundance of each species individually.

Please note: the user is responsible for making sure that the abundances given in Species.in are consistent with the element abundances in Elements.in (e.g. the sum of all carbon in Species.in should give the element abundance of carbon in Elements.in). There is no automatic way to do that at the moment. An example model using such a setup is in examples/StandardObjects/ClassI (within your ProDiMo source directory).

105
H            -4.3010
H+           -60.000000000
H-           -60.000000000
H2           -0.3011
H2+          -60.000000000
H3+          -60.000000000
H2exc        -60.000000000
He           -1.046
He+          -60.000000000
C            -60.000000000
C+           -60.000000000
CH           -60.000000000
CH+          -60.000000000
CH2          -60.000000000
CH2+         -60.000000000
CH3          -60.000000000
CH3+         -60.000000000
CH4          -60.000000000
CH4+         -60.000000000
CH5+         -60.000000000
CN           -60.000000000
CN+          -60.000000000
HCN          -60.000000000
HCN+         -60.000000000
HCNH+        -60.000000000
CO            -4.2868
CO+          -60.000000000
HCO          -60.000000000
HCO+         -60.000000000
...

Using DVODEF90 as an alternative chemical solver

Besides the standard LIMEX solver, one can also use the DVODEF90 solver. By setting

.true.    ! DVODEF90   : use DVODEF90 solver instead of LIMEX

DVODEF90 also works in parallel mode and can be faster in some cases (depending on the model). Tests have shown that both solvers give nearly identical results (also true for the gas temperature).