Calculating quasinormal modes of Schwarzschild anti-de Sitter black holes
using the continued fraction method

Ramin G. Daghigh1, Michael D. Green2, and Jodin C. Morey3

1 Natural Sciences Department, Metropolitan State University, Saint Paul,
Minnesota, USA 55106

2 Mathematics and Statistics Department, Metropolitan State University, Saint
Paul, Minnesota, USA 55106

3 School of Mathematics, University of Minnesota, Minneapolis, Minnesota, USA
55455

###### Abstract

We investigate the scalar, gravitational, and electromagnetic quasinormal mode
spectra of Schwarzschild anti-de Sitter black holes using the numerical
continued fraction method. The spectra have similar, almost linear structures.
With a few exceptions, the low overtone quasinormal modes are consistent with
previously obtained results in the literature that use other numerical
techniques. The intermediate and high overtone quasinormal modes, in
comparison to the Schwarzschild case, converge very quickly to the asymptotic
formulas previously obtained by analytic monodromy techniques. We also
determine one of the purely imaginary modes in the gravitational spectrum to
be algebraically special. In addition, we find a connection between the
analytic asymptotic formulas and the purely imaginary modes. In particular,
these formulas can be used to predict the bifurcation of the lowest damped
electromagnetic modes. Finally, we find no high overtone quasinormal modes
with high oscillation frequency and low damping, which had been previously
predicted.

1\. Introduction

Black holes in anti-de Sitter (AdS) spacetime have attracted a great deal of
attention for multiple reasons. It was shown by Hawking and Page [1] that,
unlike black holes in flat spacetime, large black holes in AdS spacetime have
positive specific heat and can be in stable equilibrium with thermal radiation
at a fixed temperature. Also, according to the AdS/CFT correspondence proposed
by Maldacena [2], a large static black hole in AdS spacetime corresponds to a
thermal state of a system in conformal field theory (CFT). More specifically,
the dynamical time scale for the return to thermal equilibrium in CFT is equal
to the decay rate of the perturbation of a large black hole in AdS spacetime.

It is, therefore, important to determine the natural vibrational modes of
perturbations for black holes in AdS spacetime. These vibrational modes are
called quasinormal modes (QNMs), which are discrete and complex. The imaginary
part of the frequency signals the presence of damping, a necessary consequence
of boundary conditions that require energy to be carried away from the system.

The significance and stability of QNMs, when discontinuities are introduced to
the Regge-Wheeler or QNM potential, were investigated initially by Nollert in
[3]. Later, we refined and expanded upon Nollert’s results in [4]. This
investigation was followed by multiple papers [5, 6, 7, 8, 9, 10] on the
stability of the QNM spectrum. It has been shown that QNMs (including the
least damped mode) are unstable when discontinuities are introduced to the
Regge-Wheeler potential. However, the ringdown waveform stays stable in the
presence of small discontinuities. Therefore, one can still use the QNM
spectrum of the smooth potential to describe the behavior of the ringdown
waveform. Consequently, the stability issue is not a deterrent for us in this
paper. In addition, since our universe is not an AdS space, the importance of
AdS lies mainly in the theoretical and mathematical realm where we do not need
to introduce physical discontinuities (such as the presence of matter in the
vicinity of a black hole) to the QNM potential.

The low overtone QNMs of Schwarzschild-AdS black holes for scalar
perturbations were first calculated numerically by Horowitz and Hubeny [11].
They computed these modes in four, five, and seven spacetime dimensions, which
are of interest in the context of the AdS/CFT correspondence. This work was
completed by Zhu et al. [12] and Konoplya [13], who analyzed the scalar QNMs
for small Schwarzschild-AdS black holes in detail. Cardoso and Lemos [14]
extended these calculations to include electromagnetic and gravitational
perturbations. In addition, Moss and Norman [15] calculated roughly twenty QNM
roots for the gravitational perturbations of the Schwarzschild-AdS black hole
and showed the duality between axial and polar perturbations can be preserved
if we apply Dirichlet boundary conditions to axial perturbations, and Robin
boundary conditions to polar perturbations. Berti and Kokkotas confirmed all
the above results in [16] and extended their numerical calculations to
Reissner-Nordström black holes in AdS spacetime. Finally, the low overtone
QNMs of Dirac spinors are addressed by Giammatteo and Jing in [17]. Motivated
by AdS/CFT correspondence, many authors have extended the QNM calculations to
higher spacetime dimensions (some examples are given in [21]). All these
authors used either the method developed by Horowitz and Hubeny [11] or some
variation of the Fröbenius method. As far as we know, no one has applied
Leaver’s continued fraction method [18] to calculate the QNMs of
Schwarzschild-AdS black holes. The continued fraction method has proven to be
the most reliable numerical technique to calculate these modes. Also, the
continued fraction method with Nollert’s improvement [19] is able to extend
the numerical calculations to high overtone modes with very good precision.

The high overtone QNMs of Schwarzschild-AdS black holes were calculated
numerically for scalar, electromagnetic, and gravitational perturbations in
four spacetime dimensions by Cardoso, Konoplya, and Lemos in [20]. Analytic
calculations of the infinitely high overtone (asymptotic) QNMs of four
dimensional Schwarzschild-AdS black holes were done by Cardoso et al in [22]
using a method based on the monodromy technique developed by Motl and Neitzke
[23, 24]. Natario and Schiappa, in [25], generalized the analytic results to
include Schwarzschild-AdS black holes in dimensions greater than four and
Ghosh et al [26] generalized the analytic results further by studying the
infinitely high overtone QNMs of asymptotically non-flat black holes in a
generic way.

Two of us in [27] adapted the analytic technique based on the complex
coordinate WKB method developed by Andersson and Howls in [28] to calculate
the high overtone QNMs of Schwarzschild-AdS black holes in dimensions greater
than three. In addition to confirming the previously obtained results by
Natario and Schiappa, for large black holes in certain spacetime dimensions,
we showed that the analytic monodromy technique implies the existence of an
asymptotic region of the QNM spectrum where the real part of the QNM
frequencies approaches infinity while the damping rate approaches a finite
value. The interesting thing about these “highly real” modes is that they have
a damping rate less than the least damped QNM found by Horowitz and Hubeny in
[11]. Therefore, assuming they exist, these modes will be the dominant modes
in the context of the AdS/CFT correspondence. It was also shown by one of us
in [29] that the same highly real modes can also be obtained using the
monodromy technique developed by Natario and Schiappa in [25]. The authors of
[30] searched for these highly real modes using two different numerical
methods, power series and time evolution, but they were not able to find them.
Since the continued fraction method with Nollert’s improvement has proven to
be the most effective in searching for high overtone QNMs, we use it here to
search for these modes.

The paper is organized as follows. In section 2, we lay down the general
formalism. In section 3, we describe our numerical process. In section 4, we
present the results. In section 5 we compare the numerical results with
previously found analytic asymptotic formulas. Finally, we summarize our
findings in section 6.

2\. General Formalism

For a spherically symmetric black hole in four spacetime dimensions, all
perturbations are governed generically by a Schrödinger wave-like equation of
the form

${d^{2}\psi\over dr_{*}^{2}}+\left[\omega^{2}-V(r)\right]\psi=0~{},$ (2.1)

where $r$ is the radial coordinate, $V(r)$ is the QNM potential and $\omega$
turns out to be the QNM frequency. The Tortoise coordinate $r_{*}$ is defined
by

$dr_{*}={dr\over f(r)}~{},$ (2.2)

where $f(r)$ is the metric function of the form

$f(r)=1-{2M\over r}+\frac{r^{2}}{R^{2}}~{}.$ (2.3)

Here, $M$ is the ADM mass of the black hole and $R$ is the AdS radius.111The
value of the cosmological constant, $\Lambda$, is given by $\Lambda=-3/R^{2}.$
In this paper, we assume the perturbations depend on time as $e^{-i\omega t}$.
Consequently, in order to have damping, the imaginary part of $\omega$ must be
negative.

The effective potential in Eq. (2.1) is given in [14, 20] as

$V(r)=f(r)\left[\frac{\ell(\ell+1)}{r^{2}}-\frac{2M(s^{2}-1)}{r^{3}}+\frac{\xi}{R^{2}}\right],$
(2.4)

where $s$ is the spin of the perturbation, $\ell$ is the multipole number, and
the index $\xi$ is $2$ for scalar and $0$ for electromagnetic and
gravitational perturbations. Note that we only consider odd parity (axial)
gravitational perturbations, since there exists a symmetry that relates odd
and even parity perturbations as shown by Moss and Norman in [15].

The effective potential is zero at the event horizon ($r\rightarrow r_{+}$).
We choose our boundary conditions so that the asymptotic behavior of the
solutions is

$\psi(r)\approx\left\\{\begin{array}[]{ll}e^{-i\omega r_{*}}&\mbox{as
$r\rightarrow r_{+}$ ($r_{*}\rightarrow-\infty$) },\\\ 0&\mbox{as
$r\rightarrow\infty$ ($r_{*}\rightarrow 0$) ,}\end{array}\right.$ (2.5)

which represents an in-going wave at the event horizon and no waves at
infinity.

The wave equation (2.1) in the radial coordinate is

$f{d^{2}\psi\over dr^{2}}+f^{\prime}{d\psi\over
dr}-\left[\frac{\rho^{2}}{f}+\frac{\ell(\ell+1)}{r^{2}}-\frac{2M(s^{2}-1)}{r^{3}}+\frac{\xi}{R^{2}}\right]\psi=0~{},$
(2.6)

where $\rho=-i\omega$. It is easy to show that for scalar perturbations ($s=0$
and $\xi=1$) the asymptotic behaviors of the wavefunction at the boundaries
are

$\psi\overset{r\rightarrow r_{+}}{\,\,\,\,\hbox
to30.0pt{\rightarrowfill}\,\,\,\,}(r-r_{+})^{b\rho}~{}~{}\mbox{and}~{}~{}\psi\overset{r\rightarrow\infty}{\,\,\,\,\hbox
to30.0pt{\rightarrowfill}\,\,\,\,}r^{-2}~{},$ (2.7)

where

$b=\frac{R^{2}r_{+}}{3r_{+}^{2}+R^{2}}.$ (2.8)

We can now write the solution to the wave equation (2.6), with the desired
behavior at the boundaries, in the form

$\psi(r)=\left(\frac{r-r_{+}}{r}\right)^{b\rho}\left(\frac{r_{+}}{r}\right)^{2}\sum_{n=0}^{\infty}a_{n}\left(\frac{r-r_{+}}{r}\right)^{n}.$
(2.9)

For electromagnetic ($s=1$) and gravitational ($s=2$) perturbations, where
$\xi=0$, the asymptotic behaviors of the wavefunction at the boundaries are

$\psi\overset{r\rightarrow r_{+}}{\,\,\,\,\hbox
to30.0pt{\rightarrowfill}\,\,\,\,}(r-r_{+})^{b\rho}~{}~{}\mbox{and}~{}~{}\psi\overset{r\rightarrow\infty}{\,\,\,\,\hbox
to30.0pt{\rightarrowfill}\,\,\,\,}r^{-1}~{},$ (2.10)

where the constant $b$ is given in Eq. (2.8). We can now write the solution to
the wave equation (2.6), with the desired behavior at the boundaries, in the
form

$\psi(r)=\left(\frac{r-r_{+}}{r}\right)^{b\rho}\left(\frac{r_{+}}{r}\right)\sum_{n=0}^{\infty}a_{n}\left(\frac{r-r_{+}}{r}\right)^{n}.$
(2.11)

Substituting either (2.9) or (2.11) into the wave equation (2.6) leads to a
seven-term recurrence relation

$\alpha_{1}a_{1}+\beta_{1}a_{0}=0$ (2.12)
$\alpha_{2}a_{2}+\beta_{2}a_{1}+\gamma_{2}a_{0}=0$ (2.13)
$\alpha_{3}a_{3}+\beta_{3}a_{2}+\gamma_{3}a_{1}+\delta_{3}a_{0}=0$ (2.14)
$\alpha_{4}a_{4}+\beta_{4}a_{3}+\gamma_{4}a_{2}+\delta_{4}a_{1}+\zeta_{4}a_{0}=0$
(2.15)
$\alpha_{5}a_{5}+\beta_{5}a_{4}+\gamma_{5}a_{3}+\delta_{5}a_{2}+\zeta_{5}a_{1}+\eta_{5}a_{0}=0$
(2.16)
$\alpha_{n}a_{n}+\beta_{n}a_{n-1}+\gamma_{n}a_{n-2}+\delta_{n}a_{n-3}+\zeta_{n}a_{n-4}+\eta_{n}a_{n-5}+\kappa_{n}a_{n-6}=0$
(2.17)

where $n=6,7,8,\dots$ and $a_{0}$ is a constant which we can take to be 1. The
coefficients in the above recurrence relation for scalar perturbations are
given in Appendix A. The coefficients for gravitational and electromagnetic
perturbations are given in Appendix B.

The recurrence relation (2.17) can be reduced to a six-term recurrence
relation

$\alpha_{n}^{\prime}a_{n}+\beta_{n}^{\prime}a_{n-1}+\gamma_{n}^{\prime}a_{n-2}+\delta_{n}^{\prime}a_{n-3}+\zeta_{n}^{\prime}a_{n-4}+\eta_{n}^{\prime}a_{n-5}=0$
(2.18)

by eliminating the $\kappa_{n}$ term using Gaussian elimination for $n\geq 6$
as follows:

$\displaystyle\alpha_{n}^{\prime}$ $\displaystyle=$ $\displaystyle\alpha_{n}$
$\displaystyle\beta_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\beta_{n}-\frac{\beta_{n-1}\kappa_{n}}{\eta_{n-1}}$
$\displaystyle\gamma_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\gamma_{n}-\frac{\gamma_{n-1}\kappa_{n}}{\eta_{n-1}}$
$\displaystyle\delta_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\delta_{n}-\frac{\delta_{n-1}\kappa_{n}}{\eta_{n-1}}$
$\displaystyle\zeta_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\zeta_{n}-\frac{\zeta_{n-1}\kappa_{n}}{\eta_{n-1}}$
$\displaystyle\eta_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\eta_{n}-\frac{\eta_{n-1}\kappa_{n}}{\eta_{n-1}}.$

Repeating this procedure four times results in a three-term recurrence
relation

$\alpha_{n}^{(4)}a_{n}+\beta_{n}^{(4)}a_{n-1}+\gamma_{n}^{(4)}a_{n-2}=0.$
(2.19)

Following Leaver [18] we obtain a continued fraction:

$\displaystyle\beta_{1}^{(4)}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{1}^{(4)}\gamma_{2}^{(4)}}{\beta_{2}^{(4)}-\frac{\alpha_{2}^{(4)}\gamma_{3}^{(4)}}{\beta_{3}^{(4)}-\cdots}}$
(2.20) $\displaystyle=$
$\displaystyle\frac{\alpha_{1}^{(4)}\gamma_{2}^{(4)}}{\beta_{2}^{(4)}-}\frac{\alpha_{2}^{(4)}\gamma_{3}^{(4)}}{\beta_{3}^{(4)}-}\frac{\alpha_{3}^{(4)}\gamma_{4}^{(4)}}{\beta_{4}^{(4)}-}\cdots$

