-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnumerical-chemistry.tex
204 lines (190 loc) · 10.9 KB
/
numerical-chemistry.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
\subsection{Chemistry}
\label{sec.num.chemistry}
While it is often safe to assume that species (both chemical and
ionization) within a gas can be treated as being in equilibrium, in
some regimes that are found in astrophysics this assumption leads to
considerable error. For example, the cooling and collapse of
primordial gas in Population III star formation is dominated by
molecular hydrogen, which in the absence of dust forms via an
inefficient pair of collisional processes that depend heavily on the
local, highly non-equilibrium population of free electrons. As a
result, when modeling primordial star formation it is critical to
follow the non-equilibrium evolution of the chemical species of
hydrogen, including molecular hydrogen and deuterium.
The primordial non-equilibrium chemistry routines used in \enzo\ were
first described by \citet{abel97} and \citet{anninos97}, but have
since been extended with updated reaction rates and the inclusion of
deuterium species \citep{McGreer2008, 2009PhDT.........5T}. These
routines follow the non-equilibrium chemistry of a gas of primordial
composition with 12 total species: H, \Hp, He, \Hep, \Hepp, \Hm, \HHp,
\HH, e$^-$, D, \Dp, and HD. \enzo\ also computes the radiative
heating and cooling of the gas from atomic and molecular line
excitation, recombination, collisional excitation, free-free
transitions, Compton scattering of the cosmic microwave background, as
well as several models for a metagalactic UV background that heat the
gas via photoionization and photodissociation (see
Section~\ref{sec.num.cooling} for more details). The chemical and
thermal states of the gas can be updated either at the same
hydrodynamical timestep (i.e., decoupled and operator-split) or
through the same subcycling system (i.e., a coupled chemical and
thermal system). The default behavior of \enzo\ is to couple these
two systems at subcycles of the hydrodynamic timestep; this results in
updates to both the chemical and thermal states of the gas (which also
inform the temperature, the reaction rate coefficients and the cooling
functions of the gas) on timescales that are faster than those of the
gas dynamics. The gamma used by \enzo\ to compute the temperature of
the gas from the energy and density characteristics utilizes a
variable gamma that includes effects of the rotational state of
molecular hydrogen, enabling it to vary from $5/3$ (fully-atomic) to
$7/5$ (fully molecular). This further coupling of the chemical and
thermal states of the gas underscores the need for coupled chemistry
and radiative cooling solutions.
Input parameters to \enzo\ govern the chemical species that are
updated during the course of the simulation. This can include only
the atomic species (H, \Hp, He, \Hep, \Hepp, and e$^-$), those species
relevant for molecular hydrogen formation (\HH, \HHp, and \Hm), and
can further include deuterium and its species (D, \Dp, and HD). A
total of 9 kinetic equations are solved for the 12 species mentioned
above, including 29 kinetic and radiative processes. See
Table~\ref{table.collisional} for the collisional processes and
Table~\ref{table.radiative} for the radiative processes solved.
The chemical reaction equation network is technically challenging to
solve due to the huge range of reaction timescales involved. The
characteristic times for creation and destruction of the various
species and reactions can differ by many orders of magnitude and are
often very sensitive to the chemical and thermal state of the gas.
This makes a fully-implicit scheme, with convergence criteria and
error tolerance, strongly preferable for such a set of equations.
However, most implicit schemes require an iterative procedure to
converge, and for large networks (such as this one) an iterative,
fully-implicit method can be very time-consuming and computationally
costly for a relatively small increase in accuracy. At the present
time, this makes fully-implicit methods somewhat undesirable for
large, three-dimensional simulations.
\enzo\ solves the rate equations using a method based on a
semi-implicit formulation in order to provide a stable, positive
definite and first-order accurate solution. The update discretization
splits chemical changes into formation components and destruction
components and updates with a mixed set of time states, as described
in \citet{anninos97}. The formation components of species $S_i$ are
computed at the current subcycle time, where the contribution of
species $S_i$ to its own destruction components are computed at the
updated time; all other contributions to the destruction component are
computed at the current time. This mixed state improves accuracy,
ensures species values are positive definite, and is equivalent to one
Jacobi iteration of an implicit Euler solve. This technique is
optimized by taking the chemical intermediaries \Hm and \HHp, which
have large rate coefficients and low concentrations, and grouping them
into a separate category of equations. Due to their fast reactions,
these species are very sensitive to small changes in the more abundant
species and are (at almost all times in astrophysical calculations)
close to their equilibrium values. Attempting to resolve their
formation and destruction times would necessitate extremely small
timesteps. Therefore, reactions governing these two species can be
decoupled from the rest of the network and treated independently
through analytic solutions for equilibrium values. This allows a
significant speedup in the solution speed, as the timestepping scheme
is applied only to the slower 7- or 10-species network (depending on
whether deuterium is included or not), which will be much closer to
the overall hydrodynamic timestep of the simulation.
Even so, the accuracy and stability of the scheme is maintained by
subcycling the rate solver within a single hydrodynamic timestep.
These subcycle timesteps are determined so that the estimated
fractional change in the electron concentration is limited to no more
than $10\%$ per timestep; additional criteria may be applied based on
the expected change in internal energy from radiative cooling and from
chemical heating due to the formation of molecular hydrogen.
It is important to note the regime in which this model is valid.
According to \citet{abel97} and \citet{anninos97}, the reaction
network is valid for temperatures between $10^0 - 10^8$ K. The
original model discussed in these two references was only claimed to
be valid up to $n_H \sim 10^4$~cm$^{-3}$. However, addition of the
3-body \HH~formation process (equation 20 in
Table~\ref{table.collisional}) allows correct modeling of the
chemistry of the gas up until the point where collisionally-induced
emission from molecular hydrogen becomes an important cooling process,
which occurs at $n_{\rm H} \sim 10^{14}$~cm$^{-3}$. A further concern
is that the optically thin approximation for radiative cooling
eventually breaks down, which occurs before $n_{\rm H} \sim 10^{16} -
10^{17}$~cm$^{-3}$ in gas of primordial composition. Beyond this
point, modifications to the cooling function that take into account
the non-negligible opacity in the gas must be made, as discussed by
\citet{2004MNRAS.348.1019R}, and was put into \enzo\ for the work
published in \citep{2009Sci...325..601T,2009PhDT.........5T}. The
formation of molecular hydrogen as catalyzed by dust was recently
added to \enzo\ to enable studies of low-metallicity gas, as well as
the inclusion of appropriate timestepping criteria to account for the
input of ionizing radiation. Even with these modifications, a
completely correct description of the cooling of primordial gas at
very high densities requires some form of radiation transport, which
will greatly increase the cost of the simulations. Furthermore, at
very high densities, the stiffness of the molecular hydrogen reaction
rates may require better than a first-order accurate solution; as
such, the transition to this regime will likely necessitate a
fully-implicit, iterative solver.
%---------------- table of collisional processes
\begin{table}
\begin{center}
{\bfseries Collisional Processes}\\[1ex]
\begin{tabular}{llllllll}
(1) & H & + & e$^-$ & $\rightarrow$ & H$^+$ &+& 2e$^-$ \\
(2) & H$^+$ &+ &e$^-$ & $\rightarrow$ & H &+ &$\gamma$ \\
(3) & He &+& e$^-$ & $\rightarrow$ & He$^+$ &+& 2e$^-$ \\
(4) & He$^+$ &+& e$^-$ & $\rightarrow$ & He &+ &$\gamma$ \\
(5) & He$^{+}$ &+& e$^-$ & $\rightarrow$ & He$^{++}$ &+& 2$e^-$ \\
(6) & He$^{++}$ &+& e$^-$ & $\rightarrow$ & He$^+$ &+& $\gamma$ \\
\hline
(7) & H &+& e$^-$ &$\rightarrow$& H$^-$ &+& $\gamma$ \\
(8) & H$^-$ &+& H &$\rightarrow$ & H$_2$ & +& e$^-$ \\
(9) & H &+ &H$^+$ &$\rightarrow$ &H$_2^+$ &+ &$\gamma$ \\
(10) & H$_2^+$ &+ &H &$\rightarrow$ &$H_2$ &+ &$H^+$ \\
(11) & H$_2$ &+ &H$^+$ &$\rightarrow$ &H$_2^+$ & +& H \\
(12) & H$_2$ &+ &e$^-$ & $\rightarrow$ & 2H & + & e$^-$ \\
(13) & H$_2$ & + & H & $\rightarrow$ & 3H & & \\
(14) & H$^-$ & + & e$^-$ & $\rightarrow$ & H & + & 2e$^-$ \\
(15) & H$^-$ & + & H & $\rightarrow$ & 2H & + & e$^-$ \\
(16) & H$^-$ & + & H$^+$ & $\rightarrow$ & 2H & & \\
(17) & H$^-$ & + & H$^+$ & $\rightarrow$ & H$_2^+$ & + & e$^-$ \\
(18) & H$_2^+$ & + & e$^-$ & $\rightarrow$ & 2H & & \\
(19) & H$_2^+$ & + & H$^-$ & $\rightarrow$ & H$_2$ & + & H \\
(20) & 2H & + & H$_2$ & $\rightarrow$ & 2H$_2$ & & \\
(21) & 2H & + & H & $\rightarrow$ & H$_2$ & + & H \\
(22) & H$_2$ & + & H$_2$ & $\rightarrow$ & H$_2$ & + & 2H \\
(23) & 3H & & & $\rightarrow$ & H$_2$ & + & H \\
\hline
(24) & D & + & e$^-$ & $\rightarrow$ & D$^+$ &+& 2e$^-$ \\
(25) & D$^+$ &+ &e$^-$ & $\rightarrow$ & D &+ &$\gamma$ \\
(26) & H$^+$ &+ &D & $\rightarrow$ & H &+ &D$^+$ \\
(27) & H &+ &D$^+$ & $\rightarrow$ & H$^+$ &+ &D \\
(28) & H$_2$ &+ &D$^+$ & $\rightarrow$ & HD &+ &H$^+$ \\
(29) & HD &+ &H$^+$ & $\rightarrow$ & H$_2$ &+ &D$^+$ \\
(30) & H$_2$ &+ &D & $\rightarrow$ & HD &+ &H \\
(31) & HD &+ &H & $\rightarrow$ & H$_2$ &+ &D \\
(32) & H$^-$ &+ &D & $\rightarrow$ & HD &+ &e$^-$ \\
\end{tabular}
\caption[]{Collisional processes solved in the \enzo\ nonequilibrium
primordial chemistry routines.}
\label{table.collisional}
\end{center}
\end{table}
\begin{table}
\begin{center}
{\bfseries Radiative Processes}\\[1ex]
\begin{tabular}{llllllll}
(33) & H & + & $\gamma$ & $\rightarrow$ & H$^+$ & + & e$^-$ \\
(34) & He & + & $\gamma$ & $\rightarrow$ & He$^+$ & + & e$^-$ \\
(35) & He$^+$ & + & $\gamma$ & $\rightarrow$ & He$^{++}$ & + & e$^-$ \\
(36) & H$^-$ & + & $\gamma$ & $\rightarrow$ & H & + & e$^-$ \\
(37) & H$_2$ & + & $\gamma$ & $\rightarrow$ & H$_2^+$ & + & e$^-$ \\
(38) & H$_2^+$ & + & $\gamma$ & $\rightarrow$ & H & + & H$^+$ \\
(39) & H$_2^+$ & + & $\gamma$ & $\rightarrow$ & 2H$^+$ & + & e$^-$ \\
(40) & H$_2$ & + & $\gamma$ & $\rightarrow$ & H$_2^*$ & $\rightarrow$ & 2H \\
(41) & H$_2$ & + & $\gamma$ & $\rightarrow$ & 2H & & \\
(42) & D & + & $\gamma$ & $\rightarrow$ & D$^+$ & + & e$^-$ \\
\end{tabular}
\caption[]{Radiative processes solved in the \enzo\ nonequilibrium
primordial chemistry routines.}
%\red{Missing rates for HD destruction/ionization?}}
\label{table.radiative}
\end{center}
\end{table}