which is an equation in $\rho=-i\omega$. The solutions $\omega$ are the QNMs.
To evaluate this numerically we must truncate the continued fraction at some
depth. We use a technique developed by Nollert [19] to approximate the value
of the tail end of the continued fraction.

In principle, we should be able to find all the QNMs by finding the roots to
(2.20), but in practice our numerical methods only find a few roots using this
single equation. For this reason it is useful to “invert” the equation. The
$n$-th inversion, for $n=2,3,\dots$, of this continued fraction is defined to
be [18] :

$\beta_{n}^{(4)}-\frac{\alpha_{n-1}^{(4)}\gamma_{n}^{(4)}}{\beta_{n-1}^{(4)}-}\frac{\alpha_{n-2}^{(4)}\gamma_{n-1}^{(4)}}{\beta_{n-2}^{(4)}-}\cdots\frac{\alpha_{1}^{(4)}\gamma_{2}^{(4)}}{\beta_{1}^{(4)}}=\frac{\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}}{\beta_{n+1}^{(4)}-}\frac{\alpha_{n+1}^{(4)}\gamma_{n+2}^{(4)}}{\beta_{n+2}^{(4)}-}\frac{\alpha_{n+2}^{(4)}\gamma_{n+3}^{(4)}}{\beta_{n+3}^{(4)}-}\cdots$
(2.21)

We will call (2.20) the $1$st inversion. Different roots are easier to find
numerically in different inversions. Therefore, by changing inversions we are
able to find as many roots as we want. However, inverting (2.20) introduces
spurious roots (meaning they are not actual QNMs). What is common among these
spurious roots is that they satisfy the equation
$\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}=0$. One can show that if $\omega$ is a
root that appears or disappears between the $n$th and $(n+1)$th inversion then
it must be the case that $\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}=0$. Conversely,
if $\omega$ satisfies $\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}=0$ then it will not
be a root of any inversion past the $n$th (because the left-hand-side of Eq.
(2.21) which is zero, will appear in the denominator in any subsequent
inversion). Thus spurious roots can only appear for one inversion. It is
possible for QNMs to satisfy $\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}=0$, although
we know of only one example, which we discuss below.

Numerically, when solving for $\omega$, some roots will disappear between
inversions simply because they have a smaller basin of attraction in the new
inversion. As far as we can tell, however, QNMs numerically appear across
multiple inversions, whereas the roots of
$\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}=0$ only appear for one inversion. We do
not believe these roots are QNMs for the following reasons. They do not appear
in other numerical studies that use different methods, such as [11, 13, 14,
15, 16], or in analytic studies, such as [27, 33]. In fact, for each QNM found
in previous numerical studies we find an infinite family of spurious roots
branching from it. For the roots on these branches that we have checked, they
all satisfy the equation $\alpha_{n}^{(4)}\gamma_{n+1}^{(4)}=0$. For the above
reasons, we disregard any root that does not persist for at least two
inversions. The one exception is the first inversion, since the roots of Eq.
(2.20) are the QNMs. For example, the purely imaginary root $\omega=-2i$ for
gravitational waves with $\ell=2$, $r_{+}=1$ found in [14, 15, 20] (but not in
[16]) is only visible in the 1st inversion, where $\alpha_{1}^{(4)}=0$ and
$\beta_{1}^{(4)}=0$. This $\omega$ is not found in higher inversions which are
undefined when $\beta_{1}^{(4)}=0$.

3\. Numerical Process

For the purpose of brevity, in the rest of this paper, we choose units in
which $R=1$. That means $r$, $r_{*}$, $r_{+}$ and $t$ are expressed in units
of $R$. The QNM frequency $\omega$ is in units of $R^{-1}$ and the units for
the QNM potential $V(r)$ are $R^{-2}$. All of our numerical calculations
follow these units.

We can summarize our numerical process as follows. We begin by looking at
computer generated complex plots of the lowest inversions near the origin.
From this, we note the region containing the lowest roots. We then generate a
grid of points in this region, which we use as initial guesses for a root
finder to solve (2.21) over several inversions and two different depths.

Within each inversion, we keep the roots which exist in both depths (we
consider roots to be the same if the distance between them is less than
$10^{-8}$). From those roots, we then keep the ones which persist between
inversions. What remains are roots clustered around the QNMs. We throw out any
clusters whose radius is greater than $10^{-8}$. We then identify the
mathematical mean of each cluster, and keep it as a QNM.

With the first few QNMs in hand, we use the slope defined by the last two
(most damped) QNMs to predict the location of the next QNM. We then use a root
finder to locate that mode. We repeat the process until the desired number of
roots are found.

Finally, in the gravitational QNM spectrum, we observed some “special” roots
along the imaginary axis. These roots are numerically very difficult to find.
Therefore, we also generate thousands of points along the imaginary axis,
which we use as initial guesses for a root finder over several inversions and
depths. We repeat the process explained above to ensure any roots found are
QNMs. Purely imaginary QNMs appear to go to negative infinity as either $\ell$
increases or $r_{+}$ decreases. Our method gets slower as the absolute value
of $\omega$ increases because it requires a much higher depth in the continued
fraction and thus higher precision in each calculation. This makes them
increasingly difficult to find.

4\. Results

Using the above techniques we have generated a large number of QNMs given in
the tables below. Tables 1 and 2 show the QNMs for gravitational and scalar
perturbations, while tables 3 and 4 show the QNMs for electromagnetic
perturbations.

Table I. Gravitational QNMs for $\ell=2$ with various $r_{+}$ values
---
$n$ | $r_{+}=0.4$ | $r_{+}=1$ | $r_{+}=50$
| $-8.6206055i$ | $-2.0000000i$ | $-0.0266703i$
$1$ | $3.1620904-0.4309248i$ | $3.0331142-2.4042343i$ | $92.5047566-133.1901656i$
$2$ | $4.4090815-1.5134459i$ | $4.9607292-4.8981936i$ | $158.1091839-245.8221619i$
$3$ | $5.8818960-2.6637516i$ | $6.9053582-7.2897272i$ | $223.2778817-358.3814086i$
$4$ | $7.4543991-3.7971516i$ | $8.8546996-9.6604241i$ | $288.3378802-470.9112893i$
$5$ | $9.0724986-4.9039109i$ | $10.8078366-12.0234382i$ | $353.3567010-583.4291864i$
$6$ | $10.7109440-5.9893001i$ | $12.7638395-14.3826649i$ | $418.3561115-695.9413097i$
$7$ | $12.3584850-7.0599788i$ | $14.7219942-16.7396924i$ | $483.3450430-808.4502934i$
$8$ | $14.0102590-8.1209090i$ | $16.6817901-19.0953039i$ | $548.3277639-920.9574123i$
$9$ | $15.6641429-9.1753252i$ | $18.6428630-21.4499373i$ | $613.3065437-1033.4633473i$
$10$ | $17.3191697-10.2252823i$ | $20.6049491-23.8038599i$ | $678.2826879-1145.9684908i$
$25$ | $42.1793880-25.8188101i$ | $50.0886404-59.0817570i$ | $1652.8356481-2833.5188349i$
$50$ | $83.6542363-51.7010061i$ | $99.2863968-117.8490937i$ | $3277.0366624-5646.0863655i$
$100$ | $166.6415739-103.4216166i$ | $197.7192949-235.3674961i$ | $6525.4254523-11271.2171329i$
$250$ | $415.6592307-258.5370623i$ | $493.0642611-587.9038157i$ | $16270.5850765-28146.6068515i$
$500$ | $830.7198929-517.0407262i$ | $985.3294717-1175.4550896i$ | $32512.5165965-56272.2556955i$
$1000$ | $1660.8602827-1034.0357622i$ | $1969.8735020-2350.5524365i$ | $64996.3795527-112523.5531555i$

Table II. Scalar QNMs for $\ell=0$ with various $r_{+}$ values
---
$n$ | $r_{+}=0.4$ | $r_{+}=1$ | $r_{+}=50$
$1$ | $2.3629278-1.0064708i$ | $2.7982232-2.6712058i$ | $92.4936786-133.1932905i$
$2$ | $3.9785674-2.0727637i$ | $4.7584892-5.0375691i$ | $158.1007632-245.8243897i$
$3$ | $5.6170486-3.1195485i$ | $6.7192677-7.3944927i$ | $223.2708309-358.3832592i$
$4$ | $7.2637270-4.1610833i$ | $8.6822268-9.7485168i$ | $288.3316967-470.9129117i$
$5$ | $8.9144995-5.2002234i$ | $10.6466697-12.1012475i$ | $353.3511287-583.4306502i$
$6$ | $10.5676627-6.2380243i$ | $12.6121072-14.4532764i$ | $418.3509995-695.9426546i$
$7$ | $12.2223635-7.2749838i$ | $14.5782377-16.8048758i$ | $483.3402936-808.4515447i$
$8$ | $13.8781214-8.3113728i$ | $16.5448710-19.1561903i$ | $548.3233098-920.9585873i$
$9$ | $15.5346423-9.3473535i$ | $18.5118811-21.5073044i$ | $613.3023358-1033.4644586i$
$10$ | $17.1917350-10.3830299i$ | $20.4791821-23.8582713i$ | $678.2786896-1145.9695477i$
$25$ | $42.0743049-25.9043386i$ | $50.0028337-59.1163010i$ | $1652.8330984-2833.5195137i$
$50$ | $83.5707795-51.7596984i$ | $99.2235879-117.8736016i$ | $3277.0348548-5646.0868482i$
$100$ | $166.5781067-103.4627550i$ | $197.6738412-235.3848628i$ | $6525.4241725-11271.2174753i$
$250$ | $415.6166454-258.5630109i$ | $493.0349397-587.9148155i$ | $16270.5842665-28146.6070684i$
$500$ | $830.6889195-517.0590689i$ | $985.3085361-1175.4628718i$ | $32512.5160236-56272.2558490i$
$1000$ | $1660.8379527-1034.0487325i$ | $1969.8585979-2350.5579408i$ | $64996.3791475-112523.5532641i$

Table III. Electromagnetic QNMs for $\ell=1$ with various $r_{+}$ values
---
$n$ | $r_{+}=0.4$ | $r_{+}=1$ | $r_{+}=10$
$1$ | $2.3298090-0.4411910i$ | $2.1630231-1.6990926i$ | $-15.5057909i$
|  |  | $-28.4729072i$
$2$ | $3.6810663-1.4970895i$ | $3.8438195-4.1519361i$ | $8.2213479-48.1622512i$
$3$ | $5.1989859-2.5870273i$ | $5.6734729-6.5764558i$ | $20.6800750-71.4498997i$
$4$ | $6.7740513-3.6687077i$ | $7.5537242-8.9805378i$ | $33.2738238-94.3751851i$
$5$ | $8.3760254-4.7415670i$ | $9.4583847-11.3723844i$ | $45.9659842-117.1903869i$
$6$ | $9.9929806-5.8077365i$ | $11.3772229-13.7563338i$ | $58.7231435-139.9407128i$
$7$ | $11.6192611-6.8689450i$ | $13.3052644-16.1348222i$ | $71.5263567-162.6486670i$
$8$ | $13.2518153-7.9264047i$ | $15.2397363-18.5093289i$ | $84.3638491-185.3268708i$
$9$ | $14.8888361-8.9809611i$ | $17.1789429-20.8808126i$ | $97.2278283-207.9830821i$
$10$ | $16.5291767-10.0332155i$ | $19.1217742-23.2499273i$ | $110.1128787-230.6223968i$
$25$ | $41.2862642-25.6954654i$ | $48.4442954-58.6628581i$ | $304.5346705-569.3640074i$
$50$ | $82.7032719-51.6530851i$ | $97.5218549-117.5343929i$ | $630.0334920-1132.9423059i$
$100$ | $165.6386458-103.4553474i$ | $195.8333278-235.1581149i$ | $1282.2614028-2259.3279839i$
$250$ | $414.5891139-258.6834199i$ | $491.0151820-587.8349288i$ | $3241.0199331-5637.2325149i$
$500$ | $829.5981416-517.2744444i$ | $983.1552394-1175.4931925i$ | $6507.1087972-11266.1907295i$
$998$ | $1656.3653866-1032.2898448i$ | $1963.6351037-2345.9972772i$ | $13014.3923657-22478.3521464i$

Table IV. Electromagnetic QNMs for $\ell=1$ with various $r_{+}$ values
---
$n$ | $r_{+}=20$ | $r_{+}=50$
$1$ | $-30.2425917i$ | $-75.0959881i$
| $-59.1642407i$ | $-149.6543652i$
$2$ | $-98.5793625i$ | $-226.9561884i$
| $-102.4094897i$ | $-291.7433132i$
$3$ | $30.8947884-148.2650727i$ | $43.3602597-388.2720375i$
$4$ | $56.0241892-194.1418962i$ | $106.4001605-503.6209656i$
$5$ | $81.3308344-239.7811308i$ | $169.7166198-617.7889596i$
$6$ | $106.7526402-285.2820029i$ | $233.2722272-731.5912326i$
$7$ | $132.2585800-330.6927893i$ | $297.0113691-845.1528498i$
$8$ | $157.8283078-376.0405837i$ | $360.8940471-958.5463019i$
$9$ | $183.4479338-421.3419990i$ | $424.8912893-1071.8162100i$
$10$ | $209.1076026-466.6079166i$ | $488.9817705-1184.9918150i$
$25$ | $596.2180499-1143.8340830i$ | $1455.6342381-2877.9895678i$
$50$ | $1244.2718735-2270.5138181i$ | $3073.7307600-5694.3937547i$
$100$ | $2542.8189534-4522.3165071i$ | $6315.9664176-11323.2404591i$
$250$ | $6442.5908521-11275.2135685i$ | $16052.9560069-28203.4504456i$
$500$ | $12945.1855388-22528.2790603i$ | $32288.6938102-56332.7117793i$
$998$ | $25900.8192803-44942.9434135i$ | $64636.4399310-112362.5921925i$

In the tables, purely imaginary modes are numbered differently than the
others. In the gravitational case we do not assign these modes a number. In
the electromagnetic case, they appear in groups of two to which we assign the
same number. This numbering will be explained in the next section where we
discuss asymptotic behavior.

We also determined the QNM spectrum for higher values of $\ell$. Figures 1, 2,
and 3 show how these modes compare to the ones in the tables.

Figure 1: QNM spectrum for gravitational odd parity perturbations with
different values of $\ell$ and $r_{+}$.

Figure 2: QNM spectrum for electromagnetic perturbations with different values
of $\ell$ and $r_{+}$. For $r_{+}=50$ the roots for $\ell=1,2,3,4$ have
bifurcated into four roots along the imaginary axis consistent with [31, 32].

Figure 3: QNM spectrum for scalar perturbations with different values of
$\ell$ and $r_{+}$.

Where our results overlap with those previously published, we have agreement
except for the items listed below:

  * •

For electromagnetic modes with $\ell=2$ and $r_{+}=1$, the authors of [14]
find $2.496-1.579i$ as the first mode, where we find $3.223-1.384i$. Our
result agrees with [16].

  * •

For electromagnetic modes with $\ell=2$ and $r_{+}=50$ we find $-75.269i$
whereas [14] finds $-75.139i$. Our result agrees with [16].

  * •

For gravitational modes with $\ell=2$ and $r_{+}=1$, there is some confusion
over whether $\omega=-2i$ is a QNM. It first appears in [14] with a question
mark, it then appears in a graph in [15]. In [16], the authors say they find
no evidence for it. Finally it appears again in [20] (this time without the
question mark). We can verify that $-2i$ is indeed a QNM. In fact this mode is
“algebraically special,”meaning it is purely imaginary and can be identified
analytically. It is easy to verify that $\omega=-2i$ satisfies
$\alpha_{1}^{(4)}=\alpha_{1}=0$ and $\beta_{1}^{(4)}=\beta_{1}=0$, and thus
satisfies Eq. (2.20).

  * •

For gravitational modes with $\ell=2$, the authors of [16] found
$\omega=-1.335i/r_{+}$ is a good fit for the relation between $r_{+}$ and the
purely imaginary $\omega$ for $r_{+}\geq 10$. For $r_{+}\leq 1$ we obtain a
good fit with $\omega=2.4416i-4.4097i/r_{+}$. To get this fit we calculated
the purely imaginary $\omega$ for $r_{+}=0.3,0.4,...,1.0$.

  * •

We do not find the “highly real” QNMs that were calculated analytically in
[27, 29] for gravitational modes with $\ell=2$ and scalar modes with $\ell=0$
for large black holes. We searched in areas of the complex plane suggested by
the asymptotic formula in [27], but found no QNMs. However, there are issues
that make searching for these modes troublesome. Our search was limited to
$\omega$ with a relatively small real part ($\omega_{R}<65000,n\approx 500$)
because larger values of $\omega$ take too long to compute. The modes only are
expected to appear for large $r_{+}$ which gives rise to QNMs with large real
part. Our method gets slower as the absolute value of $\omega$ increases
because it requires a much higher depth in the continued fraction and thus
higher precision in each calculation. Considering the difficulty we had
finding purely imaginary modes with high absolute value, it is possible the
guesses need to be very close to the actual roots to find the highly real
modes. We searched for these modes using $r_{+}=50$. We also searched for them
by rewriting the QNM equation in the large $r_{+}$ limit (taking
$f\approx-{2M\over r}+\frac{r^{2}}{R^{2}}$). In both cases we were unable to
find the highly real QNMs.

5\. Comparison with the Asymptotic Spectrum

For scalar and gravitational perturbations, the asymptotic (high overtone)
values of $\omega$ have been calculated analytically in equation (3.39) of
[27] (and also equation (3.39) of [25]) as222We have changed the indexing in
[27], which is completely arbitrary, to match that in our tables.:

$\omega\eta=(n-1)\pi+\frac{5\pi}{4}+\frac{i}{2}\ln 2.$ (5.1)

for $n=1,2,\dots$ where $\eta$ is given by

$\eta=\sum_{j=0}^{2}\frac{1}{2k_{j}}\ln\left(-\frac{1}{R_{j}}\right).$ (5.2)

$R_{j}$ are the zeros of the metric function $f(r)$ given in Eq. (2.3) and
$k_{j}=f^{\prime}(R_{j})/2$.

The value of $\eta$ and the explicit asymptotic formulas for $\omega$ for
various $r_{+}$ values are given in Table V. Although these formulas were
analytically calculated as asymptotic approximations, the actual scalar and
gravitational QNMs agree with the formulas remarkably well even for small
values of $n$. This agreement can be seen in Figure 4. There, we only show the
gravitational modes for the $r_{+}=0.4$ case, but the agreement is similar for
all $r_{+}$ values for both scalar and gravitational modes.

Figure 4: Comparison of the asymptotic formula (5.1) with the actual QNMs. The
agreement is even better for larger $r_{+}$ values.

The asymptotic values of $\omega$ for electromagnetic perturbations in the
$4$-dimensional Schwarzschild-AdS spacetime can be derived from Eq. (2.90) of
[33] by replacing $j=1$ and $\beta=\pi/2$ (see Table 6 in [33]). This leads to
the expression:

$\omega\eta=(n-1)\pi+\frac{5\pi}{4}-\frac{i}{2}\ln 3.$ (5.3)

This gives the correct asymptotic behavior. However, unlike the scalar and
gravitational cases, this formula does not give good results for smaller $n$
values. The approximation becomes worse for larger values of $r_{+}$. Based on
our calculated values, we can modify Eq. (5.3) to

$\omega\eta=(n-1)\pi+\frac{4\pi}{5}-i\ln[3.75(r_{+}+0.3)],$ (5.4)

which gives better results for a wide range of $r_{+}$. We show the difference
in these two asymptotic formulas in the figures below.

Figure 5: Comparison of the asymptotic approximation (5.3) and the improved
asymptotic approximation (5.4) with the actual QNMs with $r_{+}=50$. Figure 6:
Comparison of the asymptotic approximation (5.3) and the improved asymptotic
approximation (5.4) with the actual QNMs with $r_{+}=0.4$.

To find the improved asymptotic formula (5.4) we computed
$\omega_{n}\eta-(n-1)\pi$ for $r_{+}=0.4,1,10,20,50$ and the largest value of
$n$ we have in Tables III & IV ($n=998$). We found the real part was
approximately equal to $4\pi/5$ and we fit $\omega_{I}(r_{+})$ with a function
of the form $\ln[a(r_{+}+b)]$.

Table V. Values of $\eta$ and the corresponding asymptotic formulas. |
---|---
Grav/Scalar, EM1, EM2 refer to Eqs. (5.1), (5.3), (5.4), respectively. |
$r_{+}$ | $\eta$ |  Asymptotic Formula for $\omega$ | Type
$0.4$ | $1.3634+0.849079i$ | $(2.18944-1.10931i)+(1.6603-1.03398i)(n-1)$ | Grav/Scalar
|  | $(1.89459-1.58277i)+(1.6603-1.03398i)(n-1)$ | EM1
|  | $(1.01061-1.33722i)+(1.6603-1.03398i)(n-1)$ | EM2
$1$ | $0.658045+0.785398i$ | $(2.72065-2.72051i)+(1.9691-2.35019i)(n-1)$ | Grav/Scalar
|  | $(2.05045-3.28203i)+(1.9691-2.35019i)(n-1)$ | EM1
|  | $(0.390221-2.87305i)+(1.9691-2.35019i)(n-1)$ | EM2
$10$ | $0.0605905+0.104372i$ | $(18.8202-26.6994i)+(13.0693-22.5129i)(n-1)$ | Grav/Scalar
|  | $(12.4003-30.4263i)+(13.0693-22.5129i)(n-1)$ | EM1
|  | $(-15.7287-33.2109i)+(13.0693-22.5129i)(n-1)$ | EM2
$20$ | $0.0302464+0.0523163i$ | $(37.4904-53.3876i)+(26.0203-45.0065i)(n-1)$ | Grav/Scalar
|  | $(40.3947-51.7085i)+(26.0203-45.0065i)(n-1)$ | EM1
|  | $(-41.2494-71.8882i)+(26.0203-45.0065i)(n-1)$ | EM2
$50$ | $0.012093+0.0209412i$ | $(93.6207-133.461i)+(64.9677-112.503i)(n-1)$ | Grav/Scalar
|  | $(61.5386-151.988i)+(64.9677-112.503i)(n-1)$ | EM1
|  | $(-135.665-198.36i)+(64.9677-112.503i)(n-1)$ | EM2

Figures 4, 5, 6 show that formulas (5.1) and (5.4) are good approximations to
the QNM spectrum for all perturbations (gravitational, scalar,
electromagnetic). The formulas work best for the lowest multipole number
$\ell$. However, for all values of $\ell$, the QNMs are well approximated by
the values in the asymptotic formula for large enough $n$.

In the case of gravitational modes (Figure 4), each QNM corresponds to exactly
one mode given by Eq. (5.1) except for the purely imaginary one. Therefore,
the purely imaginary modes are “special” in the sense that the asymptotic
formula does not assign them a number. In Figure 5, we see that each value of
$n$ in the formula for EM2 from Table V again corresponds to exactly one
$\omega_{n}$ with the exception of the purely imaginary modes. Each pair of
purely imaginary modes corresponds to exactly one mode given by the asymptotic
formula. It turns out that as $r_{+}$ increases, the modes move toward the
imaginary axis. When a mode hits the imaginary axis it splits into two purely
imaginary modes. This is the bifurcation discussed in [31, 32]. Because of
this bifurcation, in Tables III & IV, the purely imaginary modes are grouped
in pairs, corresponding to a single value of $n$. We find it intriguing that
one could have predicted the bifurcation by matching the QNMs to the
asymptotic predictions.

6\. Summary

We conducted a comprehensive investigation of the gravitational, scalar, and
electromagnetic QNM spectra of a Schwarzschild black hole in AdS spacetime
using the numerical continued fraction method. With a few noted exceptions,
the low overtone QNMs that we found are consistent with previously obtained
results in the literature that use other numerical techniques. The
intermediate and high overtone QNMs converge quickly to the asymptotic
formulas previously obtained by analytic monodromy techniques. Some of the
highlights of our results are

  * •

For gravitational modes with $\ell=2$ and $r_{+}=1$, we verified analytically
that the purely imaginary root, $-2i$, is an algebraically special frequency.

  * •

Using our numerical results, we were able to improve the asymptotic formula
obtained in [33] for high overtone electromagnetic QNMs. Unlike the previous
formula, our asymptotic formula works well for a wide range of QNMs and black
hole sizes.

  * •

The asymptotic formula obtained using analytic monodromy techniques, in
addition to matching well to the high overtone QNMs, is able to predict the
bifurcation phenomenon at the low overtone region of the electromagnetic QNM
spectrum.

  * •

In our search, we were unable to find the “highly real” QNMs that were
calculated analytically in [27, 29] for gravitational modes with $\ell=2$ and
scalar modes with $\ell=0$ for large black holes. This makes it less likely
that these modes exist. We discussed issues that make searching for these
modes troublesome. We consequently have no conclusive remark regarding the
existence/nonexistence of these modes.

  * •

In the gravitational QNM spectrum, we observed some purely imaginary modes
that appear to go to negative infinity as either $\ell$ increases or $r_{+}$
decreases. We are surprised these modes have not appeared in the literature.
This is an interesting part of the QNM spectrum that can potentially be
explored using the analytic monodromy techniques developed to calculate highly
damped modes.

Acknowledgments

We are grateful to Kostas Kokkotas and Emanuele Berti for sharing with us
their insight on the numerical method.

Appendix A Scalar Coefficients

For scalar perturbations, the coefficients for the recursion relations (2.12
\- 2.17) are

$\displaystyle\alpha_{n}$ $\displaystyle=$
$\displaystyle-(1+3r_{+}^{2})^{2}(n+b\rho)^{2}+r_{+}^{2}\rho^{2}.$ (6.1)
$\displaystyle\beta_{n}$ $\displaystyle=$ $\displaystyle
3+\ell(\ell+1)(1+3r_{+}^{2})-6b\rho+6b^{2}\rho^{2}$ (6.2)
$\displaystyle+15r_{+}^{2}+18r_{+}^{4}+6n^{2}(1+5r_{+}^{2}+6r_{+}^{4})-27r_{+}^{2}b\rho-27r_{+}^{4}b\rho-2r_{+}^{2}\rho^{2}$
$\displaystyle+30r_{+}^{2}b^{2}\rho^{2}+36r_{+}^{4}b^{2}\rho^{2}-3n(1+3r_{+}^{2})(2+3r_{+}^{2}-4b\rho-8r_{+}^{2}b\rho).$
$\displaystyle\gamma_{n}$ $\displaystyle=$ $\displaystyle
2r_{+}^{2}(1+n+b\rho)-6+5\ell(\ell+1)+17n-13n^{2}+17b\rho-26nb\rho-13b^{2}\rho^{2}+r_{+}^{2}\rho^{2}$
(6.3)
$\displaystyle-(1+r_{+}^{2})^{2}\left[64+60n^{2}-107b\rho+60b^{2}\rho^{2}-107n+120nb\rho\right]$
$\displaystyle-2r_{+}^{2}(1+r_{+}^{2})(1+4n+4b\rho)$
$\displaystyle-(1+r_{+}^{2})\left[-49+9\ell(\ell+1)+94n-58n^{2}+94b\rho-116nb\rho-58b^{2}\rho^{2}\right].$
$\displaystyle\delta_{n}$ $\displaystyle=$
$\displaystyle-4\ell(\ell+1)+6(-1+n+b\rho)^{2}+2r_{+}^{2}(1+r_{+}^{2})(n+b\rho)$
$\displaystyle+(1+r_{+}^{2})^{2}\left[117+54n^{2}-146b\rho+54b^{2}\rho^{2}+2n(-73+54b\rho)\right]$
$\displaystyle+(1+r_{+}^{2})\left[-69+10\ell(\ell+1)+98n-40n^{2}+98b\rho-80nb\rho-40b^{2}\rho^{2}\right].$
$\displaystyle\zeta_{n}$ $\displaystyle=$
$\displaystyle-2+\ell(\ell+1)+3n-n^{2}+3b\rho-2nb\rho-b^{2}\rho^{2}$
$\displaystyle+(1+r_{+}^{2})\left[-43+5\ell(\ell+1)+48n-14n^{2}+48b\rho-28nb\rho-14b^{2}\rho^{2}\right]$
$\displaystyle-(1+r_{+}^{2})^{2}\left[107+28n^{2}-105b\rho+28b^{2}\rho^{2}+7n(8b\rho-15)\right].$
$\displaystyle\eta_{n}$ $\displaystyle=$
$\displaystyle(1+r_{+}^{2})\left[-10+\ell(\ell+1)+9n-2n^{2}+9b\rho-4nb\rho-2b^{2}\rho^{2}\right]$
(6.6)
$\displaystyle+(1+r_{+}^{2})^{2}\left[49+8n^{2}-39b\rho+8b^{2}\rho^{2}+n(16b\rho-39)\right].$
$\displaystyle\kappa_{n}$ $\displaystyle=$
$\displaystyle-(1+r_{+}^{2})^{2}(-3+n+b\rho)^{2}.$ (6.7)

Appendix B Gravitational and Electromagnetic Coefficients

For electromagnetic and gravitational perturbations, the coefficients for the
recursion relations (2.12 \- 2.17) are

$\displaystyle\alpha_{n}$ $\displaystyle=$
$\displaystyle-(1+3r_{+}^{2})^{2}(n+b\rho)^{2}+r_{+}^{2}\rho^{2}$ (6.8)
$\displaystyle\beta_{n}$ $\displaystyle=$ $\displaystyle
4+\ell(\ell+1)(1+3r_{+}^{2})-8b\rho+6b^{2}\rho^{2}-s^{2}(1+4r_{+}^{2}+3r_{+}^{4})$
(6.9)
$\displaystyle+19r_{+}^{2}+21r_{+}^{4}+6n^{2}(1+5r_{+}^{2}+6r_{+}^{4})-39r_{+}^{2}b\rho-45r_{+}^{4}b\rho-2r_{+}^{2}\rho^{2}$
$\displaystyle+30r_{+}^{2}b^{2}\rho^{2}+36r_{+}^{4}b^{2}\rho^{2}-n(1+3r_{+}^{2})(8+15r_{+}^{2}-12b\rho-24r_{+}^{2}b\rho)$
$\displaystyle\gamma_{n}$ $\displaystyle=$
$\displaystyle-22+5\ell(\ell+1)-13n^{2}-2r_{+}^{2}+33b\rho+2r_{+}^{2}b\rho$
(6.10)
$\displaystyle~{}~{}~{}~{}~{}-13b^{2}\rho^{2}+r_{+}^{2}\rho^{2}+n(33+2r_{+}^{2}-26b\rho)$
$\displaystyle-(1+r_{+}^{2})^{2}\left[124+60n^{2}-12s^{2}-161b\rho+60b^{2}\rho^{2}-161n+120nb\rho\right]$
$\displaystyle-(1+r_{+}^{2})\Big{[}-114+9\ell(\ell+1)-58n^{2}-10r_{+}^{2}+7s^{2}+154b\rho+8r_{+}^{2}b\rho$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-58b^{2}\rho^{2}+2n(77+4r_{+}^{2}-58b\rho)\Big{]}.$
$\displaystyle\delta_{n}$ $\displaystyle=$ $\displaystyle
20-4\ell(\ell+1)+6n^{2}-22b\rho+6b^{2}\rho^{2}-22n+12nb\rho$ (6.11)
$\displaystyle+(1+r_{+}^{2})^{2}\left[224+54n^{2}-19s^{2}-212b\rho+54b^{2}\rho^{2}+4n(-53+27b\rho)\right]$
$\displaystyle+(1+r_{+}^{2})\Big{[}-156+10\ell(\ell+1)-40n^{2}-4r_{+}^{2}+9s^{2}+154b\rho+2r_{+}^{2}b\rho$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-40b^{2}\rho^{2}+2n(77+r_{+}^{2}-40b\rho)\Big{]}.$
$\displaystyle\zeta_{n}$ $\displaystyle=$
$\displaystyle-6+\ell(\ell+1)+5n-n^{2}+5b\rho-2nb\rho-b^{2}\rho^{2}$
$\displaystyle-(1+r_{+}^{2})\left[-94+5\ell(\ell+1)+72n-14n^{2}+5s^{2}+72b\rho-28nb\rho-14b^{2}\rho^{2}\right]$
$\displaystyle-(1+r_{+}^{2})^{2}\left[200+28n^{2}-15s^{2}-147b\rho+28b^{2}\rho^{2}+7n(8b\rho-21)\right].$
$\displaystyle\eta_{n}$ $\displaystyle=$
$\displaystyle(1+r_{+}^{2})\left[-21+\ell(\ell+1)+13n-2n^{2}+s^{2}+13b\rho-4nb\rho-2b^{2}\rho^{2}\right]$
(6.13)
$\displaystyle+(1+r_{+}^{2})^{2}\left[89+8n^{2}-6s^{2}-53b\rho+8b^{2}\rho^{2}+n(16b\rho-53)\right].$
$\displaystyle\kappa_{n}$ $\displaystyle=$
$\displaystyle-(1+r_{+}^{2})^{2}[n^{2}-s^{2}-2n(4-b\rho)+(4-b\rho)^{2}].$
(6.14)

References

## References

  * [1] S.W. Hawking and D.N. Page, Thermodynamics of black holes in anti-de Sitter space, Comm. Math. Phys. 87 (1983) 588.
  * [2] J.M. Maldacena, The large N limit of superconformal field theories and supergravity, Adv. Theoret. Math. Phys. 2 (1998) 231; The Large-N Limit of Superconformal Field Theories and Supergravity, Int. J. Theor. Phys. 38 (1999) 1113.
  * [3] H.-P. Nollert, About the Significance of QNMs of Black Holes, Phys. Rev. D53 (1996) 4397.
  * [4] R.G. Daghigh, M.D. Green, J.C. Morey, Significance of Black Hole Quasinormal Modes: A Closer Look, Phys. Rev. D101 (2020) 104009.
  * [5] J.L. Jaramillo, R.P. Macedo, and L. Al Sheikh, Pseudospectrum and Black Hole Quasinormal Mode Instability, Phys. Rev. X 11 (2021) 031003.
  * [6] W.-L. Qian, K. Lin, C.-Y. Shao, B. Wang, R.-H. Yue, Asymptotical quasinormal mode spectrum for piecewise approximate effective potential, Phys. Rev. D103 (2021) 024019.
  * [7] K. Destounis, R.P. Macedo, E. Berti, V. Cardoso, J.L. Jaramillo, Pseudospectrum of Reissner-Nordström black holes: Quasinormal mode instability and universality, Phys. Rev. D104 (2021) 084091.
  * [8] J.L. Jaramillo, R.P. Macedo, L. Al Sheikh, Gravitational Wave Signatures of Black Hole Quasinormal Mode Instability, Phys. Rev. Lett. 128 (2022) 211102.
  * [9] E. Gasperin, J.L. Jaramillo, Energy scales and black hole pseudospectra: the structural role of the scalar product, Class. Quant. Grav. 39 (2022) 115010.
  * [10] M. Ho-Yeuk Cheung, K. Destounis, R.P. Macedo, E. Berti, V. Cardoso, Destabilizing the Fundamental Mode of Black Holes: The Elephant and the Flea, Phys. Rev. Lett. 128 (2022) 111103.
  * [11] G.T. Horowitz, V.E. Hubeny, Quasinormal modes of AdS black holes and the approach to thermal equilibrium, Phys. Rev. D62 (2000) 024027.
  * [12] J.-M. Zhu, B. Wang and E. Abdalla, Object picture of quasinormal ringing on the background of small Schwarzschild anti–de Sitter black holes, Phys. Rev. D63 (2001) 124004;
  * [13] R.A. Konoplya, Quasinormal modes of a small Schwarzschild-anti-de Sitter black hole, Phys. Rev. D66 (2002) 044009.
  * [14] V. Cardoso and J.P.S. Lemos, Quasinormal modes of Schwarzschild-anti-de Sitter black holes: Electromagnetic and gravitational perturbations, Phys. Rev. D64 (2001) 084017.
  * [15] I.G. Moss and J.P. Norman, Classical and Quantum Gravity Gravitational quasinormal modes for anti-de Sitter black holes, Class. Quant. Grav. 19 (2002) 2323.
  * [16] E. Berti and K.D. Kokkotas, Quasinormal modes of Reissner-Nordström-anti-de Sitter black holes: Scalar, electromagnetic, and gravitational perturbations, Phys. Rev. D67 (2003) 064020.
  * [17] M. Giammatteo and J. Jing, Dirac quasinormal frequencies in Schwarzschild-AdS space-time, Phys. Rev. D71 (2005) 024007.
  * [18] E.W. Leaver, An analytic representation for the quasi-normal modes of Kerr black holes, Proc. Roy. Soc. Lond. A 402 (1985) 285.
  * [19] H.-P. Nollert, Quasinormal modes of Schwarzschild black holes: The determination of quasinormal frequencies with very large imaginary parts, Phys. Rev. D47 (1993) 5253.
  * [20] V. Cardoso, R. Konoplya, and J.P.S. Lemos, Quasinormal frequencies of Schwarzschild black holes in anti-de Sitter spacetimes: A complete study of the overtone asymptotic behavior, Phys. Rev. D68 (2003) 044024.
  * [21] T. R. Govindarajan and V. Suneeta, Quasi-normal modes of AdS black holes: a superpotential approach, Class. Quant. Grav. 18 (2001) 265; A.O. Starinets, Quasinormal modes of near extremal black branes, Phys. Rev. D66 (2002) 124013; S. Musiri and G. Siopsis, Quasinormal modes of large AdS black holes, Phys. Lett. B563 (2003) 102; S. Musiri and G. Siopsis, Asymptotic form of quasi-normal modes of large AdS black holes, Phys. Lett. B576 (2003) 309; A. Nunez and A.O. Starinets, AdS/CFT correspondence, quasinormal modes, and thermal correlators in $\mathcal{N}=4$ supersymmetric Yang-Mills theory, Phys. Rev. D67 (2003) 124013; G. Siopsis, Large mass expansion of quasi-normal modes in AdS5, Phys. Lett. B590 (2004) 105; G. Siopsis, On quasi-normal modes and the AdS5/CFT4 correspondence, Nucl. Phys. B715 (2005) 483; S. Musiri, S. Ness and G. Siopsis, Perturbative calculation of quasinormal modes of AdS Schwarzschild black holes, Phys. Rev. D73 (2006) 064001; A. Lopez-Ortega, Electromagnetic quasinormal modes of D-dimensional black holes, Gen. Rel. Grav. 38 (2006) 1747; J.J. Friess, S. S. Gubser, G. Michalogiorgakis and S. S. Pufu, Expanding plasmas and quasinormal modes of anti-de Sitter black holes, JHEP 04 (2007) 080; G. Siopsis, Low frequency quasi-normal modes of AdS black holes, JHEP 05 (2007) 042.
  * [22] V. Cardoso, J. Natario, and R. Schiappa, Asymptotic quasinormal frequencies for black holes in nonasymptotically flat space-times, J. Math. Phys. 45 (2004) 4698.
  * [23] L. Motl, An analytical computation of asymptotic Schwarzschild quasinormal frequencies, Adv. Theoret. Math. Phys. 6 (2003) 1135.
  * [24] L. Motl and A. Neitzke, Advances in Theoretical and Mathematical Physics, Adv. Theoret. Math. Phys. 7 (2003) 307.
  * [25] J. Natario and R. Schiappa, On the classification of asymptotic quasinormal frequencies for d-dimensional black holes and quantum gravity, Adv. Theoret. Math. Phys. 8 (2004) 1001.
  * [26] A. Ghosh, S. Shankaranarayanan, and S. Das, High frequency quasi-normal modes for black holes with generic singularities: II. Asymptotically non-flat spacetimes, Class. Quant. Grav. 23 (2006) 1851.
  * [27] R.G. Daghigh and M.D. Green, A detailed analytic study of the asymptotic quasinormal modes of Schwarzschild-anti de Sitter black holes, Class. Quant. Grav. 26 (2009) 125017.
  * [28] N. Andersson and C. J. Howls, The asymptotic quasinormal mode spectrum of non-rotating black holes, Class. Quant. Grav. 21 (2004) 1623.
  * [29] R.G. Daghigh, The highly real quasinormal modes of Schwarzschild-Anti de Sitter black holes, JHEP 04 (2009) 045.
  * [30] J. Morgan, V. Cardoso, A.S. Miranda, C. Molina, and V.T. Zanchin, Gravitational quasinormal modes of AdS black branes in d spacetime dimensions, JHEP 09 (2009) 117.
  * [31] J. Wang, Z. Chen, X. Tong, Q. Pan, J. Jing, Bifurcation of the Maxwell quasinormal spectrum on asymptotically anti-de Sitter black holes, Phys. Rev. D103 (2021) 064079.
  * [32] S. Fortuna, I. Vega, Electromagnetic quasinormal modes of Schwarzschild-anti-de Sitter black holes: bifurcations, spectral similarity, and exact solutions in the large black hole limit, arXiv 2022.12196 [gr-qc].
  * [33] C.-H. Chen, H.-T. Cho, A. Chrysostomou, A.S. Cornell, Asymptotic quasinormal frequencies of different spin fields in d-dimensional spherically-symmetric black holes, Class. Quant. Grav. 39 (2022) 055001.


the of into a an the of on at with for and in into from by without

# Combining Machine Learning and Spectroscopy to Model Reactive Atom + Diatom
Collisions

Juan Carlos San Vicente Veliz  Department of Chemistry, University of Basel,
Klingelbergstrasse 80, CH-4056 Basel, Switzerland  Julian Arnold  Department
of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel,
Switzerland  Raymond J. Bemish  Air Force Research Laboratory, Space Vehicles
Directorate, Kirtland AFB, New Mexico 87117, USA  Markus Meuwly
<EMAIL_ADDRESS>Department of Chemistry, University of Basel,
Klingelbergstrasse 80, CH-4056 Basel, Switzerland

###### Abstract

The prediction of product translational, vibrational, and rotational energy
distributions for arbitrary initial conditions for reactive atom+diatom
collisions is of considerable practical interest in atmospheric re-entry. Due
to the large number of accessible states, determination of the necessary
information from explicit (quasi-classical or quantum) dynamics studies is
impractical. Here, a machine-learned (ML) model based on translational energy
and product vibrational states assigned from a spectroscopic, ro-vibrational
coupled energy expression based on the Dunham expansion is developed and
tested quantitatively. All models considered in this work reproduce final
state distributions determined from quasi-classical trajectory (QCT)
simulations with $R^{2}\sim 0.98$. As a further validation, thermal rates
determined from the machine-learned models agree with those from explicit QCT
simulations and demonstrate that the atomistic details are retained by the
machine learning which makes them suitable for applications in more coarse-
grained simulations. More generally, it is found that ML is suitable for
designing robust and accurate models from mixed computational/experimental
data which may also be of interest in other areas of the physical sciences.

##  1 Introduction

Atom-diatom collisions at high collision energy are complex due to the
multitude of possible ways in which the available energy can be redistributed.
Relevant processes include reactions, exchange of energy into translation,
rotation, and vibration and depending on the energy, low-lying electronic
states. The number of possible states that are accessible increases
exponentially with the available energy. To be able to understand bulk energy
transfer in high-energy processes, the probabilities of these processes need
to be quantified and understood.

For a reactive atom-diatom collision, determining all state-to-state cross
sections is computationally very challenging even when treating the process
within classical mechanics. This is due to the large number of reactant states
$(\sim 10^{4})$ that can combine with any product state (also $\sim 10^{4}$)
which gives rise to $\sim 10^{8}$ state-to-state cross sections. If such cross
sections are determined from quasi-classical trajectory (QCT) simulations,
typically $10^{4}$ to $10^{5}$ simulations are required for converged results
which yields an estimated $10^{12}$ to $10^{13}$ QCT simulations that would
have to be run at a given collision energy. Depending on the range of relevant
collision energies ($\sim 5$ eV in hypersonics) the number of required QCT
simulations can further increase by one or two orders of magnitude. This is
usually not possible nor desirable. With individual QCT calculations taking
roughly one second, even parallelization and Moore’s law scaling for
improvement in processor speed will not provide sufficient speedup for
exhaustive calculations in the next decades. Current computer technology
limits tractable calculations to $10^{8}$ meaning that any complete model for
the probabilities determined from QCT simulations would be supported by one in
$\sim 10^{7}$ outcomes. In other words: either many final state distributions
are unconverged or the final states are not covered at all, or both. This
sparse representation creates a challenge for accurately modeling the
molecular dynamics in a way that is useful to larger scale simulations.
However, because state-to-state cross sections are of paramount importance to
determine rates for the process of interest, alternative ways to address the
problem are required as this information is further used in reaction networks
to model more complicated chemistries.

One application for which this is particularly relevant, is the high
temperature, high enthalpy flow prevalent in hypersonics. It is not uncommon
to find nominal temperatures that exceed 10000 K in shocks and expansions that
cause local non-equilibrium in the rovibrational states of the chemically
active molecules. Following the dynamics and chemical development of rarefied
gas flows employs primarily two strategies. One uses computational fluid
dynamics (CFD)1 which is based on the Navier Stokes formulation of fluid
dynamics and is valid for small Knudsen number (ratio of the mean free path
length to a physical length scale), whereas direct simulation Monte Carlo
(DSMC)2 is more broadly applicable and is also valid for high Knudsen number.
In DSMC, particles (representing the physical atoms and molecules) move and
collide in physical space to which techniques from statistical mechanics can
be applied.3 Particles in DSMC carry information about their position,
velocity, mass, size, and internal state (for molecules) and move in cells
within which they can collide and exchange energy.

DSMC cycles through the following steps: 1.) moving particles over a time step
$\Delta t$ smaller than the local mean free collision time; 2.) moving
particles across cell boundaries or reflecting them at solid boundaries; 3.)
changing their internal states as a consequence of collisions or reactions;
4.) sample average particle information.3 Step 3 is where microscopic
information about thermal reaction rates, state-to-state cross sections, and
vibrational relaxation rates enters. For chemical reactions most often the
“total collision model”, based on a modified Arrhenius equation, is used,4
although more refined models are also available.5 A more detailed and accurate
description is afforded by state-to-state or state-to-distribution models.
However, experimentally, it is very challenging or even impossible to
determine the relevant quantities at sufficiently high temperature. Most cross
sections presently used are derived from chemical kinetics, many of which have
not been measured at and above 3000 K, or can not be measured but are required
at even higher temperatures. Alternatively, the essential information can also
be obtained from quasi-classical trajectory (QCT) simulations using state-of-
the-art potential energy surfaces. Such an approach provides all necessary
information but implicitly assumes that the PESs and the classical dynamics
underlying QCT simulations are meaningful. In the present work information
from experimental spectroscopy is blended into such a model.

A widely used model for describing microscopic details of hypersonic flow
including chemical and relaxation processes is due to Park.6, 7 This approach
considers separate temperatures $T_{\rm v}$ and $T_{\rm t}$ for the vibration
and rotation/translation, respectively. In application to kinetics the
temperature is taken to be the geometric mean of these, the so-called
“$T-T_{v}$ model”. “Chemistry” enters such models through forward and reverse
reaction rates in the law of mass action for interacting chemical species
which are often determined from temperature-dependent (modified) Arrhenius
expressions. Following and extending the approach from Millikan and White,8 an
important intuitive correction established a framework for including
vibrational non-equilibrium in vibrational relaxation.9

One problematic aspect of the Park model is that vibrational energy becomes
“frozen” above the translational energy because vibrational relaxation is only
included within the limits of Landau-Teller theory10, 8, 11 but the
contribution that accounts for removal of vibrational energy due to
dissociation of the product is neglected as is the movement of large amounts
of translational energy to vibrational energy via atom exchange reactions.
Including the contribution arising from dissociation of the diatomic products
was addressed and corrected in a recent kinetic model which, however,
neglected explicit coupling between intramolecular rotation and vibration.12,
13 Also, with larger computational platforms, it has been possible to
investigate the underlying physics on which the Park approach rests, namely
the two-temperature assumption and the preferential dissociation model14 which
assumes that the amount of vibrational energy removed during dissociation is
large. For N2 \+ N and N2 \+ N2 it was shown that the $T-T_{v}$ model predicts
a much faster N2 dissociation for $T\leq 20000$ K than that obtained with
direct molecular simulation (DMS) whereas for $T=30000$ K the two models
agree.15

Machine-learned (ML) models have proven to be effective in predicting product
distributions based on information about initial states $(E_{\rm trans},v,j)$.
For one, a neural network (NN) based state-to-state (STS) model was conceived
that predicts the cross section for a given transition $(E_{\rm
trans},v,j)\rightarrow(E^{\prime}_{\rm trans},v^{\prime},j^{\prime})$. This
was demonstrated for the [NNO] reactive system16 in that not only the cross
sections for transitions that were not part of the training set were correctly
predicted, but also quantities derived from the cross sections – such as the
total thermal rate $k(T)$ – are in very good agreement with those determined
directly from QCT simulations. Conversely, a distribution-to-distribution
(DTD) model is capable of describing the map between initial and final state
distributions.17

From a practical perspective the most useful model is a state-to-distribution
(STD) model from which final state distributions $P(E^{\prime}_{\rm trans})$,
$P(v^{\prime})$, and $P(j^{\prime})$ can be determined for every initial state
$(E_{\rm trans},v,j)$. This is what is required for more coarse-grained
simulations, such as DSMC. The present work presents state-to-distribution
models for the N(4S)+O2(X${}^{3}\Sigma_{\rm g}^{-}$) $\rightarrow$
NO(X${}^{2}\Pi$) +O(3P) reaction based on translational energy $E_{\rm trans}$
(instead of the diatom’s internal energy $E_{\rm int}$, see Ref.18) and final
vibrational state assignment including mechanical ro-vibrational coupling. It
has been suspected earlier18 that including such coupling may benefit model
performance for high $(v,j)$ states. The trained ML models are based on data
used for the earlier STD model18 to allow for direct comparison between the
different approaches.

The present work is structured as follows. First, the methods used are
presented. This is followed by an analysis of the data based on product
vibrational state assignments $v^{\prime}$ using semiclassical mechanics or
from a model Hamiltonian (here a truncated Dunham expansion). Then, the
performance of NN-trained models based on translational energy together with
the two possibilities for assigning final vibrational quantum numbers is
assessed and compared. Finally, thermal rates obtained via the two approaches
are compared.

##  2 Methods

###  2.1 Quasi-Classical Trajectory Simulations and Analysis

The necessary data for ML-based models characterizing atom + diatom collisions
are based on quasi-classical trajectory (QCT) simulations. The QCT data and an
earlier STD model were used and provided a means to compare the different
approaches.18 Additional QCT simulations were run as needed using the same
procedures as discussed before.18

Initial conditions for QCT simulations were generated from semiclassical
quantization.19, 20 Such an approach couples vibration and rotation in the
sense that the initial vibrational state is assigned from the numerical
solution of an integral involving the rotational barrier $j(j+1)/r^{2}$ where
$r$ is the diatomic bond length. After propagation of a specific initial
condition $(E_{\rm trans},v,j)$ for a given impact parameter $b$, the final
states $(v^{\prime},j^{\prime})$ for the diatomic need to be determined and
the final translational energy $E_{\rm trans}^{\prime}$ is obtained from the
relative velocities and reduced masses of the products. Assigning
$(v^{\prime},j^{\prime})$ quantum numbers is based on final momenta and
positions which are transformed to suitable coordinates. The internal final
angular momentum ${\bf j^{\prime}=q^{\prime}\times p^{\prime}}$ for the
product diatomic species is determined from the final position ${\bf
q^{\prime}}$ and momentum ${\bf p^{\prime}}$. Then the quadratic equation

$j^{\prime}=-\frac{1}{2}+\frac{1}{2}\left(1+4\frac{{\bf j^{\prime}}\cdot{\bf
j^{\prime}}}{\hbar^{2}}\right)^{\frac{1}{2}}$ (1)

is solved to determine the rotational quantum number ($j^{\prime}$) as a non-
integer number.19, 20

Using semiclassical (SC) mechanics, the non-integer vibrational quantum number
($v^{\prime}$) of the final diatomic species is calculated according to19, 20

$v_{\rm
SC}^{\prime}=-\frac{1}{2}+\frac{1}{\pi\hbar}\int^{r^{+}}_{r^{-}}\left\\{2\mu\left(E_{\rm
int}-V(r)-\frac{{\bf j}\cdot{\bf
j}}{2mr^{2}}\right)\right\\}^{\frac{1}{2}}dr,$ (2)

where $r$ is the diatomic bond length, $r^{+}$ and $r^{-}$ are the turning
points of the diatomic species on the effective potential with rotational
state $j^{\prime}$ for internal energy $E_{\rm int}^{\prime}$, $\mu$ is the
reduced mass, and $V(r)$ is the potential energy curve of the product diatom.

Integer ro-vibrational quantum numbers are then assigned as the nearest
integers $(v^{\prime},j^{\prime})$ using histogram binning. To ensure
conservation of total energy, the ro-vibrational energy $E_{\rm
int}^{{}^{\prime}}$ is recomputed from semiclassical quantization19, 20 using
the integer quantum numbers $(v^{\prime},j^{\prime})$ and the final
translational energy for the atom+diatom system is adjusted according to
$E_{\rm trans}^{{}^{\prime}}=E_{\rm tot}-E_{\rm int}^{{}^{\prime}}$ where the
final total internal energy is determined from the final momenta and positions
of the two atoms forming the diatomic.

For high $j^{\prime}$ quantum numbers the angular momentum causes the
potential energy surface of the system to distort, changing the characteristic
frequency of the vibrations. A first-order approximation to the energy based
on the Dunham expansion is shown in Eq. 3. To include mechanical coupling
between vibration and rotation for the products, as afforded by a model
Hamiltonian (“MH”), the assignment of the final angular momentum
$j^{{}^{\prime}}$ is retained as in Eq. 1 but $v^{\prime}$ is determined from
solving

$\displaystyle E(v_{\rm MH}^{\prime},j^{\prime})=$
$\displaystyle\omega_{e}(v_{\rm
MH}^{\prime}+1/2)-\frac{\omega_{e}^{2}}{4D_{e}}(v_{\rm
MH}^{\prime}+1/2)^{2}+B_{e}j^{\prime}(j^{\prime}+1)+D_{e}[j^{\prime}(j^{\prime}+1)]^{2}$
(3) $\displaystyle-\alpha_{e}(v_{\rm
MH}^{\prime}+1/2)j^{\prime}(j^{\prime}+1)$

for $v^{\prime}$. Here, $E(v^{\prime},j^{\prime})=E_{\rm tot}-E_{\rm
trans}^{\prime}$, $\omega_{e}$ is the harmonic frequency, $\omega_{e}x_{e}$ is
the first order correction, $B_{e}$ is the rotational constant, $D_{e}$ is the
centrifugal constant, and $\alpha_{e}$ is the ro-vibrational coupling
constant. For NO the data is $\omega_{e}=1904.20$ cm-1,
$\omega_{e}x_{e}=14.08$ cm-1, $B_{e}=1.672$ cm-1, $D_{e}=0.00000054$ cm-1, and
$\alpha_{e}=0.0171$ cm-1.21 The corresponding parameters on the MRCI/aug-cc-
pVTZ curve for NO are $\omega_{e}=1871.0$ cm-1 and $\omega_{e}x_{e}=14.04$
cm-1.

###  2.2 Neural Network

One approach to creating a model that spans all of the possible probabilities
for the outcome of a collision between an atom and a diatom is to use a neural
network (NN) representation. Such an approach was used previously in earlier
STD work.18 This NN consists of seven residual layers with two hidden layers
per residual layer, and uses 11 input and 254 output nodes corresponding to
the 11 features representing the initial reactant state and the 254 amplitudes
characterizing the product state distributions. For training, the NN inputs
were standardized which ensures that the distributions of the transformed
inputs $x_{i}^{\prime}$ over the entire training data are each characterized
by ($\bar{x}^{\prime}_{i}=0$, $\sigma^{\prime}_{i}=1$). The NN outputs were
normalized which ensures that the distributions of the transformed outputs
$x_{i}^{\prime}$ over the entire training data being characterized by
($\bar{x}^{\prime}_{i}=\bar{x}_{i}$, $\sigma^{\prime}_{i}=1$). Standardization
of the data generally yields a faster convergence.22 The loss function was the
root-mean-squared deviation between reference QCT data and model predictions.
As the NN outputs are probabilities, being non-negative even after
normalization, a softplus function is used as an activation function of the
output layer.

The weights and biases of the NN were initialized according to the Glorot
scheme23 and optimized using Adam24 with an exponentially decaying learning
rate. The NN was trained using TensorFlow 25 and the set of weights and biases
resulting in the smallest loss as evaluated on the validation set were
subsequently used for predictions. Overall, final state distributions from
2184 initial conditions on a grid defined by $(0.5\leq E_{\rm trans}\leq 5.0)$
eV with $\Delta E_{\rm trans}=0.5$ eV $(5.0\leq E_{\rm trans}\leq 8.0)$ eV
with $\Delta E_{\rm trans}=1.0$ eV, $v=[0,2,4,6,8,10,12,15,18,21,24,27,30,\\\
34,38]$, and and $0\leq j\leq 225$ with step size $\Delta j=15$ were
generated, of which 7 were excluded due to low reaction probability. This “on-
grid” set of size $N=2184$, of which 7 data points were excluded due to low
reaction probability, was randomly split into $N_{\rm train}=1700$ for
training, $N_{\rm valid}=400$ for validation, and $N_{\rm test}=77$ for
testing.18 The “on-grid” set used for training, validation and test is
distinguished from “off-grid” data which refers to initial conditions for
which at least one of the entries in $(E_{\rm trans},v,j)$ differs from an
“on-grid” initial condition. All NNs in this work were trained on a 3.60 GHz X
8 Intel Core i7-9700K CPU resulting in training times shorter than 3 minutes.
For additional technical details, see Ref.18

##  3 Results

###  3.1 Data Preparation

A first important task in designing a ML model is the selection and
preparation of the data. Among other aspects, the present work assesses
whether a reliable and predictive STD model can be conceived from training on
/ predicting of a) the translational energy $E_{\rm trans}^{\prime}$ and b)
the final vibrational quantum number $v_{\rm MH}^{\prime}$ determined from a
model Hamiltonian. The reason for this is that $E_{\rm trans}^{\prime}$
contains independent physical information from $(v^{\prime},j^{\prime})$
states about the collision system whereas $E^{\prime}_{\rm int}$ and the
$(v^{\prime},j^{\prime})$ states, as used in the previous model,18 are
somewhat redundant. Similarly, while assignment of $v^{\prime}$ from
semiclassical quantization is one possibility, such an approach neglects part
of the mechanical $(v,j)$ coupling. Furthermore, using a model Hamiltonian
(such as a Dunham expansion26 or a Watson Hamiltonian,27 based on coefficients
fit to represent spectroscopically measured line positions or transitions) for
computing $v_{\rm MH}^{\prime}$ includes valuable information from experiment
and makes the model somewhat less dependent on the level of theory at which
the intermolecular interactions have been / can be determined in practice.
This is the reason to explore changes in the final vibrational distribution if
a model Hamiltonian (here the Dunham expression up to first order in coupling
vibration and rotation) is used. Retaining higher order terms is also possible
to further refine the approach but is not expected to fundamentally change the
findings.

$E^{\prime}_{\rm int}$ versus $E_{\rm trans}^{\prime}$: The reason to employ
$E^{\prime}_{\rm int}$ for training the original STD model18 was that the
number of grid points required to faithfully represent $P(E^{\prime}_{\rm
int})$ was smaller than for $P(E_{\rm trans}^{\prime})$ due to the smaller
span of energies and smoother features. Figures 1A and B report $P(E_{\rm
trans}^{\prime})$ and $P(E^{\prime}_{\rm int})$ for all 2184 initial reactant
conditions (“on-grid”) considered. It appears that $P(E_{\rm trans}^{\prime})$
extends to higher energies than $P(E^{\prime}_{\rm int})$. This is illustrated
by considering a few select final state distributions, see panels E and F in
Figure 1.

Figure 1: Final state distributions $P(E_{\rm trans}^{\prime})$ (panel A),
$P(E_{\rm int}^{{}^{\prime}})$ (panel B) $P(v^{{}^{\prime}})$ (panel C) and
$P(j^{{}^{\prime}})$ (panel D) from semiclassical quantization for NO2 from
QCT simulations for all 2184 initial reactant states. Panels E to H report
selected distributions $P(E_{\rm trans}^{\prime})$, $P(E_{\rm
int}^{{}^{\prime}})$, $P(v^{{}^{\prime}})$, and $P(j^{{}^{\prime}})$ to
highlight their different shapes.

Vibrational Quantum Number from Semiclassical Assignment and Ro-vibrational
Energy Expression: Final state distributions $P(v^{\prime})$ and
$P(j^{\prime})$ from using semiclassical quantization in the final state
analysis are reported in Figure 1. Panels C and D show product state
distributions corresponding to all 2184 initial reactant states considered.
For $P(v^{\prime})$ they are found to extend out to $v^{\prime}\leq 45$
whereas for $P(j^{{}^{\prime}})$ the highest final state is $j^{\prime}\sim
240$ with maximum values $v^{\prime}_{\rm max}=47$, $j^{\prime}_{\rm
max}=240$. Figures 1G and H show individual final state distributions and
highlight the various shapes of these distributions encountered depending on
the initial condition.

Next, the distributions from using SC and MH final states $v_{\rm
SC}^{\prime}$ and $v_{\rm MH}^{\prime}$ are presented for all initial
conditions, see Figure 2A. The black symbols denote $P(v_{\rm SC}^{\prime})$
and extend up to $v_{\rm SC}^{\prime}=46$. This compares with a final state
distribution $P(v_{\rm MH}^{\prime})$ from using the mechanically coupled
energy expression for the assignment of the final state (red symbols) which
only extends up to $v_{\rm MH}^{\prime}=36$. The difference (green) between
the two assignment schemes is reported in Figure 2B. For $v^{\prime}\leq 20$
the difference in the population is $\sim 10$ % which increases to
considerably higher values for $v^{\prime}\sim 30$ and above.

Figure 2: Complete distribution for $P(v_{\rm SC}^{\prime})$ (black) and
$P(v_{\rm MH}^{\prime})$ (red) for the set of 2184 initial conditions (Panel
A) and for selected individual contributions (Panel C). The absolute
difference between $P(v_{\rm SC}^{\prime})$ and $P(v_{\rm MH}^{\prime})$ is
given in panels B (for complete set) $\&$ D (for selected individual
distributions). Selected test initial conditions: $E_{\rm trans}=8.0,v=0,j=0$
(solid lines) and $E_{\rm trans}=1.0,v=30,j=0$ (dashed lines).

Panels C and D of Figure 2 compare final state distributions $P(v_{\rm
SC}^{\prime})$ and $P(v_{\rm MH}^{\prime})$ for two specific initial states.
For $E_{\rm trans}=8.0,v=0,j=0$ (solid lines) the two final vibrational state
distributions are almost indistinguishable, see also their absolute difference
(green) in Figure 2D. On the other hand, for initial $E_{\rm
trans}=1.0,v=30,j=0$ (dashed lines) the maximum for $P(v_{\rm SC}^{\prime})$
(black) is at $v_{\rm SC}^{\prime}=30$ which shifts to $v_{\rm
MH}^{\prime}=26$ for $P(v_{\rm MH}^{\prime})$. Evidently, the two
distributions also differ, see blue dashed line in Figure 2D. The general
finding is that including mechanical ro-vibrational coupling in assigning the
final vibrational state $v_{\rm MH}^{\prime}$ populates lower states compared
with an analysis based on a semiclassical approach.

Assignment of the final state to integer values $(v_{\rm
MH}^{\prime},j^{\prime})$ leads to differences in the corresponding internal
energy $E_{\rm int}^{\prime}(v_{\rm MH}^{\prime},j^{\prime})$ from the value
obtained when considering $E_{\rm int}^{\prime}=E_{\rm tot}-E^{\prime}_{\rm
trans}$ based on energy conservation. Hence, the difference $E_{\rm
int}^{\prime}(v_{\rm MH}^{\prime},j^{\prime})-E_{\rm int}^{\prime}$ has to be
redistributed into $E^{\prime}_{\rm trans}$ which leads to $E_{\rm
trans,MH}^{\prime}$ from including mechanical coupling in the rovibrational
energy. Figure 3A compares the final translational energy distributions
$P(E_{\rm trans,SC}^{\prime})$ (black) and $P(E_{\rm trans,MH}^{\prime})$
(red). The overall shapes of the two distributions are comparable but when
considering the absolute difference between the two distributions (Figure 3B,
green) variations up to 50 % are found. For high translational energy the
absolute differences decay to zero.

Figure 3: Complete distribution for $P(E_{\rm trans,SC}^{\prime})$ (black) and
$P(E_{\rm trans,MH}^{\prime})$ (red) for the set of 2184 initial conditions
(Panel A) and for selected individual contributions (Panel C). The absolute
difference between $P(E_{\rm trans,SC}^{\prime})$ and $P(E_{\rm
trans,MH}^{\prime})$ is given in panels B (for complete set) $\&$ D (for
selected individual distributions) for initial conditions $E_{\rm
trans}=8.0,v=0,j=0$ (solid lines) and $E_{\rm trans}=8.0,v=24,j=120$ (dashed
lines).

Focussing on two specific initial conditions ($E_{\rm trans}=8.0,v=0,j=0$ and
$E_{\rm trans}=8.0,v=24,j=120$, see Figure 3C) indicate that the overall
shapes of the final translational energy distributions depend little on
whether mechanical coupling was included or not to determine $v^{\prime}$. For
initial ($E_{\rm trans}=8.0,v=0,j=0$) nonzero probability at low translational
energy starts at a higher value for an analysis based on MH (red) compared
with SC (black). A similar behaviour is seen for the initial condition given
by ($E_{\rm trans}=8.0,v=24,j=120$) (dashed lines in Figure 3C). Absolute
differences are reported in Figure 3D and reach a maximum of $\sim 20$ %.

###  3.2 Trained STD Models

Next, full STD models were trained based on assignment of the final
vibrational state from semiclassical analysis ($v_{\rm SC}^{\prime}$) or using
a first-order Dunham expansion ($v_{\rm MH}^{\prime}$) together with $E_{\rm
trans}$, respectively. In the following, “on-grid” values refer to initial
conditions that were used for training the NN (training, validation, test data
- see Methods) and “off-grid” corresponds to initial conditions that differed
from “on-grid” values in at least one of the initial conditions $(E_{\rm
trans},v,j)$.

Figure 4: Comparison of product state distributions from QCT (symbols) and STD
model evaluation (red solid line) for SC, for initial conditions from the “on-
grid” set (excluding $E_{\rm trans}=0.5$ eV). STD predictions are for the
initial condition for which the prediction is best (top, largest $R^{2}$), is
closest to the average $R^{2}$ (middle), and worst (bottom, lowest $R^{2}$),
respectively. The corresponding $R^{2}$ values for ($P(E^{\prime}_{\rm
trans})$, $P(v^{\prime})$, $P(j^{\prime})$) are: [0.9985,0.9985,0.9995]
(best); [0.9880,0.9954,0.9895] (closest to mean) and [0.9307,0.4082,0.3578]
(worst). The corresponding reaction probabilities $P_{r}=N_{r}/N_{\rm tot}$
are 0.378, 0.217, and 0.064 for the best, average and worst prediction,
respectively. Here, $N_{\rm tot}$ and $N_{r}$ the total number and number of
reactive trajectories, respectively.

The performance of the trained models in terms of $R^{2}$ and RMSD is
summarized in Tables 1 and 2. The analysis is carried out for all on-grid-
values, excluding seven distributions with low probabilities
($P(v^{\prime})<0.005)$), and excluding the lowest translational energy
$E_{\rm trans}=0.5$ eV because the number of QCT simulations may be
insufficient to fully converge these final state distributions. In addition,
this analysis was also carried out for initial conditions off-grid in $(v,j)$
and at fixed $E_{\rm trans}=4.0$ eV. For all degrees of freedom and all test
sets $R^{2}>0.97$ which indicates reliable statistical models. Such $R^{2}$
measures are also comparable to the performance if instead of using $E_{\rm
trans}$ the models are trained on $E_{\rm int}$ for which $R^{2}>0.98$ was
found.18 If simulations at the lowest translational energy (0.5 eV) are
excluded, the performance is somewhat improved for all degrees of freedom.
Finally, whether the final vibrational state $v^{\prime}$ was assigned from a
semiclassical treatment or a mechanically coupled energy expression has
negligible influence on the performance of the models. When using the RMSD as
the measure to compare reference distributions with predictions from the
trained STD models (see Table 2) similar conclusions are drawn as for $R^{2}$.

Table 1: $R^{2}$ between QCT results and the trained STD model based on $E_{\rm trans}$ and either semiclassical or mechanically coupled determination of $v^{\prime}$. On-grid (All) is for all initial conditions from the training, validation and test set for all translational energies considered. The numbers in bracket are from the test set (77 initial conditions) only. On-grid* and Off-grid* are for initial conditions with $E_{\rm trans}=4.0$ eV only and Off-grid* contains all initial conditions for which either $v$, $j$, or both $(v,j)$ differed from the values used for the training, validation and test sets. | On-grid (All) | On-grid (test) | On-grid* | Off-grid*
---|---|---|---|---
$P(E_{\rm trans,SC}^{\prime})$ | 0.9878 | 0.9844 | 0.9949 | 0.9918
$P(E_{\rm trans,MH}^{\prime})$ | 0.9737 | 0.9840 | 0.9959 | 0.9942
$P(v_{\rm SC}^{\prime})$ | 0.9881 | 0.9835 | 0.9949 | 0.9961
$P(v_{\rm MH}^{\prime})$ | 0.9834 | 0.9817 | 0.9966 | 0.9972
$P(j_{\rm SC}^{\prime})$ | 0.9845 | 0.9816 | 0.9940 | 0.9950
$P(j_{\rm MH}^{\prime})$ | 0.9879 | 0.9882 | 0.9950 | 0.9951
“SC” (overall) | 0.9868 | 0.9832 | 0.9946 | 0.9943
“MH” (overall) | 0.9817 | 0.9846 | 0.9959 | 0.9955

An explicit comparison of final state distributions from models based on SC
(Figure 4) and MH (Figure 5) for assigning the final vibrational state is
provided for the best-performing, for an average-performing, and for the least
performing distribution. The QCT final state distributions are shown as green
open circle and the STD predictions as the red solid line. For the best-
performing prediction the STD model closely follows the target data. This is
still the case for a prediction that represents an average performance. It is
also notable that the shape of the distributions for the same degree of
freedom can change appreciably, depending on the initial condition. For the
worst-performing prediction it is found that SC (see Figure 4) leads to a
visually inferior model compared to using MH (Figure 5). In both worst-
performing cases, the number of reactive trajectories $N_{r}$ as a fraction of
the total number of trajectories $N_{\rm tot}$ is only $\sim 5$ % and
increasing the sampling may also lead to better reference data.

Figure 5: Comparison of product state distributions from QCT (symbols) and STD
model evaluation (red solid line) for MH, for initial conditions from the “on-
grid” set (excluding $E_{\rm trans}=0.5$ eV). STD predictions are for the
initial condition for which the prediction is best (top, largest $R^{2}$), is
closest to the average $R^{2}$ (middle), and worst (bottom, lowest $R^{2}$),
respectively. The corresponding $R^{2}$ for ($P(E^{\prime}_{\rm trans})$,
$P(v^{\prime})$, $P(j^{\prime})$) are: [0.9987,0.9987,0.9991](best);
[0.9911,0.9971,0.9911](closest to mean) and [0.9741,0.5043,0.8305](worst). The
corresponding reaction probabilities $P_{r}=N_{r}/N_{\rm tot}$ are 0.365,
0.308, and 0.044 for the best, average and worst prediction, respectively.
Here, $N_{\rm tot}$ and $N_{r}$ the total number and number of reactive
trajectories, respectively.

For a more global characterization of the final states for the entire range of
possible initial states $(v,j)$, 5 independent models were trained from
reference data for SC (Figure 6A) and MH (Figure 6B). The averaged $R^{2}$
over all 5 models using “on-grid” (see Methods; black circles) and “off-grid”
(see caption Figure 6; crosses) data for each of the initial states considered
shows a uniformly high performance for both approaches with typical $R^{2}\sim
0.98$ or better. For the highest $v-$values the absolute difference between
using MH and SC expressions for $v^{\prime}$ becomes most apparent. There,
using a “MH” model appears to slightly outperform an analysis based on a “SC”
energy expression. But the differences are insignificant.

Figure 6: 2D map of the averaged (over 5 independently trained STD models)
$R^{2}$ values between QCT data ($P(E^{\prime}_{\rm
trans}),P(v^{\prime}),P(j^{\prime})$) and STD evaluation for initial
conditions “on-grid” (see Methods, circles) and “off-grid”
($[v=1,3,5,7,9,11,13,17,22,28,32,36],\\\
j=[10,20,25,35,40,50,55,65,70,80,85,100,110,115,125,130,140,145,155,160,170,175,\\\
185,200,205,215,220]$, crosses). Panel A for “SC” and panel B for “MH”
assignment of the final vibrational state. All QCT simulations were carried
out at $E_{\rm trans}=4.0$ eV.

###  3.3 Thermal Rates

It is also of interest to determine thermal rates from the trained models. In
general, such a rate is determined from the reaction probability $P_{r}$
according to

$k(T)=g(T)\sqrt{\frac{8k_{\rm B}T}{\pi\mu}}\pi b^{2}_{\rm max}P_{r}.$ (4)

For QCT simulations based on stratified sampling of the impact parameter $b$,
$P_{r}=\sum_{k=1}^{K}V_{k}\frac{N_{k}^{r}}{N^{\rm tot}_{k}}$ where $N^{r}_{k}$
and $N^{\rm tot}_{k}$ are the number of reactive and total trajectories,
respectively, for stratum $k$ with impact parameter $b_{k}$. The fractional
volumes $V_{k}=\frac{b_{k}^{2}-b_{k-1}^{2}}{b_{\rm max}^{2}}$ of stratum $k$
obey $\sum_{k=1}^{K}V_{k}=1$. For the STD model the reaction probability can
be obtained as $P_{r}=\int_{E=0}^{E_{\rm max}}P(E)dE$ where $E=E_{\rm
trans}^{\prime}$ and $E_{\rm max}=15.0$ eV. For the forward N(4S) +
O2(X${}^{3}\Sigma_{\rm g}^{-}$) $\rightarrow$ NO(X${}^{2}\Pi$) + O(3P)
reaction in the 4A′ electronic state the degeneracy factor $g(T)=1/3$ and
$\mu$ is the reduced mass of the reactants.28 The two approaches are compared
in Figure 7 and favourable agreement is found over a wide temperature range.
Hence, the STD model can also be used to determine macroscopic quantities such
as realistic reaction rates which is essential.

Figure 7: The thermal forward rate $k_{f}$ calculated from QCT (open red
circle) and STD model (solid black line) for the 4A${}^{{}^{\prime}}$ state of
the N(4S) + O2(X${}^{3}\Sigma_{g}^{-}$ ) $\rightarrow$ NO(X${}^{2}\Pi$) +
O(3P) reaction between 1000 and 20000 K. Additional simulations were made with
the STD model based on $E_{\rm trans}$ and “MH” (blue dashed line). The
present rates agree quantitatively with those directly obtained from QCT
simulations, in particular considering the reduction of required compute time
by 7 orders of magnitude or more. Rates for “MH” are expected to change
slightly if higher-order terms in the Dunham expansion are included.
Experimental total forward reaction rate $k_{\rm f}$ (including contributions
from the doublet and the quartet states) are also shown for comparison: (red
triangle)29, (orange triangle)30 and (magenta circle)31.

##  4 Discussion and Conclusion

The present work discusses a rapid, accurate, and physics-based model to
determine final state distributions for given initial conditions $(E_{\rm
trans},v,j)$ for atom+diatom reactions. It is shown that training on $E_{\rm
trans}^{\prime}$ and $v^{\prime}$ from an energy expression including
mechanical coupling (Dunham expression) yields a model performance of
$R^{2}\sim 0.98$ or better compared with rigorous QCT reference simulations.
Models based on $(E_{\rm trans},v,j)$ are more physically meaningful than
those using $(E_{\rm int},v,j)$ because $E_{\rm int}$ and $(v,j)$ are partly
redundant whereas $E_{\rm trans}$ contains new, complementary, and physically
relevant information. The gain in speed compared to explicit QCT simulations
is about 7 orders of magnitude required for inference (i.e., excluding
training the NN which is minutes, and generating the required QCT reference
data for training) assuming that $10^{5}$ QCT trajectories are sufficient for
converged final state distributions. Hence, this is a promising avenue to be
used in more coarse-grained simulations of reaction networks relevant to
hypersonics and combustion.

Including coupling in determining the final $v^{\prime}$ state of the diatomic
by means of a truncated Dunham expansion as a typical model Hamiltonian (MH)
leads to population of lower $v^{\prime}$ states than assignment from
semiclassical mechanics (SC). The overall shapes of the final state
distributions $P(v^{\prime}_{\rm SC})$ and $P(v^{\prime}_{\rm MH})$ do not
change appreciably, see for example Figure 2. However, final vibrational
distributions from SC extend to higher values $v^{\prime}$ than those from
using a Dunham expansion. The influence of this effect is illustrated in
Figure 8 which compares final internal energy distributions from the two
assignment schemes (SC and MH) for all on-grid initial conditions. The
relevant NO dissociation energies from the SC and MH models are also indicated
and allow to determine the fraction of final states that are above this
threshold to be 21 % and 14 %, respectively, for SC and MH.

If larger products than diatomics are produced in such reactions, it is
expected that the MH approach for assigning final vibrational states is more
convenient than an assignment based on semiclassical mechanics. Furthermore,
using effective Hamiltonians to determine internal final states of molecular
fragments opens the possibility to blend accurate spectroscopic information
into such models which makes them less dependent on the level of quantum
chemical theory at which the intermolecular interactions can be practically
described. For example, multi-reference configuration interaction (MRCI)
calculations with extended basis sets (triple zeta and larger) for molecules
with more than 3 atoms can quickly become computationally prohibitive.
However, it should also be noted that both schemes have their advantages and
drawbacks. For experimental observables, such as thermal rates, both models
perform on par, see Figure 7.

Figure 8: Final internal energy distribution $P(E_{\rm int}^{\prime})$
depending on the assignment scheme (SC, MH) used. Distributions are generated
from the complete data set (on-grid and off-grid) used as initial conditions
in the QCT simulations. The dissociation energy for NO (dashed black) is shown
as well (6.27 eV and 6.29 eV for SC and MH, respectively, compared with 6.50
eV from experiment.21) and the percentage of population above dissociation is
21 % and 14 % for SC and MH, respectively. The $D_{e}$ value for MH was
determined from using a scaling factor between experimental $D_{e}$, and
$D_{e}=\omega_{e}^{2}/(4\omega_{e}x_{e})$ using spectroscopic constants.21

The actual population of highly excited states in the product diatom is
particularly relevant for hypersonics because such states can easily
dissociate during the next collision and the ensuing products (free atoms at
high translational energy) induce rich chemistry. In general, highly
vibrationally excited states are correlated with facile dissociation whereas
high $j-$values prevent and delay breakup of the diatomic. Furthermore, the
vibrational coordinate can also lead to electronic state relaxation rate
differences in subsequent reactions.32 Hence, overall, depending on the
assignment scheme used for $(v,j)$ the energy partitioning into other degrees
of freedom (translation, electronic) will be affected, too.

In summary, the present work demonstrates that it is possible to construct
machine-learned models for final state distributions given state-specific
information for the reactants for atom + diatom reactions. The trained models
are accurate, rapid to evaluate, and can be extended as needed for other
applications, e.g. combustion, by supplying suitable training data. It is
anticipated that such approaches are beneficial to perform more coarse-grained
modeling to applications in hypersonic flow and computational combustion
studies.

## Data and Code Availability

Exemplary data sets and code for training and evaluating STD models can be
found at https://github.com/MMunibas/STDMH.

## Supporting Information

The supporting information contains one additional Table.

##  5 Acknowledgments

This work was supported by AFOSR, the Swiss National Science Foundation
through grants 200021-117810, 200020-188724, the NCCR MUST, and the University
of Basel (to MM). J.A. acknowledges financial support from the Swiss National
Science Foundation individual grant (Grant No. 200020-200481). Discussions
with Dr. D. Koner are greatly appreciated.

## References

  * Knight et al. 2012 Knight, D.; Longo, J.; Drikakis, D.; Gaitonde, D.; Lani, A.; Nompelis, I.; Reimann, B.; Walpot, L. Assessment of CFD capability for prediction of hypersonic shock interactions. _Progr. Aerospace Sci._ 2012, _48-49_ , 8–26
  * Boyd and Schwartzentruber 2017 Boyd, I. D.; Schwartzentruber, T. E. _Nonequilibrium Gas Dynamics and Molecular Simulation_ ; Cambridge University Press, New York, 2017
  * Boyd 2015 Boyd, I. D. Computation of hypersonic flows using the direct simulation Monte Carlo method. _J. Spacecr. Rockets_ 2015, _52_ , 38–53
  * Bird 1994 Bird, G. A. _Molecular Gas Dynamics and the Direct Simulation of Gas Flows_ ; 1994
  * Boyd et al. 1997 Boyd, I.; Bose, D.; Candler, G. Monte Carlo modeling of nitric oxide formation based on quasi-classical trajectory calculations. _Phys. Fluids_ 1997, _9_ , 1162–1170
  * Park 1993 Park, C. Review of Chemical-Kinetic Problems of Future NASA Missions. 1. Earth Entries. _J. Thermophys. Heat Transf._ 1993, _7_ , 385–398
  * Park et al. 2001 Park, C.; Jaffe, R.; Partridge, H. Chemical-kinetic parameters of hyperbolic Earth entry. _J. Thermophys. Heat Transf._ 2001, _15_ , 76–90
  * Millikan and White 1963 Millikan, R. C.; White, D. R. Systematics of Vibrational Relaxation. _J. Chem. Phys._ 1963, _39_ , 3209–3213
  * Park et al. 1994 Park, C.; Howe, J. T.; Jaffe, R. L.; Candler, G. V. Review of chemical-kinetic problems of future NASA missions. II-Mars entries. _J. Thermophys. Heat Transf._ 1994, _8_ , 9–23
  * Landau 1936 Landau, L. Theory of sound dispersion. _Physikalische zeitschrift der Sowjetunion_ 1936, _10_ , 34–43
  * Nikitin and Troe 2008 Nikitin, E.; Troe, J. 70 years of Landau–Teller theory for collisional energy transfer. Semiclassical three-dimensional generalizations of the classical collinear model. _Phys. Chem. Chem. Phys._ 2008, _10_ , 1483–1501
  * Singh and Schwartzentruber 2020 Singh, N.; Schwartzentruber, T. Non-Boltzmann vibrational energy distributions and coupling to dissociation rate. _J. Chem. Phys._ 2020, _152_ , 224301
  * Singh and Schwartzentruber 2020 Singh, N.; Schwartzentruber, T. Consistent kinetic–continuum dissociation model I. Kinetic formulation. _J. Chem. Phys._ 2020, _152_ , 224302
  * Marrone and Treanor 1963 Marrone, P. V.; Treanor, C. E. Chemical relaxation with preferential dissociation from excited vibrational levels. _Phys. Fluids_ 1963, _6_ , 1215–1221
  * Valentini et al. 2016 Valentini, P.; Schwartzentruber, T. E.; Bender, J. D.; Candler, G. V. Dynamics of nitrogen dissociation from direct molecular simulation. _Phys. Rev. Fluids_ 2016, _1_
  * Koner et al. 2019 Koner, D.; Unke, O. T.; Boe, K.; Bemish, R. J.; Meuwly, M. Exhaustive state-to-state cross sections for reactive molecular collisions from importance sampling simulation and a neural network representation. _J. Chem. Phys._ 2019, _150_ , 211101
  * Arnold et al. 2020 Arnold, J.; Koner, D.; Käser, S.; Singh, N.; Bemish, R. J.; Meuwly, M. Machine Learning for Observables: Reactant to Product State Distributions for Atom-Diatom Collisions. _J. Phys. Chem. A_ 2020, _124_ , 7177–7190
  * Arnold et al. 2022 Arnold, J.; San Vicente Veliz, J. C.; Koner, D.; Singh, N.; Bemish, R. J.; Meuwly, M. Machine learning product state distributions from initial reactant states for a reactive atom–diatom collision system. _J. Chem. Phys._ 2022, _156_ , 034301
  * Karplus et al. 1965 Karplus, M.; Porter, R. N.; Sharma, R. D. Exchange Reactions with Activation Energy. I. Simple Barrier Potential for (H, H2). _J. Chem. Phys._ 1965, _43_ , 3259–3287
  * Truhlar and Muckerman 1979 Truhlar, D. G.; Muckerman, J. T. In _Atom - Molecule Collision Theory_ ; Bernstein, R. B., Ed.; Springer US, 1979; pp 505–566
  * Herzberg 1950 Herzberg, G. _Molecular spectra and molecular structure: Spectra of diatomic molecule s_ ; Van Nostrand Reinhold, New York, 1950
  * LeCun et al. 2012 LeCun, Y. A.; Bottou, L.; Orr, G. B.; Müller, K.-R. _Neural networks: Tricks of the trade_ ; Springer, 2012; Vol. 7700; p 948
  * Glorot and Bengio 2010 Glorot, X.; Bengio, Y. Understanding the difficulty of training deep feedforward neural networks. Proceedings of the 13th International Conference on Artificial Intelligence and Statistics. 2010; pp 249–256
  * Kingma and Ba 2014 Kingma, D.; Ba, J. Adam: A method for stochastic optimization. _arXiv preprint arXiv:1412.6980_ 2014,
  * Abadi et al. 2016 Abadi, M. et al. Tensorflow: Large-scale machine learning on heterogeneous distributed systems. _arXiv preprint arXiv:1603.04467_ 2016,
  * Dunham 1932 Dunham, J. The Wentzel-Brillouin-Kramers method of solving the wave equation. _Phys. Rep._ 1932, _41_ , 713
  * Watson 1968 Watson, J. K. Simplification of the molecular vibration-rotation Hamiltonian. _Mol. Phys._ 1968, _15_ , 479–490
  * San Vicente Veliz et al. 2020 San Vicente Veliz, J. C.; Koner, D.; Schwilk, M.; Bemish, R. J.; Meuwly, M. The N(4S) + O2(X ${}^{3}\Sigma_{\rm g}^{-}$) $\leftrightarrow$ O(3P) + NO(X ${}^{2}\Pi$) reaction: thermal and vibrational relaxation rates for the 2A’, 4A’ and 2A” states. _Phys. Chem. Chem. Phys._ 2020, _22_ , 3927–3939
  * Kaufman and Decker 1959 Kaufman, F.; Decker, L. J. _7th Symp. (Int.) Combustion_ 1959, 57
  * Livesey et al. 1971 Livesey, J. B.; Roberts, A. L.; Williams, A. The Formation of Oxides of Nitrogen in some Oxy-Propane Flames. _Combust. Sci. Technol._ 1971, _4_ , 9–15
  * Wilson 1967 Wilson, W. Rate constant for the reaction N + O2 $\rightarrow$ NO + O. _J. Chem. Phys._ 1967, _46_ , 2017–2018
  * Andrienko and Boyd 2016 Andrienko, D.; Boyd, I. D. Vibrational relaxation and dissociation in O2–O mixtures. 46th AIAA Thermophysics Conference. 2016; p 4021

##  6 Table of Contents Graphics

##  7 Supporting Information:

Table 2: RMSD between QCT results and the trained STD model based on $E_{\rm trans}$ and either semiclassical or mechanically coupled determination of $v^{\prime}$. On-grid (All) is for all initial conditions from the training, validation and test set for all translational energies considered. The numbers in bracket are from the test set (77 initial conditions) only. On-grid* and Off-grid* are for initial conditions with $E_{\rm trans}=4.0$ eV only and Off-grid* contains all initial conditions for which either $v$, $j$, or both $(v,j)$ differed from the values used for the training, validation and test sets. | On-grid (All) | On-grid (test) | On-grid* | Off-grid*
---|---|---|---|---
$P(E_{\rm trans,SC}^{\prime})$ | 0.0110 | 0.0134 | 0.0053 | 0.0080
$P(E_{\rm trans,MH}^{\prime})$ | 0.0143 | 0.0151 | 0.0056 | 0.0081
$P(v_{\rm SC}^{\prime})$ | 0.0015 | 0.0019 | 0.0011 | 0.0011
$P(v_{\rm MH}^{\prime})$ | 0.0019 | 0.0025 | 0.0011 | 0.0010
$P(j_{\rm SC}^{\prime})$ | 0.0002 | 0.0003 | 0.0002 | 0.0002
$P(j_{\rm MH}^{\prime})$ | 0.0002 | 0.0003 | 0.0002 | 0.0002
“SC” (overall) | 0.0042 | 0.0052 | 0.0022 | 0.0031
“MH” (overall) | 0.0055 | 0.0059 | 0.0028 | 0.0031


$\lor$ !Olive marble (statement 2) 4\. Eliminate olive marble: !Olive marble
(from statement 2) 5\. Given: Gray marble $\leftrightarrow$ Yellow marble
(statement 1) 6\. Eliminate gray marble: !Gray marble (from statement 1) 7\.
Combine 4 and 6: Yellow marble $\rightarrow$ !Gray marble (from statements 1
and 2) 8\. Combine 7 and 3: Yellow marble $\rightarrow$ Black marble (from
statements 2, 3, and 7) 9\. Conclusion: True. These steps show that the
conclusion follows logically from the given statements. Annotator Comments:
Placeholder Placeholder Reformulation of statement 1. Placeholder Placeholder
Reformulation of statement 2. Placeholder Placeholder Reformulation of
statement 3. Placeholder Placeholder Repetition of conclusion. Placeholder
Chain construction (invalid). Placeholder Placeholder Y $\rightarrow\neg$ O
Placeholder $\neg$ O $\rightarrow$ B (invalid inference). Placeholder Y
$\rightarrow$ B Final answer (incorrect). Symbolic Strategy (invalid).
Translate statement 1. Translate statement 2. (invalid) Translate statement 3.
Translate conclusion. Placeholder Repeat formal conclusion. Repeat formal
statement 3. Repeat formal statement 2. Inference. Repeat formal statement 1.
Invalid inference. Placeholder Placeholder Final answer (incorrect).
Conclusion (invalid). Figure 19: The response (lower left box) of
Zephyr-7B-$\beta$ to problem 6 (top box) of the problem set, illustrating
yellow chain construction and the symbolic strategy. The model first
constructs a chain of conditionals to prove the validity of the conclusion,
linking relevant entities in premise two and three of the problem statement.
Subsequently, the model “explains” its reasoning by employing the symbolic
strategy, converting statements into formal logic and operating on them. Note
that the model makes several logical errors on its way to prove the logical
validity of the final conclusion.