# Graph spatial sampling

_Li-Chun Zhang_ _Statistisk sentralbyrå, Norway_ _University of Southampton,
UK (L.Zhang@soton.ac.uk)_

###### Abstract

We develop lagged Metropolis-Hastings walk for sampling from simple undirected
graphs according to given stationary sampling probabilities. It is explained
how the technique can be applied together with designed graphs for sampling of
units-in-space. We illustrate that the proposed graph spatial sampling
approach can be more flexible for improving the design efficiency compared to
the existing spatial sampling methods.

_Key words:_ graph sampling, random tessellation, local pivotal method,
spatial trend

##  1 Introduction

Denote by $U=\\{1,...,N\\}$ the population of units-in-space, or simply
_units_. By _graph spatial sampling (GSS)_ , one would first _design_ a graph
$G=(U,A)$ and then sample from $G$ – hence its node set $U$ – by graph
sampling methods (Zhang, 2022; Zhang and Patone, 2017). The key idea is to
sensibly introduce the edge set $A$, for which we consider only simple
undirected graphs in this paper, in order to achieve certain desired spatial
properties.

For instance, many spatial sampling methods aim to drastically reduce the
chance of sampling contiguous (or nearby) units, compared to directly sampling
from $U$ by non-spatial methods. To illustrate the idea in terms of GSS, three
graphs are given in Figure 1 for sampling 2 out of 9 spatial units given as
the nodes in $G$, where the edges defining the adjacency among the nodes are
introduced in various ways. Depending on a chosen graph, one can employ
different means for reducing the chance of selecting two contiguous units.

Figure 1: $G_{1}$, adjacency among all contiguous units; $G_{2}$, circular
adjacency mostly among contiguous units; $G_{3}$, all contiguous units are
non-adjacent

First, in the popular generalised random tessellation stratified (GRTS) method
(Stevens and Olsen, 2004), systematic sampling is applied to the units in $U$
arranged on a path of $N-1$ edges which is a special graph $G=(U,A)$. More
generally, the nodes can be arranged in a circle, as in the graph $G_{2}$ in
Figure 1. One can select a systematic sample of size 2 along the circle either
clockwise or anti-clockwise. One can obtain an $(N-1)$-path for GRTS design by
deleting one edge from $G_{2}$; however, it would then be impossible to select
a systematic sample that is always of the size 2. Thus, the approach of GSS
encompasses GRTS sampling.

Moreover, instead of tessellation, one can also consider the graph $G_{3}$ in
Figure 1, where none of the contiguous units are adjacent. As it will be
explained later, one can apply random walk without backtracking in $G_{3}$ and
take as the sample the nodes visited by two successive steps of the walk at
equilibrium, which are never contiguous.

The local pivotal method (LPM) by Grafström et al. (2012) is another popular
spatial sampling method, which can be applied to the graph $G_{1}$ in Figure
1, where the contiguous units (adjacent in $G_{1}$) constitute the nearest
units. As will be shown later, the LPM greatly reduces the chance of selecting
an adjacent pair of nodes in $G_{1}$, compared to random sampling from $U$
directly; whereas GSS from $G_{3}$ avoids this altogether as described above.

Notice that all the designs above are immeasurable. The GSS design using
$G_{3}$ can be made measurable by allowing ‘random jumps’ in addition, as will
be explained later. However, contiguous units can then be selected, e.g. by a
random jump from 1 to 2. Generally, measurability is not considered a priority
in spatial sampling for the sake of improved efficiency, but it does create a
problem for variance estimation when the sampling is without replacement.

Below, we first develop a general technique of walk sampling from graphs in
Section 2, which achieves the desired sampling probabilities. In Section 3, we
explain and illustrate how graph sampling can provide a more flexible approach
to spatial sampling, compared to the existing popular methods. A discussion of
some future topics is given in Section 4.

##  2 Lagged Metropolis-Hastings walk

Given $G=(U,A)$, let $a_{ij}=1$ if $(ij)\in A$ and 0 otherwise. Let
$d_{i}=\sum_{j\in U}a_{ij}$ be the degree of node $i$. We assume $d_{i}\geq 2$
for all the nodes in $G$ and there are no loops, such as the case with all the
graphs in Figure 1. At discrete time step $t>1$, let $X_{t}$ denote the state
(i.e. the current node) of a lagged Metropolis-Hastings walk (LMHW) in $G$,
where $(X_{0},X_{1})$ are the two initial states, and the LMHW transition
probability is given by

$\displaystyle p_{(ih)j}\coloneqq\Pr(X_{t+1}=j\mid X_{t}=h,X_{t-1}=i)$
$\displaystyle=\frac{ru_{j}}{d_{h}+r}+\mathbb{I}(j=i)\frac{wa_{hj}}{d_{h}+r}\min\Big{(}\frac{u_{j}}{u_{h}},1\Big{)}$
$\displaystyle +\mathbb{I}(j\neq
i)\frac{d_{h}-wa_{ih}}{d_{h}+r}\big{(}\frac{a_{hj}}{d_{h}-a_{ih}}\big{)}\min\Big{(}\frac{u_{j}}{u_{h}},1\Big{)}$
(1)

where $\bm{u} =(u_{1},...,u_{N})$ is a positive _preference vector_ satisfying
$\sum_{i\in U} u_{i}=1$. That is, the walk either jumps randomly from
$X_{t}=h$ to any node $j$ (in $U$) with the probability $u_{j}r/(d_{h}+r)$, or
it moves to an adjacent node $j$ with the probability $d_{h}/(d_{h}+r)$. In
the latter case, it can either backtrack to the _previous_ $X_{t-1}=i$ (if
adjacent) with a probability regulated by $w$ or move to _another_ adjacent
node, both of which are subject to a Metropolis-Hastings (MH) acceptance
mechanism, hence the term LMHW. There would be no backtracking under LMHW if
$w=0$, and random jumps are disallowed if $r=0$ as long as $G$ is connected.

The LMHW (1) generalises the lagged random walk (LRW) proposed by Zhang (2021)
where $u_{i}\equiv 1$, i.e. without the MH mechanism. Moreover, in the case of
$w=1$, the LRW reduces to targeted random walk (TRW) of Avrachenkov et al.
(2010), under which there is no difference between the previous state
$X_{t-1}$ (if adjacent) and the other nodes adjacent to $X_{t}$. Thompson
(2006) considers random walk (not lagged) subject to MH acceptance mechanism.

The process $\\{X_{t}:t\geq 0\\}$ is non-Markovian if $w<1$. Let
$\bm{x}_{t}=(X_{t-1},X_{t})$ for $t\geq 1$. Given any initial
$\bm{x}_{1}=(X_{0},X_{1})$, LMHW (1) generates a Markov chain
$\\{\bm{x}_{t}:t\geq 1\\}$, since

$\Pr(\bm{x}_{t+1}|\bm{x}_{t},...,\bm{x}_{1})=\Pr(\bm{x}_{t+1}|\bm{x}_{t})=\Pr(X_{t+1}|\bm{x}_{t})$
(2)

It is irreducible, if $G$ is connected or if random jumps are allowed
generally, such that there exists a unique stationary distribution,
$\Pr\big{(}\bm{x}_{t}=(h,j)\big{)}$, which is given by

$p_{(hj)}=\sum_{i\in U}p_{(ih)}p_{(ih)j}$ (3)

A unique stationary distribution of $X_{t}$ follows, which satisfies the
_mixed_ equation

$p_{h}\coloneqq\Pr(X_{t}=h)=\sum_{i\in
U}\Pr\big{(}\bm{x}_{t}=(i,h)\big{)}=p_{h}p_{hh}
+\sum_{\begin{subarray}{c}\bm{x}=(i,h)\\\
i\in\nu_{h}\end{subarray}}p_{\bm{x}}+\sum_{\begin{subarray}{c}i\not\in\nu_{h}\\\
i\neq h\end{subarray}}\frac{p_{i}ru_{h}}{d_{i}+r}$ (4)

where $p_{hh}\coloneqq\Pr(X_{t}=h|X_{t-1}=h)$ at equilibrium, and
$\nu_{h}=\\{i:a_{ih}=1,i\in U\\}$ is the _neighbourhood_ of $h$ (containing
its adjacent nodes), and a transition from $h$ to any node outside $\nu_{h}$
can only be accomplished by a random jump. Notice that $\bm{x}_{t}=(h,h)$ for
$t>1$ is possible if a random jump from $h$ lands on $h$ itself, or if a
proposed move to an adjacent node is rejected. Appendix A gives a proof that
the stationary probability is given by

$p_{h}\propto(d_{h}+r)u_{h}$ (5)

We have $p_{h}=\pi_{h}/n$ if $(d_{h}+r)u_{h}\propto\pi_{h}$ for all $h\in U$,
where $\\{\pi_{i}:i\in U\\}$ are the given sample inclusion probabilities by
$\pi$ps sampling without replacement from $U$, where $\sum_{i\in U}\pi_{i}=n$.

##  3 GSS by LMHW

###  3.1 Equal-probability spatial sampling without replacement

_Equal-probability spatial sampling without replacement (EpSSWoR)_ of sample
size $n$ has the same sample inclusion probability $n/N$ as simple random
sampling without replacement (SRSWoR) from $U$ directly, but the second-order
inclusion probability of SRSWoR (which is the same for any pair of distinct
units) can be modified to achieve desired spatial properties.

LMHW can yield a GSS method for EpSSWoR. To ensure sampling without
replacement, it is necessary to set $(r,w)=(0,0)$. In addition, to achieve
equal probability (5) at equilibrium and to remove the possibility of
rejecting any proposed transition to an adjacent node, devise a connected
_2-regular_ graph $G$, where $d_{i}\equiv 2$, and set $u_{i}\equiv 1$. Now
that $d_{i}\equiv 2$ and there is no backtracking, $X_{t+1}$ must be an
adjacent node to $X_{t}$ which is not visited by the walk in the previous
$N-1$ steps. It follows that any $n$-sequence of states
$(X_{t},X_{t+1},...,X_{t+n-1})$ is a sample of $n$ distinct units from $U$,
where $p_{h}\equiv 1/N$ for any node $h$ in the sequence.

####  3.1.1 Illustration

Now, one can devise the 2-regular graph $G$ according to the desirable spatial
sampling properties. Suppose one would like to reduce the probability of
selecting contiguous units, denoted by $\xi$. Consider below the spatial
population $U$ in Figure 1 for an illustration.

First, let the sample size be 2. There are 12 contiguous pairs (as can be seen
in $G_{1}$) out of 36 possible pairs of distinct units, such that $\xi=1/3$
under SRSWoR from $U$ directly. Simulations of the LPM1 version of LPM
(Grafström et al., 2012) from $G_{1}$ yields $\xi=0.116$. The GRTS method
cannot ensure the sample size is always 2 in this case. For GSS by clockwise
systematic sampling from $G_{2}$ in Figure 1, there are 9 systematic samples
of size 2, i.e. $\\{1,6\\}$, $\\{5,9\\}$, …, $\\{7,2\\}$ and $\\{4,3\\}$,
where only $\\{8,5\\}$ contains contiguous units, such that $\xi=1/9=0.111$.

Meanwhile, for EpSSWoR by LMHW (1), one can use a 2-regular graph that does
not contain any edge connecting two contiguous units in $U$. There are many
such graphs, two of which are shown in Figure 2. Using such a 2-regular graph,
we obtain $\xi=0$ by construction.

Figure 2: Two 2-regular graphs for EpSSWoR from 9 units-in-space

Next, let the sample size be 3. There are $84$ distinct samples by SRSWoR from
$U$, where 22 of them do not contain any contiguous units, such that
$\xi=62/84=0.738$. Simulations of the LPM1 from $G_{1}$ yield $\xi=0.478$.
There are 3 systematic samples by the GRTS method, because $N/n=3$. For
instance, let the 8-path be given by $(1,4,7,8,9,6,3,2,5)$, i.e. removing the
edge between 1 and 5 in $G_{2}$, the three samples are $\\{1,8,3\\}$,
$\\{4,9,2\\}$ and $\\{7,6,5\\}$, such that $\xi=1/3$. The same holds for GSS
by systematic sampling from $G_{2}$ in Figure 1.

Meanwhile, for EpSSWoR by LMHW from $G_{4}$ (Figure 2), there are 9 distinct
samples, where contiguous units are present in the four samples $\\{3,9,2\\}$,
$\\{9,2,8\\}$, $\\{4,6,7\\}$ and $\\{6,7,5\\}$, such that $\xi=4/9=0.444$.
However, suppose one instead adopts $G_{5}$ in Figure 2, then only the sample
$\\{6,2,9\\}$ contains contiguous units, such that $\xi=1/9=0.111$.

####  3.1.2 Implementation

One can construct 2-regular graphs by means of the recursive partitions used
for GRTS design. The example of Stevens and Olsen (2004) is given in Figure 3
(left), containing 64 units divided into 16 parts. Instead of connecting the
nearby units as in the GRTS design, one can connect the more distant units, as
illustrated for the 0-units in Figure 3 (right), which are non-contiguous due
to the other units 1, 2, 3. The starting and end nodes are underlined in
Figure 3. Without loss of generality, suppose the one in the bottom-left
corner is the end node. One can connect it to one of the 1-nodes that is not
contiguous to the starting 0-node. Similarly for the other units 1, 2, 3.
Finally, since the 3-nodes and 0-nodes are never contiguous here, connecting
the end 3-node and the starting 0-node yields a non-contiguous 2-regular graph
$G$.

Figure 3: Recursive partition for GRTS design (left), 0-nodes in a 2-regular
graph (right)

Numerous 2-regular graphs can be devised like this; denote the collection of
them by $\mathbb{G}$. For each $G$ in $\mathbb{G}$, let $\Omega_{G}$ contain
the $N$ possible samples of the given size $n$. One can either calculate or
simulate some _design measure_ over $\Omega_{G}$, denoted by $\tau_{G}$, such
as $\xi$ above or the expected sample spatial balance (SSB) measure of Stevens
and Olsen (2004) or the sampling variance given $\\{y_{i}:i\in U\\}$ generated
by a suitable spatial population model. One can explore $(G,\tau_{G})$ over
$G\in\mathbb{G}$ and choose the graph $G$ that has the best design measure
$\tau_{G}$.

Given EpSSWoR by GGS-LMHW, a design-unbiased estimator of the population total

$Y=\sum_{i\in U}y_{i}$

is the Horvitz-Thompson estimator. However, unbiased estimation of its
sampling variance is impossible as long as the sampling design is
immeasurable.

###  3.2 Unequal-probability spatial sampling

Let the graph $G$ for GSS be connected so that random jumps are unnecessary
and set $r=0$. To further reduce the chance of selecting the same node more
than once, set $w=0$ so that $\Pr(X_{t+1}=X_{t-1})=0$. Finally, set the
preference vector $\bm{u}$, such that

$p_{i}=d_{i}u_{i}\eta=\frac{\pi_{i}}{n}\qquad\text{and}\qquad
u_{i}=\frac{\pi_{i}}{nd_{i}\eta}\qquad\text{and}\qquad\eta=\frac{1}{n}\sum_{i\in
U}\frac{\pi_{i}}{d_{i}}$

at equilibrium for any $i\in U$. Now, as long as $u_{i}$ is not a constant
over $U$, one cannot avoid selecting some node more than once due to the
rejections. Since the inclusion probability of any given node in
$\\{X_{t+1},X_{t+2},...,X_{t+m}\\}$ becomes intractable as $m$ increases, we
use the stationary sampling probabilities $p_{i}$ for unbiased estimation of
the total $Y$.

####  3.2.1 Based on $m$-sequence at equilibrium

Let $(X_{t+1},X_{t+2},...,X_{t+m})$ be a sequence of $m$ states from the LMHW
at equilibrium. An unbiased estimator of $Y$ can be given in various forms as

$\hat{Y}_{W}=\frac{1}{m}
\sum_{j=1}^{m}\frac{y_{X_{t+j}}}{p_{X_{t+j}}}=\frac{n}{m}
\sum_{j=1}^{m}\frac{y_{X_{t+j}}}{\pi_{X_{t+j}}}=\frac{1}{m}\sum_{j=1}^{m}
\sum_{i\in U}\frac{y_{X_{t+j}}}{p_{X_{t+j}}}\mathbb{I}(X_{t+j}=i)$ (6)

This includes $\pi_{i}\equiv n/N$ as a special case. We have
$E(\hat{Y}_{W}|m)=Y$ because $\Pr\left[\mathbb{I}(X_{t+j}=i)\right]=p_{i}$
given any time step $t+j$ under LMHW at equilibrium. Notice that the $m$
states in (6) are not independent of each other, although the first expression
of $\hat{Y}_{W}$ looks the same as the Hansen-Hurwitz estimator under sampling
with replacement.

Let $(X_{t+a},...,X_{t+b})$ be a subsequence of $(X_{t+1},...,X_{t+m})$, where
$1\leq a\leq b\leq m$. It is said to a _tie_ of order $b-a+1$ for some node
$h\in U$, denoted by $\kappa_{a,b}^{h}$, if

$X_{t+a-1}=i\neq h,~{} X_{t+a}=\cdots=X_{t+b}=h,~{} X_{t+b+1}=j\neq h$

Appendix B gives an unbiased estimator of $Y$ based on all the ties in an
$m$-sequence. However, simulations suggest that it is typically less efficient
than the simpler estimator (6).

####  3.2.2 Illustration

Consider the following stylised examples of spatial populations $\bm{y}_{U}$
for Figure 1:

$\displaystyle\left[\begin{array}[]{ccc}1&2&1\\\ 2&3&2\\\
1&2&1\end{array}\right]\qquad\left[\begin{array}[]{ccc}3&2.5&2\\\ 2.5&2&1.5\\\
2&1.5&1\end{array}\right]\qquad\left[\begin{array}[]{ccc}3&2&1\\\ 2&1&2\\\
1&2&3\end{array}\right]\qquad\left[\begin{array}[]{ccc}3&2&3\\\ 2&1&2\\\
3&2&3\end{array}\right]$ $\displaystyle\hskip 11.38109pt \text{Centre}\hskip
62.59605pt \text{Corner}\hskip 65.44133pt \text{Polar}\hskip 54.06023pt
\text{Vortex}$

Let $\bm{\pi}_{U}$ be all equal if $\pi_{5}/\pi_{i}\equiv 1$, or unequal if
$\pi_{5}/\pi_{i}\equiv 2$ or $\pi_{5}/\pi_{i}\equiv 0.5$, for all $i\neq 5$.
We apply the LPM1 (Grafström et al., 2012) to select a sample of size 2, as
well as LMHW sampling from each of $G_{1}$ \- $G_{5}$ with $m=2$ and $r=w=0$.
Simulations yield the relative efficiency (RE) of a given sampling method
against SRSWoR with $n=2$.

Table 1: RE of LPM1 ($n=2$) or GSS ($m=2$) from $G_{1}$ \- $G_{5}$

$\bm{y}_{U}$ | $\pi_{5}/\pi_{i}$ | LPM1 | $G_{1}$ | $G_{2}$ | $G_{3}$ | $G_{4}$ | $G_{5}$
---|---|---|---|---|---|---|---
Centre | 1 | 0.92 | 1.27 | 0.58 | 0.57 | 0.87 | 0.85
|  | (0.18) |  |  |  |
2 | 0.62 | 0.32 | 0.08 | 0.56 | 0.84 | 0.82
|  | (0.20) | (0.1) | (0.1) | (0.1) | (0.1)
Corner | 1 | 0.71 | 1.91 | 1.72 | 0.77 | 0.66 | 0.68
Polar | 1 | 1.04 | 1.40 | 0.88 | 0.89 | 0.65 | 0.67
0.5 | 0.81 | 1.08 | 0.96 | 0.81 | 0.66 | 0.67
|  | (0.27) | (0.07) | (0.07) | (0.07) | (0.07)
Vortex | 1 | 0.91 | 1.27 | 0.58 | 0.58 | 0.88 | 0.84
0.5 | 0.58 | 0.45 | 0.18 | 0.46 | 0.74 | 0.73

Note: Positive $\Pr(n=1)$ by GSS given in parentheses

The results in Table 1 are based on $10^{4}$ simulations of each sampling
method given $\bm{y}_{U}$. It is possible here that $\Pr(n=1)=\sum_{h\in
U}p_{(hh)}>0$ under LMHW sampling due to the rejected moves, where the
probability depends only on $(G,\bm{\pi}_{U})$ but not $\bm{y}_{U}$. Setting
$\pi_{5}/\pi_{i}\equiv 2$ can only be plausible for the centre $\bm{y}_{U}$,
similarly as setting $\pi_{5}/\pi_{i}\equiv 0.5$ for the polar or vortex
$\bm{y}_{U}$. Given $\pi_{5}/\pi_{i}\equiv 2$ for the centre $\bm{y}_{U}$, the
two GSS methods using $G_{1}$ or $G_{2}$ select mostly contiguous units, both
of which are actually more efficient than the other methods that aim to avoid
selecting contiguous units; similarly given $\pi_{5}/\pi_{i}\equiv 0.5$ for
the vortex $\bm{y}_{U}$. This serves as a reminder not to treat any particular
sample spatial balance property as a panacea for design efficiency, without
taking into account the spatial distribution of $\bm{y}_{U}$.

For equal-probability sampling across the 4 populations, although the LPM1
improves upon SRSWoR except in one case, it is always dominated by some (or
all) of the GSS methods using $G_{3}$ \- $G_{5}$. Among these GSS methods,
using $G_{4}$ or $G_{5}$ yields essentially the same RE here, using $G_{3}$ is
more efficient for the centre and vortex $\bm{y}_{U}$ but not otherwise. It is
thus important to consider different graph designs for different spatial
distributions of $\bm{y}_{U}$.

###  3.3 Comparison of designs by simulation

Grafström et al. (2012) suggest the LPM can yield large gains over the GRTS
method for populations with smooth spatial trends, particularly in their
Example 5 with $400$ units evenly spread over the unit square and
$\pi_{i}\equiv n/N$ for $i\in U$, where the $y$-values are given by

$\text{sinTrend:}\quad y(x_{1},x_{2})=3(x_{1}+x_{2})+\sin\\{6(x_{1}+x_{2})\\}$

and $(x_{1},x_{2})$ are the coordinates. We consider also the four types of
$\bm{y}_{U}$ in Section 3.2.2 for this $U$, where $0.5\leq y_{i}\leq 5$ for
$i\in U$, which is about the same range as the sinTrend $\bm{y}_{U}$ above.

Table 2: RE and ESSB of LPM1, $G_{6}$SS or $G_{7}$SS |  | RE |
---|---|---|---
Sample | Method | sinTrend | Centre | Corner | Polar | Vortex | ESSB
$n=16$ | LPM1 | 0.151 | 0.248 | 0.127 | 0.221 | 0.244 | 0.080
$G_{6}$SS | 0.561 | 0.025 | 0.801 | 0.060 | 0.025 | 0.079
$G_{7}$SS | 0.044 | 1.371 | 0.016 | 1.047 | 1.362 | 0.192
$n=32$ | LPM1 | 0.090 | 0.147 | 0.072 | 0.132 | 0.150 | 0.074
$G_{6}$SS | 0.925 | 0.020 | 1.288 | 0.077 | 0.020 | 0.111
$G_{7}$SS | 0.027 | 1.489 | 0.009 | 1.227 | 1.543 | 0.238
$n=48$ | LPM1 | 0.067 | 0.111 | 0.053 | 0.098 | 0.114 | 0.079
$G_{6}$SS | 1.138 | 0.015 | 1.595 | 0.085 | 0.015 | 0.154
$G_{7}$SS | 0.022 | 1.421 | 0.007 | 1.321 | 1.375 | 0.238

Two 2-regular graphs are used for GSS here. The graph $G_{6}$ follows the
description in Section 3.1.2 (Figure 3), with the $4\times 4$-partition of $U$
and 25 nodes in each part. The graph $G_{7}$ uses the $2\times 2$-partition as
follows. First, index each unit $(x_{1},x_{2})$ as $(r_{1},r_{2})$, where
$r_{1}$ is the rank of $x_{1}$ and $r_{2}$ that of $x_{2}$. Next, each pair of
units $(r_{1},r_{2})$ and $(20-r_{1}+1,20-r_{2}+1)$ are made adjacent, for
$r_{1},r_{2}=1,...,20$, i.e. between top-left and bottom-right parts as well
as between top-right and bottom-left parts. Finally, the units in the top-left
and bottom-left parts are randomly paired to be adjacent, likewise for the
top-right and bottom-right parts.

Table 2 gives the RE-results (each by $10^{4}$ simulations) and the expected
sample spatial balance (ESSB), where the sample size $n\in\\{16,32,48\\}$ as
in Grafström et al. (2012). For any $n$, $G_{6}$SS improves greatly over LPM1
for the Centre and Vortex $\bm{y}_{U}$, whereas $G_{7}$SS does so for the
Corner and sinTrend $\bm{y}_{U}$. For the Polar $\bm{y}_{U}$, the RE is seen
to become closer between LPM1 and $G_{6}$SS as $n$ increases, while both are
considerably more efficient than SRSWoR. Notice that, since the ESSB is a
constant given $n$ here, whichever the spatial population $\bm{y}_{U}$, one
cannot anticipate the design efficiency _only_ based on such a measure.

There exists a trend along $x_{1}+x_{2}$ in both the Corner and sinTrend
$\bm{y}_{U}$, apart from a sinus undulation in the latter. The results suggest
that the merits of $G_{7}$SS vs. LPM for the sinTrend $\bm{y}_{U}$ can be
anticipated based on the Corner $\bm{y}_{U}$. Due to the structural similarity
between the Centre and Vortex $\bm{y}_{U}$, the merits of $G_{6}$SS vs. LPM
for one population can be anticipated from that for the other. The results for
the Polar $\bm{y}_{U}$ suggest there may be room for improving the graph
design for GSS as $n$ increases for this and similar spatial populations.

##  4 Some future topics

Random walk has numerous applications (e.g. Masuda et al., 2017; Brin and
Page, 1998). LMHW offers a more flexible technique, which allows one to choose
the desired stationary probabilities via the preference vector $\bm{u}$ while
controlling the probability of back-tracking by $w$. It can be considered for
many problems beyond spatial sampling.

Both the GRTS method and the LPM can be motivated from the perspective of
improving the expected SSB compared to sampling from $U$ directly. GSS
provides a flexible approach to accommodate the anticipated spatial
distribution of $\bm{y}_{U}$ in addition. It encompasses the GRTS method and,
as illustrated above, suitable graph designs can yield large gains over the
LPM. To facilitate the practice of GSS, one should develop suitable graph
design algorithms that scale as the population size increases, and investigate
their properties for various typical spatial distributions of $\bm{y}_{U}$ in
a more systematic manner.

For spatial sampling without replacement from $U$, variance estimation does
not admit a theoretical solution. For GSS that allows for repeated selection
of a given unit by LMHW, one can initiate multiple independent walks, each
yielding an unbiased estimator (6) — one can use the mean of them to estimate
$Y$ and use the between-walk variance of them for unbiased variance
estimation, which is a standard technique in MCMC.

##  Appendix A Proof of (5)

Under LMHW (1), balanced flows between $\bm{x}_{t}=(i,h)$ and
$\bm{x}_{t+1}=(h,j)$ are the flows over $(X_{t-1},X_{t},X_{t+1})$ in either
direction. To show the values $\\{p_{h}=d_{h}+r:h\in U\\}$ satisfy the
balanced flows at equilibrium, one needs to consider the following situations
I - V.

##### I.

$i=h=j$, which is balanced at equilibrium, where the probability of either
flow (in the opposite order) is equal to
$\Pr\big{(}X_{t+1}=h|\bm{x}_{t}=(h,h)\big{)}$.

##### II.

$i=j\neq h$ and $i\in\nu_{h}$, i.e. $(i=j)$ —$~{}h$. Since, since $i=j$, both
the flows $(i,h,j)$ and $(j,h,i)$ are backtracking (in either direction), the
probability of which is the same by (1), so that they are always balanced.

##### III.

$i\neq j\in\nu_{h}$, i.e. $i$ —$~{}h$ —$~{}j$. Since $i\neq j$, neither
$(i,h,j)$ nor $(j,h,i)$ is backtracking, given which (1) yields

$\displaystyle\sum_{i\in\nu_{h}}\sum_{\begin{subarray}{c}j\in\nu_{h}\\\ j\neq
i\end{subarray}}\Pr\big{(}\bm{x}_{t}=(i,h)\big{)}\Big{(}\frac{d_{h}-w}{d_{h}-1}\Big{)}\frac{\min(u_{j},u_{h})}{(d_{h}+r)u_{h}}$
$\displaystyle\sum_{j\in\nu_{h}}\sum_{\begin{subarray}{c}i\in\nu_{h}\\\ i\neq
j\end{subarray}}\Pr\big{(}\bm{x}_{t}=(j,h)\big{)}\Big{(}\frac{d_{h}-w}{d_{h}-1}\Big{)}\frac{\min(u_{i},u_{h})}{(d_{h}+r)u_{h}}$

as the sums of probabilities in either direction, which are balanced by
symmetry.

##### IV.

$\\{i,j\\}\not\in\nu_{h}$, including $i=h$ or $j=h$, i.e.

$i\quad h\quad j$ or $(i=h)\quad j$ or $i\quad(h=j)$

In the first case, where $i,h,j$ are distinct, any flow $(i,h,j)$ and
$(j,h,i)$ can only take place by random jumps, which are balanced on setting
$p_{i}\equiv(d_{i}+r)u_{i}$, since

$p_{i}\frac{r}{d_{i}+r}u_{h}\frac{r}{d_{h}+r}u_{j}=p_{j}\frac{r}{d_{j}+r}u_{h}\frac{r}{d_{h}+r}u_{i}$

For the other two cases, on noting
$p_{h}p_{hh}=\Pr\big{(}\bm{x}_{t}=(h,h)\big{)}$ by definition and setting
$p_{h}=(d_{h}+r)u_{h}$, we obtain

$\sum_{j\not\in\nu_{h}}p_{h}p_{hh}\frac{ru_{j}}{d_{h}+r}=\sum_{i\not\in\nu_{h}}p_{h}p_{hh}\frac{ru_{i}}{d_{h}+r}$

where the left-hand side is the sum of probabilities in the 2nd case (i.e.
$i=h$) in the direction $(i,h,j)$, and the right-hand side is the sum of
probabilities in the 3rd case (i.e. $h=j$) in the opposite direction
$(j,h,i)$. It follows that these two cases balance out each other.

##### V.

One of $(i,j)$, say, $i$ belongs to $\nu_{h}$ but not the other, including
when $j=h$, i.e.

$i$ —$~{}h\qquad j$ or $i$ —$~{}(h=j)$

A flow in the direction $(i,h,j)$ consists of two parts: (i) $\bm{x}=(i,h)$
where $i\in\nu_{h}$, and (ii) either a random jump from $h$ to any node
outside of $\nu_{h}$ (including $h$) or a proposed move into $\nu_{i}$ is
rejected. Summing the stationary probabilities of all such flows, we have

$\big{(}\sum_{\begin{subarray}{c}\bm{x}=(i,h)\\\
i\in\nu_{h}\end{subarray}}p_{\bm{x}}\big{)}\left(\frac{r}{d_{h}+r}\sum_{j\not\in\nu_{h}}u_{j}+\frac{\sum_{i\in\nu_{h}}u_{h}-\min(u_{i},u_{h})}{(d_{h}+r)u_{h}}\right)$

where the 2nd term in the parentheses corresponding to (ii) is the sum over
all possible moves, including backtracking to $X_{t-1}=i$ and forwarding to
$l\neq X_{t-1}$ for all $i\in\nu_{h}$, which is a constant of $w$. Meanwhile,
a flow in the opposite direction $(j,h,i)$ consists of two parts: (a)
$\bm{x}=(j,h)$ where $j\not\in\nu_{h}$, including $j=h$, and (b) a transition
from $h$ to an adjacent node. Summing the stationary probabilities of all such
flows, we have

$\big{(}p_{h}-\sum_{\begin{subarray}{c}\bm{x}=(i,h)\\\
i\in\nu_{h}\end{subarray}}p_{\bm{x}}\big{)}\left(\frac{r}{d_{h}+r}
\sum_{i\in\nu_{h}}u_{i}+\frac{\sum_{i\in\nu_{h}}\min(u_{i},u_{h})}{(d_{h}+r)u_{h}}\right)\coloneqq\big{(}p_{h}-\sum_{\begin{subarray}{c}\bm{x}=(i,h)\\\
i\in\nu_{h}\end{subarray}}p_{\bm{x}}\big{)}\Delta_{h}$

where the expression in the first pair of parentheses corresponding to (a) is
by definition. To balancing the two groups of flows, we require

$p_{h}\Delta_{h}=\big{(}\sum_{\begin{subarray}{c}\bm{x}=(i,h)\\\
i\in\nu_{h}\end{subarray}}p_{\bm{x}}\big{)}\left(\frac{r}{d_{h}+r}
+\frac{d_{h}u_{h}}{(d_{h}+r)u_{h}}\right)=\sum_{\begin{subarray}{c}\bm{x}=(i,h)\\\
i\in\nu_{h}\end{subarray}}p_{\bm{x}}$ (7)

since $\sum_{i\in U} u_{i}=1$, where the left-hand side in (7) can now be
rewritten as

$\displaystyle
p_{h}\Big{(}1-p_{hh}-\frac{r}{d_{h}+r}\sum_{\begin{subarray}{c}i\not\in\nu_{h}\\\
i\neq
h\end{subarray}}u_{i}\Big{)}=p_{h}-p_{h}p_{hh}-ru_{h}\sum_{\begin{subarray}{c}i\not\in\nu_{h}\\\
i\neq h\end{subarray}}u_{j}$

and the expression in the parentheses is
$1-\sum_{j\not\in\nu_{h}}\Pr(X_{t+1}=j|X_{t}=h)$ by definition, whereas
$p_{h}=(d_{h}+r)u_{h}$ is used to obtain the last term on the right-hand side,
which is equal to the last term on the right-hand side of (4) on setting
$p_{i}\equiv(d_{i}+r)u_{i}$. In other words, setting
$p_{h}\equiv(d_{h}+r)u_{h}$ reduces (7) to (4), by which all the flows are
balanced.

Thus, balanced flows are achieved in all the situations above. This completes
the proof.

##  Appendix B Estimation based on $n_{m}$ ties in $m$-sequence

Let the $m$-sequence $(X_{t+1},...,X_{t+m})$ consist of $n_{m}$ ties, denoted
by $H_{m}=\\{h_{1},...,h_{n_{m}}\\}$, where $n_{m}$ is random and $1\leq
n_{m}\leq m-2$. The stationary probability of a tie $\kappa_{a,b}^{h}$ is
given by

$p_{(\kappa_{ab}^{h})}=\begin{cases} \sum\limits_{\begin{subarray}{c}i,j\in
U\\\ i\neq h,j\neq h\end{subarray}}p_{(ih)}p_{(ih)j}&\text{if } a=b\\\
\sum\limits_{\begin{subarray}{c}i\in U\\\ i\neq
h\end{subarray}}p_{(ih)}p_{(ih)h}p_{(hh)h}^{b-a-1}(1-p_{(hh)h})&\text{if }
a<b\end{cases}$

where $p_{(ih)}$ is given by (3) and all the transition probabilities by (1).
Let $\delta_{ab}=1$ if $(X_{t+a},...,X_{t+b})$ is a tie, and 0 otherwise. The
conditional probability of $\kappa_{ab}^{h}$ given $I_{ab}=1$ is

$\bar{p}_{(\kappa_{ab}^{h})}=p_{(\kappa_{ab}^{h})}/\sum_{g\in
U}p_{(\kappa_{ab}^{g})}$

Provided $\bar{p}_{(\kappa_{ab}^{h})}>0$ for all $h\in U$, an unbiased
estimator of $Y$ based on $H_{m}$ can be given by

$\hat{Y}_{H}=\frac{1}{n_{m}} \sum_{h\in
H_{m}}y_{h}/\bar{p}_{(\kappa_{ab}^{h})}$

We set $w=0$ to reduce the chance of selecting the same node by LMHW; we can
allow for a small positive $r$ to ensure $\bar{p}_{(\kappa_{ab}^{h})}>0$ for
any $a<b$. We have then

$p_{(hh)h}=\frac{ru_{h}}{d_{h}+r}+\frac{1}{d_{h}+r}
\sum_{j\in\nu_{h}}\\{1-\min\Big{(}\frac{u_{j}}{u_{h}},1\Big{)}\\}=\frac{ru_{h}+d_{h}-A_{h}}{d_{h}+r}$

where $A_{h}=\sum_{j\in\nu_{h}}\min\big{(}u_{j}/u_{h},1\big{)}$, and for
$i\neq h$,

$p_{(ih)h}=\frac{ru_{h}}{d_{h}+r}+\begin{cases}\frac{d_{h}-A_{h}}{d_{h}+r}
&\text{if } i\not\in\nu_{h}\\\
\frac{d_{h}}{d_{h}+r}\big{\\{}1-\frac{1}{d_{h}-1}
\\{A_{h}-\min\big{(}\frac{u_{i}}{u_{h}},1\big{)}\\}\big{\\}}&\text{if }
i\in\nu_{h}\end{cases}$

and for $i\neq h$ and $j\neq h$,

$p_{(ih)j}=\frac{ru_{j}}{d_{h}+r}+\begin{cases}0&\text{if }
j\not\in\nu_{h}\text{ or }i=j\in\nu_{h}\\\
\frac{1}{d_{h}+r}\min\big{(}\frac{u_{j}}{u_{h}},1\big{)}&\text{if }
i\not\in\nu_{h},j\in\nu_{h}\\\
\frac{1}{d_{h}+r}\,\frac{d_{h}}{d_{h}-1}\min\big{(}\frac{u_{j}}{u_{h}},1\big{)}
&\text{if } i\neq j\in\nu_{h}\end{cases}$

## References

  * [1] Avrachenkov, K., Ribeiro, B., and Sreedharan, J. K. (2010). Improving Random Walk Estimation Accuracy with Uniform Restarts. _Research report_ , RR-7394, INRIA. inria-00520350
  * [2] Brin, S. and Page, L. (1998). The anatomy of a large-scale hypertextual web search engine. Computer Networks and ISDN Systems, 30:107-117.
  * [3] Grafström, A., Lundström, N.L.P. and Schelin, L. (2012) Spatially balanced sampling through the pivotal method. Biometrics, 66:514-520.
  * [4] Masuda, N., Porter, M.A. and Lambiotte, R. (2017) Random walks and diffusion on networks. _Physics Reports_ , 716-717: 1-58. http://dx.doi.org/10.1016/j.physrep.2017.07.007
  * [5] Stevens Jr., D.L. and Olsen, A.R. (2004). Spatially balanced sampling of natural resources. Journal of the American Statistical Association, 99:262-278.
  * [6] Thompson, S.K. (2006). Targeted random walk designs. Survey Methodology, 32, 11–24.
  * [7] Zhang, L.-C. (2022). _Graph sampling._ CRC Press.
  * [8] Zhang, L.-C. (2021). Graph sampling by lagged random walks. Stat, https://onlinelibrary.wiley.com/doi/abs/10.1002/sta4.444
  * [9] Zhang, L.-C. and Patone, M. (2017). Graph sampling. _Metron_ , 75:277-299.


#  Deuterated Polycyclic Aromatic Hydrocarbons in the Interstellar Medium: The
Aliphatic C–D Band Strengths
DRAFT:

X.J. Yang11affiliation: Department of Physics, Xiangtan University, 411105
Xiangtan, Hunan Province, China<EMAIL_ADDRESS>22affiliation: Department
of Physics and Astronomy, University of Missouri, Columbia, MO 65211, USA;
<EMAIL_ADDRESS>and Aigen Li22affiliation: Department of Physics and
Astronomy, University of Missouri, Columbia, MO 65211, USA<EMAIL_ADDRESS>
###### Abstract

Deuterium (D) was exclusively generated in the Big Bang and the standard Big
Bang Nucleosynthesis (BBN) model predicts a primordial abundance of D/H
$\approx$ 26 parts per million (ppm). As the Galaxy evolves, D/H gradually
decreases because of astration. The Galactic chemical evolution (GCE) model
predicts a present-day abundance of D/H $\gtrsim$ 20$\,{\rm ppm}$. However,
observations of the local interstellar medium (ISM) have revealed that the
gas-phase interstellar D/H varies considerably from one region to another and
has a median abundance of D/H $\approx$ 13$\,{\rm ppm}$, substantially lower
than predicted from the BBN and GCE models. It has been suggested that the
missing D atoms of D/H $\approx$ 7$\,{\rm ppm}$ could have been locked up in
deuterated polycyclic aromatic hydrocarbon (PAH) molecules. However, we have
previously demonstrated that PAHs with aromatic C–D units are insufficient to
account for the missing D. Here we explore if PAHs with aliphatic C–D units
could be a reservoir of D. We perform quantum chemical computations of the
vibrational spectra of “superdeuterated” PAHs (in which one D and one H share
an C atom) and PAHs attached with D-substituted methyl group, and derive the
band strengths of the aliphatic C–D stretch ($A_{4.65}$). By applying the
computationally derived $A_{4.65}$ to the observed aliphatic C–D emission at
$\sim\,$4.6–4.8$\,{\rm\mu m}$, we find that PAHs with aliphatic C–D units
could have tied up a substantial amount of D/H and marginally account for the
missing D. The possible routes to generate PAHs with aliphatic C–D units are
also discussed.

dust, extinction — ISM: lines and bands — ISM: molecules

##  1 Introduction

A set of infrared (IR) emission features at 3.3, 6.2, 7.7, 8.6, and
11.3$\,{\rm\mu m}$, commonly attributed to polycyclic aromatic hydrocarbon
(PAH) molecules, are ubiquitously seen in a wide variety of astrophysical
environments. We should emphasize that astronomical PAHs are not necessarily
pure aromatic compounds as strictly defined by chemists. Instead, PAH
molecules in space may include ring defects, substituents, partial
dehydrogenation, or superhydrogenation (see Yang et al. 2017 and references
therein). Astronomical PAHs could also be deuterated (e.g., see Allamandola et
al. 1989; Hudgins et al. 2004; Peeters et al. 2004; Draine 2006; Onaka et al.
2014, 2022; Mori et al. 2022; Yang et al. 2020, 2021). Indeed, it has been
suggested that PAHs could be a reservoir of deuterium (D) in the interstellar
medium (ISM; Draine 2006).

According to the standard Big Bang Nucleosynthesis (BBN) model, D was
exclusively created in the Big Bang with a primordial D/H abundance of
$\left[\rm D/H\right]_{\rm prim}$ $\approx$ 26$\,{\rm ppm}$ (see Yang et al.
2020 and references therein). Due to astration (i.e., D is converted to 3He,
4He, and heavier elements in stellar interiors), the present-day Galactic
interstellar D/H abundance should be reduced from the primordial abundance to
$\left[\rm D/H\right]_{\rm ISM}$ $\gtrsim$ $20\,{\rm ppm}$. However, the gas-
phase $\left[\rm D/H\right]_{\rm gas}$ abundance of the local ISM derived from
high resolution ultraviolet (UV) spectroscopic observations varies
considerably from one region to another, and ranges from $\sim\,$5$\,{\rm
ppm}$ to $\sim\,$22$\,{\rm ppm}$, with a median abundance of $\left[\rm
D/H\right]_{\rm gas}\approx 13\,{\rm ppm}$ (see Draine 2006). This implies
that there must exist a sink to accommodate the missing D/H, typically of
$\sim\,$7$\,{\rm ppm}$. The exact sink for the missing D atoms remains
unknown. If PAHs are responsible for the missing D atoms, they should
typically have tied up an amount of $\left[\rm
D/H\right]_{\scriptscriptstyle\rm PAH}\approx 7\,{\rm ppm}$. Here $\left[\rm
D/H\right]_{\scriptscriptstyle\rm PAH}$ refers to the amount of D locked up in
PAHs (relative to the total H in the ISM).

If interstellar PAHs are deuterated through replacing one or more of their
periphal hydrogen (H) atoms by D atoms, one expects to see an emission band at
$\sim\,$4.4$\,{\rm\mu m}$, arising from the aromatic C–D stretch, in analogy
to the 3.3$\,{\rm\mu m}$ C–H stretch. Indeed, the 4.4$\,{\rm\mu m}$ C–D
stretching emission band has been detected in photodissociated regions (PDRs),
reflection nebulae, and H ii regions both in the Milky Way and in the Large
and Small Magellanic Clouds (Peeters et al. 2004; Onaka et al. 2014, 2022;
Doney et al. 2016). The intrinsic band strengths of the 4.4$\,{\rm\mu m}$ C–D
stretch have been computed by Yang et al. (2020, 2021) and Buragohain et al.
(2015, 2016, 2020) and have been applied to astronomical observations to
estimate the degree of D enrichment in PAHs. Yang et al. (2020, 2021) have
found that small interstellar PAHs could indeed be D-enriched, exceeding the
interstellar D/H abundance by a factor of $\gtrsim$1000\. Nevertheless, this
is still not sufficient to account for the missing D atoms.

As will be discussed in §4.3, PAHs could also be deuterated by attaching
D-containing aliphatic sidegroups (e.g., –CH2D) to their edge carbon (C)
atoms, with the asymmetric and symmetric aliphatic C–D stretching modes
occuring at $\sim\,$4.6–4.8$\,{\rm\mu m}$ (see Tielens 1997). Indeed, these
C–D emission bands have also been seen in interstellar sources, together with
the 4.4$\,{\rm\mu m}$ aromatic C–D stretch (Peeters et al. 2004; Onaka et al.
2014, 2022; Doney et al. 2016). In regions rich in H and D atoms, PAHs could
be superhydrogenated by inserting an extra H atoms to an edge C atom so that
two H atoms share one C atom. This converts aromatic C–H stretchs to
aliphatic. In analogy, an D atom could also be squeezed in the periphal,
originally aromatic C–H unit to create an aliphatic C–H unit and an aliphatic
C–D unit.

The goal of this work is to evaluate the role of PAHs containing aliphatic C–D
units as a potential reservoir of D atoms. To this end, we first perform
quantum chemical computations of the vibrational spectra of PAHs with
deuterated aliphatic sidegroups and “superdeuterated” PAHs in which one H atom
and one D atom share one C atom. The computational methods and target PAH
molecules are summarized in §2. The computed spectra as well as the aliphatic
C–D band strengths derived from the computed spectra are presented in §3. In
§4.1 we compile all the astronomical data for the aliphatic C–D stretch. In
§4.2 we apply the aliphatic C–D band strengths derived here to astronomical
observations. We quantitatively determine the amount of D atoms possibly
locked up in PAHs with aliphatic C–D units. The possible routes for
aliphatically deuterating PAHs in the ISM are discussed in §4.3. Finally, we
summarize our major results in §5. This paper is largely concerned with the
quantum chemical computations of PAHs containing aliphatic C–D units. Readers
who are interested only in the aliphatic C–D band intensities and the
astrophysical applications may wish to proceed directly to §3.4.

##  2 Computational Methods and Target Molecules

We use the Gaussian16 software (Frisch et al. 2016) to calculate the IR
vibrational spectra for different parent PAH molecules containing aliphatic
C–D units. We consider seven parent aromatic molecules, i.e., benzene (C6H6),
naphthalene (C10H8), anthracene (C14H10), phenanthrene (C14H10), pyrene
(C16H10), perylene (C20H12), and coronene (C24H12). Here we focus on small
PAHs with fewer than 24 C atoms (i.e., $N_{\rm C}\lesssim 24$) since deuterium
enrichment in PAHs of $N_{\rm C}\gtrsim 40$ is not expected. Large PAHs have a
sufficiently large number of internal degrees of freedom to accommodate the
maximum energy of typical UV photons without subsequent photolytic bond
cleavage occurring (Hudgins et al. 2004). Nevertheless, if the deuteration of
PAHs occurs in D-enriched ice mantles through D-atom addition (see below),
large PAHs could also be deuterated.

We consider PAHs with two categories of aliphatic C–D units, i.e., PAHs
attached with one deuterated methyl group (hereafter PAH$\\_$CH2D; see Figure
1) and “superdeuterated” PAHs (hereafter PAH$\\_$HD; see Figure 2) in which
one edge C atom is shared by one H atom and one D atom. For both categories,
various isomers are considered. As H is much more abundant than D, if a PAH
molecule is “superdeuterated”, it is very likely that it will also be
superhydrogenated (i.e., two H atoms share one C atom). Therefore, we consider
“superdeuterated” PAHs to be superhydrogenated as well.

In dense regions where dust grains are coated by an ice mantle, upon
condensation onto the ice mantle and exposed to UV irradiation, PAHs could be
deuterated if interstellar ice is appreciably deuterated. Laboratory studies
have shown that UV photolysis of PAHs in D-enriched ices results in rapid D
enrichment of PAHs through aromatic D$\rightarrow$H exchange, D-atom addition
(“superdeuteration”), and exchange through keto-enol tautomerism (see Figure 4
of Sandford et al. 2000). The latter two reactions could also occur in the
absence of UV irradiation. For illustrative purpose, we therefore also
consider a deuterated molecule with a ketone (C=O) bond and a deuterated
molecule with a deuteroxyl (OD) sidegroup, using coronene as an example (see
Figure 2).

To examine the effects of ionization on the vibrational spectra and aliphatic
C–D band strengths of deuterated PAHs containing aliphatic C–D units, we also
consider the cationic counterparts of “superdeuterated” PAHs (i.e.,
PAH$\\_$HD). For simplicity, we only consider PAHs with a single H and D pair
(see Figure 3). They are essentially the radical counterparts of those
molecules shown in Figure 2, with one H atom stripped off from their H and H
pair.

In view of the fact that, observationally, the aliphatic C–D emission at
$\sim\,$4.6–4.8$\,{\rm\mu m}$ is always seen together with the aromatic C–D
emission at $\sim\,$4.4$\,{\rm\mu m}$ (Peeters et al. 2004; Onaka et al. 2014,
2022; Doney et al. 2016), we also consider PAHs with both aromatic and
aliphatic C–D bonds, taking pyrene as an example. We consider pyrenes in which
an D atom substitutes an H atom and a deuterated methyl group (–CH2D) is also
attached (hereafter PyreD$\\_$CH2D; see Figure 4). We also consider mono-
deuterated pyrenes which are “superdeuterated” as well (hereafter
PyreD$\\_$HD; see Figure 5). Again, as D is much less abundant than H,
“superdeuterated” pyrenes are very likely also superhydrogenated. Different
isomers are considered, with D atoms and functional groups attached at
different positions.

We will refer PAHs with deuterated methyl sidegroups (see Figure 1) by the
abbreviation of the first four letters of the names of their parental PAH
molecules, followed by the position where the methyl group is attached and
then “MD” (e.g., Naph1$\\_$MD refers to naphthalene with a deuterated-methyl
group attached at position 1). For “superdeuterated” PAHs with an H and D pair
as well as an H and H pair (see Figure 2), we use the abbreviation as
following: the first four letters of the names of their parental molecules
followed by “2H” and then the order of the isomer followed by “HD” (e.g.,
Naph2H$\\_$2HD refers to the second isomer of superhydrogenated,
superdeuterated naphthalene). For “superdeuterated” PAH cations (see Figure
3), we will also refer them by the abbreviation of the first four letters of
the names of their parental molecules, followed by the position where the D
atom is inserted and then “HD+” (e.g., Naph1$\\_$HD+ refers to naphthalene
cation with the H and D pair attached at position 1). For deuterated pyrenes
attached with CH2D (see Figure 4), we use “Pyre” followed by the position
where the methyl group is attached and then a letter of “D” followed by the
position of the aromatic C–D bonds (e.g., Pyre1$\\_$D3 refers to pyrene with
the deuterated-methyl group attached at position 1 and another D atom attached
at position 3). For deuterated pyrenes which are also “superdeuterated” (see
Figure 5), we use “Pyre2H” followed by the position where the H and D pair is
attached, followed by “HD” and then “D” which is further followed by the
position of the aromatic C–D bonds (e.g., Pyre2H$\\_$1HD$\\_$D3 refers to
superdeuterated pyrene with the H and D pair attached at position 1 and
another D atom attached at position 3).

We employ the hybrid density functional theoretical method (B3LYP) at the
6-311+G∗∗ level, which gives sufficient calculational accuracies with operable
computer time (see Yang et al. 2017 and reference therein). The standard
scaling is applied to the frequencies by employing a scale factor of
$\sim\,$0.9688 (Borowski 2012). We should note that the aromatic C–D band
strengths of deuterated PAHs derived by Yang et al. (2020, 2021) from
B3LYP/6-311+G∗∗ computations are consistent with recent laboratory
measurements of deuterated quenched carbonaceous composites (Mori et al.
2022).111The DFT-computed intensity ratio of the 4.4$\,{\rm\mu m}$ aromatic
C–D stretch to the 3.3$\,{\rm\mu m}$ aromatic C–H stretch is $\approx 0.56\pm
0.19$ for mono-deuterated PAHs (Yang et al. 2020) and $\approx 0.56\pm 0.03$
for multi-deuterated PAHs (Yang et al. 2021), while the experimentally
measured ratio is $\approx 0.56\pm 0.04$ (Mori et al. 2022).

##  3 Vibrational Spectra and Aliphatic C–D Band Strengths of Deuterated PAHs

###  3.1 PAHs Attached with a D-substituted Methyl Group

For PAHs attached with a D-substituted methyl group (i.e., PAH$\\_$CH2D), we
calculate the vibrational spectra of 18 isomers of seven parent PAH molecules
(see Figure 1), including all the possible positions where the methyl group is
attached. Since each of the three H atoms in the methyl group (–CH3) could be
substituted with an D atom, therefore in total we consider 54 isomers.

For each parent molecule, we obtain the mean spectrum averaged over all the
isomers and show in Figures 6a–g. For each vibrational transition, we assign a
width of 4$\,{\rm cm}^{-1}$. A feature at $\sim\,$2170$\,{\rm cm}^{-1}$ (i.e.,
$\sim\,$4.61$\,{\rm\mu m}$), originating from aliphatic C–D stretches, is
clearly seen in all these mean spectra. Also prominent are aliphatic and
aromatic C–H stretches at $\sim\,$2980$\,{\rm cm}^{-1}$ (i.e.,
$\sim\,$3.36$\,{\rm\mu m}$) and 3060 $\,{\rm cm}^{-1}$ (i.e.,
$\sim\,$3.30$\,{\rm\mu m}$), respectively. A close inspection of the aliphatic
C–D stretches shown in Figures 6a–g reveals two or three peaks, attributed to
different isomers. Nevertheless, these peaks are very close to each other,
indicating that the centroids of the C–D stretches do not vary much among
different isomers. As a matter of fact, the other vibrational modes are not
significantly blended either, indicating that the transition frequencies for
different isomers and different molecules are very similar, i.e, neither the
position of the methyl group attached nor the position of D atoms in the
methyl group affects the centroids of the vibrational features much.

Figure 6h shows the overall mean spectrum obtained by averaging over the mean
spectra of all seven parent molecules (normalized by $N_{\rm C}$). Again, the
aliphatic C–D and C–H stretches are clearly seen, while the aromatic C–H
stretch is the most prominent. The aliphatic C–D stretch also narrowly peaks,
implying that the centroids of the C–D stretches do not vary much among
different molecules. As the C–D in-plane and out-of-plane bending modes are
often well mixed with the C–H out-of-plane bending modes and the C–C–C
skeletal vibrational modes and therefore difficult to identify (Yang et al.
2020, 2021), we shall focus on the C–D and C–H stretching modes.

We tabulate the peak wavelengths and intensities of the C–H and C–D stretch
features in Table 1. The aliphatic C–D stretches generally occur at
$\sim\,$4.60$\,{\rm\mu m}$ and the deviation is negligible. In Figure 7 we
show the intensities of the aliphatic C–D stretches ($A_{4.65}$), the
aliphatic C–H stretches ($A_{3.4}$), and the aromatic C–H stretches
($A_{3.3}$) for all 54 isomers calculated here. The intensities of these
stretches, especially the aromatic C–H stretches, show small scatters among
different isomers. For aliphatic C–D stretches, the scatters mostly arise from
the different positions of D atoms in the methyl group attached to the parent
molecule. As D atoms substitute H atoms in the methyl group in a random
manner, we should average the intensities over that for the three positions of
the D atoms in the methyl group. This should also be the case for the methyl
group attached at different positions of a parent molecule. Nevertheless, the
scatters in the average intensities among the different isomers of the same
parent molecule are small (see Table 1). The overall mean band intensities, on
a per C–D or C–H bond basis, are $\langle A_{4.65}\rangle\approx 12.9\pm
0.89\,{\rm km}\,{\rm mol}^{-1}$, $\langle A_{3.4}\rangle\approx 23.1\pm
1.3\,{\rm km}\,{\rm mol}^{-1}$, and $\langle A_{3.3}\rangle\approx 14.3\pm
1.1\,{\rm km}\,{\rm mol}^{-1}$, respectively for the aliphatic C–D stretch,
the aliphatic C–H stretch, and the aromatic C–H stretch.

###  3.2 “Superdeuterated” PAHs

The superhydrogenation of PAHs often starts from the broken of a C=C bond,
followed by the addition of a pair of two H atoms in the neighbouring C atoms.
By “superdeuterated” PAHs (PAH$\\_$HD), we refer to superhydrogenated PAHs but
with one of the H atoms in the H and H pair substituted by an D atom, not two
D atoms share one C atom. As mentioned earlier, H is much more abundant than D
in space, if a PAH molecule is “superdeuterated”, it is very likely that it
will also be superhydrogenated. Therefore, for “superdeuterated” PAHs we
consider them to have a pair of two H atoms sharing an C atom, and meanwhile a
pair of H and D atoms sharing another C atom. In total, we consider 24 isomers
for seven parent molecules (see Figure 2).

For each parent molecule, we obtain the mean spectrum by averaging over the
isomers (see Figures 8a–g). We also obtain the overall mean spectrum by
averaging over all seven parent molecules (see Figure 8h). The aliphatic C–D
stretches are clearly seen at $\sim\,$2150$\,{\rm cm}^{-1}$ (i.e.,
$\sim\,$4.65$\,{\rm\mu m}$) and exhibit multiple peaks, indicating that the
wavelengths of the aliphatic C–D stretches somewhat differ from one isomer to
another.

The wavelengths and intensities of the C–H and C–D stretching transitions
calculated for superdeuterated PAHs are tabulated in Table 2 and shown in
Figure 9. Similar to PAHs attached with a D-substituted methyl group
(PAH$\\_$CH2D), the scatters in intensities are quite small among isomers and
parent molecules, especially for the aromatic C–H stretch. The aliphatic C–D
stretches of superdeuterated PAHs are more intense and peak at longer
wavelengths than that of PAH$\\_$CH2D, with $\langle\lambda_{\rm
peak}\rangle\approx 4.66\,{\rm\mu m}$ and $\langle A_{4.65}\rangle\approx
16.2\,{\rm km}\,{\rm mol}^{-1}$ for superdeuterated PAHs, while
$\langle\lambda_{\rm peak}\rangle\approx 4.60\,{\rm\mu m}$ and $\langle
A_{4.65}\rangle\approx 12.9\,{\rm km}\,{\rm mol}^{-1}$ for PAH$\\_$CH2D.

To investigate the effects of ionization on the vibrational spectra, we
calculate the spectra for “superdeuterated” PAH cations (PAH$\\_$HD+). We
consider the situation that only one C atom on the edge is superhydrogenated
and attached with a pair of H and D atoms. We consider seven parent molecules
and in total 18 isomers, with the H and D pair attached at different positions
(see Figure 3). For each parent molecule, the mean spectrum averaged over all
the isomers is given in Figures 10a–g. The overall mean spectrum averaged over
seven parent molecules is shown in Figure 10h. It is apparent that, upon
ionization, the C–H and C–D stretches are substantially depressed, while the
C–C strecthes at $\sim\,$6–9$\,{\rm\mu m}$ are appreciably enhanced. The
wavelengths and intensities of the C–H and C–D stretches are tabulated in
Table 3 and shown in Figure 11. With $\langle\lambda_{\rm peak}\rangle\approx
4.76\pm 0.05\,{\rm\mu m}$ and $\langle A_{4.65}\rangle\approx 1.64\pm
0.75\,{\rm km}\,{\rm mol}^{-1}$, the aliphatic C–D stretches of
“superdeuterated” PAH cations peak at significantly longer wavelengths and are
$\sim\,$10 times weaker than that of neutrals.

###  3.3 PAHs with Both Aromatic and Aliphatic C–D Bonds

Considering that the aliphatic C–D bands at $\sim\,$4.6–4.8$\,{\rm\mu m}$ are
always seen in emission together with the aromatic C–D bands at
$\sim\,$4.4$\,{\rm\mu m}$ (Peeters et al. 2004; Onaka et al. 2014, 2022; Doney
et al. 2016), we also calculate the vibrational spectra of PAHs with both
aliphatic and aromatic C–D bonds, with pyrene as an illustrative example. We
consider two kinds of aliphatic C–D units, with the aliphatic D atom either in
a methyl group (hereafter PyreD$\\_$CH2D, see Figure 4) or in a
“superdeuterated” pair (i.e., one D atom and one H atom share an C atom;
hereafter PyreD$\\_$HD, see Figure 5). We consider 23 isomers for
PyreD$\\_$CH2D, with the deuterated-mehtyl group attached at position 1, 2 or
4, and the other D atom bonding with one of the remaining edge C atoms in
pyrene. For PyreD$\\_$HD, we consider 24 isomers, with the H+H pair and the
H+D pair attached at positions 1 and 2, 2 and 3, or 4 and 5, while the other D
atom, again, bonding with one of the remaining edge C atoms in the rings.

By categorizing the isomers according to the positions of the aliphatic D
atom, we show in Figure 12 the spectra of the isomers in each category. We
also obtain the mean spectrum by averaging over all the somers and show in
Figure 14. The aliphatic and aromatic C–D stretches respectively at
2170$\,{\rm cm}^{-1}$ (i.e., 4.61$\,{\rm\mu m}$) and 2250$\,{\rm cm}^{-1}$
(i.e., 4.44$\,{\rm\mu m}$) are clearly seen in the spectra of isomers as well
as in the mean spectrum, with the aliphatic C–D stretch being appreciably more
intense than the aromatic C–D stretch. It is interesting to note that the
spectra of different isomers are quite similar (see Figures 12, 14), implying
that the positions of the D atoms do not significantly affect the spectra.

The wavelengths and intensities of the C–H and C–D stretches of PyreD$\\_$CH2D
and PyreD$\\_$HD are respectively tabulated in Tables 4, 5, and the
intensities are shown in Figures 13, 15. With $\langle\lambda_{\rm
peak}\rangle\approx 4.41\,{\rm\mu m}$ and $\langle A_{4.65}\rangle\approx
7.9\,{\rm km}\,{\rm mol}^{-1}$, the wavelengths and intensities of the
aromatic C–D stretch in PyreD$\\_$CH2D are essentially consistent with that of
PyreD$\\_$HD. This is not surprising because, according to the Born-
Oppenheimer approximation, the presence of an D-substituted methyl sidegroup
or a pair of H and D atoms sharing one C atom would not affect the aromatic
C–D stretches. On the other hand, the intensities of the aliphatic C–D
stretches of PyreD$\\_$CH2D ($\langle A_{4.65}\rangle\approx 12.4\,{\rm
km}\,{\rm mol}^{-1}$) differ appreciably from that of PyreD$\\_$HD ($\langle
A_{4.65}\rangle\approx 19.5\,{\rm km}\,{\rm mol}^{-1}$). This is
understandable since they arise from different types of aliphatic C–D units.

###  3.4 Aliphatic C–D Band Intensities

In Figures 16–20, we show $A_{4.65}/A_{3.3}$—the band-intensity ratio of the
4.65$\,{\rm\mu m}$ aliphatic C–D stretch to that of the 3.3$\,{\rm\mu m}$
aromatic C–H stretch—calculated for each set of our target molecules.
Generally speaking, $A_{4.65}/A_{3.3}$ does not vary much from one molecule to
another in each set. We obtain $\langle A_{4.65}/A_{3.3}\rangle\approx 0.90\pm
0.07$ for PAHs attached with a D-substituted methyl sidegroup (PAH$\\_$CH2D),
and $\langle A_{4.65}/A_{3.3}\rangle\approx 1.10\pm 0.32$ for
“superdeuterated” PAHs (PAH$\\_$HD) in which an extra D atom shares an C atom
with an H atom. For pyrenes attached with a D-substituted methyl sidegroup
(PyreD$\\_$CH2D), we obtain $\langle A_{4.65}/A_{3.3}\rangle\approx 0.85$ with
a standard deviation of $\sim\,$0.12; for “superdeuterated” pyrenes
(PyreD$\\_$HD), we obtain $\langle A_{4.65}/A_{3.3}\rangle\approx 1.31$ with a
standard deviation of $\sim\,$0.17.

In Table 6 we summarize the mean band strengths of the aliphatic C–D and
aromatic C–H stretches, obtained by averaging over all the species considered
here. The recommended band intensities are $\langle A_{4.65}\rangle\approx
15.3\pm 3.3\,{\rm km}\,{\rm mol}^{-1}$, $\langle A_{3.3}\rangle\approx 14.8\pm
0.33\,{\rm km}\,{\rm mol}^{-1}$, and $\langle A_{4.65}/A_{3.3}\rangle\approx
1.04\pm 0.21$ for deuterated neutral PAHs, and $\langle A_{4.65}\rangle\approx
1.64\pm 0.75\,{\rm km}\,{\rm mol}^{-1}$, $\langle A_{3.3}\rangle\approx
2.14^{+3.08}_{-2.14}\,{\rm km}\,{\rm mol}^{-1}$, and $\langle
A_{4.65}/A_{3.3}\rangle\approx 1.67\pm 1.00$ for deuterated PAH cations.222The
large scatter in $\langle A_{3.3}\rangle$ for PAH cations mainly arises from
the large $A_{3.3}$ intensity of Benz$\\_$HD+ (see Table 3), which is
substantially larger than that of the other molecules. If we exclude
Benz$\\_$HD+, we obtain $\langle A_{3.3}\rangle\approx 1.00\pm 0.59\,{\rm
km}\,{\rm mol}^{-1}$ and $\langle A_{4.65}/A_{3.3}\rangle\approx 1.92\pm
0.94$. However, this does not affect the degree of deuteration of PAHs derived
below in §4.2 since the C–D and C–H stretches are dominated by neutral PAHs
Since the intensities for the C–H or C–D stretches are significantly depressed
for PAH cations, when estimating the D enrichment from astronomical
observations, we recommend to adopt $A_{4.65}/A_{3.3}\approx 1.04\pm 0.21$,
the band-intensity ratio of deuterated neutral PAHs.

##  4 Astrophysical Applications

###  4.1 Aliphatic C–D Bands: Astronomical Observations

Observational efforts have been made with the Short Wavelength Spectrometer
(SWS) on board the Infrared Space Observatory (ISO) and the Infrared Camera
(IRC) on board AKARI to search for signals of deuterated PAHs through their
vibrational C–D bands (Verstraete et al. 1996; Peeters et al. 2004; Onaka et
al. 2014, 2022; Doney et al. 2016). So far, the detetcion of the aliphatic C–D
bands has been reported in eight sources (Peeters et al. 2004, Onaka et al.
2014, Doney et al. 2016). We compile all the available observational data for
the 3.3$\,{\rm\mu m}$ aromatic C–H stretch, the 3.4$\,{\rm\mu m}$ aliphatic
C–H stretch, the 4.4$\,{\rm\mu m}$ aromatic C–D stretch, and the
“4.65”$\,{\rm\mu m}$ aliphatic C–D stretch.333The aliphatic C–D stretch of
PAHs, consisting of asymmetric and symmetric modes, actually spans a range of
wavelengths, from $\sim\,$4.6$\,{\rm\mu m}$ to $\sim\,$4.8$\,{\rm\mu m}$. In
this work, by “4.65$\,{\rm\mu m}$” we actually mean the aliphatic C–D stretch
over the wavelength range of $\sim\,$4.6–4.8$\,{\rm\mu m}$.  The observational
emission intensities of these stretches ($I_{3.3}$, $I_{3.4}$, $I_{4.4}$, and
$I_{4.65}$) and their ratios with respect to the 3.3$\,{\rm\mu m}$ aromatic
C–H stretch, $\left(I_{4.65}/I_{3.3}\right)_{\rm obs}$,
$\left(I_{4.4}/I_{3.3}\right)_{\rm obs}$ and
$\left(I_{3.4}/I_{3.3}\right)_{\rm obs}$ are summarized in Table 7.

We explore whether these stretching bands are somewhat related to each other.
Figure 21a compares the observed emission intensities of the aliphatic C–D
stretch ($I_{4.65}$) with that of the aromatic C–D stretch ($I_{4.4}$). No
correlation is found between $I_{4.65}$ and $I_{4.4}$. To eliminate their
common correlation with the starlight radiation field (i.e., the observed
emission intensities of all these bands are proportional to the starlight
intensity), we normalize $I_{4.65}$ and $I_{4.4}$ by $I_{3.3}$ and then
compare $I_{4.65}/I_{3.3}$ with $I_{4.4}/I_{3.3}$. As demonstrated in Figure
21b, $I_{4.65}/I_{3.3}$ and $I_{4.4}/I_{3.3}$ are not correlated. This
indicates that the deuteration of the aliphatic C atoms may not occur together
with the deuteration of the aromatic C atoms (i.e., the aliphatic and aromatic
C atoms get deuterated independently). Also, to examine the degree of
deuteration of the aliphatic C atoms and that of the aromatic C atoms, we
compare $I_{4.65}/I_{3.4}$ with $I_{4.4}/I_{3.3}$ and find no correlation
either (see Figure 21c). This also implies that the deuterations of the
aliphatic and aromatic C atoms occur independently. Finally, we compare
$I_{4.65}/I_{4.4}$ with $I_{3.4}/I_{3.3}$. This is to investigate if the
aliphaticity of the C–D bond is related to that of the C–H bond. As shown in
Figure 21d, they do not correlate.

###  4.2 Deuterium Depletion onto PAHs

We now evaluate from the observed emission intensities of the 4.65$\,{\rm\mu
m}$ aliphatic C–D stretch the number of D atoms (relative to H) that could
have been locked up in the aliphatic C–D units of PAHs. For a PAH molecule
containing $N_{\rm H}$ aromatic H atoms, and ${N_{\rm D}}$ aliphatic D atoms,
we define the degree of deuteration of the molecule as $f_{\rm
aliCD}\equiv{N_{\rm D}}/\left(N_{\rm H}+{N_{\rm D}}\right)$, the fraction of
the total (H and D) atoms in the form of D. The observed ratio of the power
emitted from the 4.65$\,{\rm\mu m}$ aliphatic C–D band to that from the
3.3$\,{\rm\mu m}$ aromatic C–H band is related to the degree of deuteration
through

$\left(\frac{I_{4.65}}{I_{3.3}}\right)_{\rm
obs}\approx\left(\frac{A_{4.65}}{A_{3.3}}\right)\times\left(\frac{{N_{\rm
C-D}}}{N_{\rm C-H}}\right)\times\left(\frac{B_{4.65}}{B_{3.3}}\right)~{}~{},$
(1)

where $N_{\rm C-H}$ and ${N_{\rm C-D}}$ are respectively the number of C–H and
C–D bonds in a deuterated PAH molecule, $B_{\lambda}(T)$ is the Planck
function at wavelength $\lambda$ and temperature $T$, and $A_{3.3}$ and
$A_{4.65}$ are the intrinsic band strengths of the aromatic C–H and aliphatic
C–D stretches (on a per C–H or C–D bond basis). Since the observed C–H and C–D
stretching emission bands arise predominantly from neutral PAHs, we will adopt
$A_{4.65}/A_{3.3}\approx 1.04\pm 0.21$ calculated for neutral PAHs containing
aliphatic C–D units (see §3.4 and Table 6).

Considering the relatively small abundance of D atoms in space, it is unlikely
that an aliphatic C atom would be doubly deuterated (i.e., two D atoms share
one C atom). Therefore, it is reasonable to assume ${N_{\rm C-D}}={N_{\rm D}}$
and $N_{\rm C-H}=N_{\rm H}$. The degree of deuteration of the molecule becomes

$f_{\rm aliCD}\approx\left\\{1\,+\,\left(\frac{I_{3.3}}{I_{4.65}}\right)_{\rm
obs}\times\left(\frac{A_{4.65}}{A_{3.3}}\right)\times\left(\frac{B_{4.65}}{B_{3.3}}\right)\right\\}^{-1}~{}~{}.$
(2)

The 3.3$\,{\rm\mu m}$ C–H stretch and 4.65$\,{\rm\mu m}$ C–D stretch are most
effectively emitted respectively by PAHs of vibrational temperatures of
$\sim\,$730$\,{\rm K}$ and $\sim\,$520$\,{\rm K}$ (see Foonote No. 10 in Yang
et al. 2020), stochastically heated by individual UV photons (Draine & Li
2001). For $520\lesssim T\lesssim 730\,{\rm K}$, we obtain
$B_{3.3}/B_{4.65}\approx 0.72\pm 0.14$. By relaxing the temperature range to
$400\lesssim T\lesssim 900\,{\rm K}$, we obtain $B_{3.3}/B_{4.65}\approx
0.78\pm 0.32$. Even if we extend the temperature range to $300\lesssim
T\lesssim 1000\,{\rm K}$, a similar ratio of $B_{3.3}/B_{4.65}\approx 0.78\pm
0.44$ is still obtained. In the following, we will adopt
$B_{3.3}/B_{4.65}\approx 0.78\pm 0.32$ that is suitable for $400\lesssim
T\lesssim 900\,{\rm K}$.

With $B_{3.3}/B_{4.65}\approx 0.78\pm 0.32$ and $A_{4.65}/A_{3.3}\approx
1.04\pm 0.21$ (see §3.4), for each of the eight sources we obtain $f_{\rm
aliCD}$ from $\left(I_{4.65}/I_{3.3}\right)_{\rm obs}$. These sources also
exhibit an aromatic C–D stretch at $\sim\,$4.4$\,{\rm\mu m}$ as well as an
aliphatic C–H stretch at 3.4$\,{\rm\mu m}$. Following Yang et al. (2020), for
each source we derive $f_{\rm aroCD}$, the deuteration degree of PAHs based on
the aromatic C–D stretch. Similarly, following Yang et al. (2013), for each
source we also derive $f_{\rm aliCH}$, the aliphatic fraction of PAHs based on
the aliphatic C–H stretch. The results are tabulated in Table 7.

It is gratifying that in some sources (e.g., NGC 3603, IRAS 12073, M17b)
$f_{\rm aliCD}$ can be as much as $\sim\,$20%, implying that PAHs could be a
plausible sink of the missing D atoms. In the ISM, PAHs require a total C/H
abundance of $\left[\rm C/H\right]_{\scriptscriptstyle\rm PAH}\approx 60\,{\rm
ppm}$ (Li & Draine 2001), where $\left[\rm C/H\right]_{\scriptscriptstyle\rm
PAH}$ refers to the amount of C locked up in PAHs (relative to the total H in
the ISM). For small PAHs like coronene, the ratio of the number of C atoms to
the number of H atoms is $N_{\rm C}/N_{\rm H}\approx 2$. Therefore, PAHs may
contain an amount of $\sim\,$30$\,{\rm ppm}$ of H (relative to the total H in
the ISM). A combination of $f_{\rm aliCD}$ and $f_{\rm aroCD}$ indicates that
in some sources up to $\sim\,$28% of the H atoms in PAHs could be deuterated
(see Table 7), i.e., $\left[\rm D/H\right]_{\scriptscriptstyle\rm PAH}\approx
8.4\,{\rm ppm}$ – the amount of D locked up in PAHs (relative to the total H
in the ISM). This exceeds the median of the missing D/H ($\sim\,$7$\,{\rm
ppm}$; see §1). We note that aliphatic C–D is often created from the
deuteration of methyl-group or superdeuteration. This naturally accompanies
with at least one aliphatic C–H bond. Therefore, one would expect a higher
$f_{\rm aliCH}$ than $f_{\rm aliCD}$. Indeed, Table 7 shows that in some
sources $f_{\rm aliCH}$ is twice as much as $f_{\rm aliCD}$.

Finally, we should also note that the above estimation of $\approx$ 8.4$\,{\rm
ppm}$ of D depleted in PAHs (relative to the total H in the ISM) is rather
generous. This is based on the assumption that both small ($N_{\rm C}\lesssim
40$) and large PAHs ($N_{\rm C}\gtrsim 100$) are deuterated. For large PAHs,
deuteration may occur in D-enriched ice mantles (see §4.3). Even if they are
deuterated, due to their large heat capacities, large PAHs are not excited
sufficiently to temperatures high enough to emit at the C–H and C–D stretching
bands.

In addition to large PAHs, hydrogenated amorphous carbon (HAC) solids could
also deplete some of the D atoms missing from the gas phase. Furton et al.
(1999) found that the interstellar 3.4$\,{\rm\mu m}$ aliphatic C–H absorption
band is most consistent with $\sim\,$80$\,{\rm ppm}$ of C (relative to the
total H in the ISM) with an H/C ratio of $\sim\,$0.5. Compared with PAHs, HAC
would provide a larger number of potential sites for D-depletion. It would be
interesting to see if the lines of sight exhibiting strong 3.4$\,{\rm\mu m}$
aliphatic C–H absorption also show detectable aliphatic C–D absorption at
$\sim\,$4.65$\,{\rm\mu m}$.

In the local diffuse ISM, the optical depth of the 3.4$\,{\rm\mu m}$ aliphatic
C–H absorption relative to the visual extinction ($A_{V}$) is
$\Delta\tau_{3.4\,{\rm\mu m}}/A_{V}\approx 1/274$ (Sandford et al. 1991,
Pendleton et al. 1994; see Gao et al. 2010 for a summary). Toward the Galactic
center, with $\Delta\tau_{3.4\,{\rm\mu m}}/A_{V}\approx 1/146$ (McFadzean et
al. 1989, Pendleton et al. 1994; also see Gao et al. 2010 for a summary), the
3.4$\,{\rm\mu m}$ absorption band is stronger (on a per magnitude of visual
extinction basis) by a factor of $\sim\,$1.9. If we take the band-intensity
ratio of the 4.65$\,{\rm\mu m}$ aliphatic C–D to the 3.4$\,{\rm\mu m}$
aliphatic C–H to be $A_{4.65}/A_{3.4}\approx 0.38$ as measured for deuterated
quenched carbonaceous composite materials (Mori et al. 2022), and assume that
HAC contains $\sim\,$40$\,{\rm ppm}$ of H and $\sim\,$7$\,{\rm ppm}$ of D
(relative to the total H in the ISM), we expect the optical depth of the
4.65$\,{\rm\mu m}$ aliphatic C–D absorption relative to the visual extinction
to be $\Delta\tau_{4.65\,{\rm\mu m}}/A_{V}\approx 1/4150$ for the local ISM
and $\Delta\tau_{4.65\,{\rm\mu m}}/A_{V}\approx 1/2200$ for the Galactic
center. Such a weak band is beyond the detection capability of ISO. JWST’s
unprecedented sensitivity will provide a unique opportunity to detect or place
an upper limit on the 4.65$\,{\rm\mu m}$ aliphatic C–D absorption band in
heavily obscured lines of sight.

###  4.3 Possible Routes for Aliphatically Deuterating PAHs

A significant fraction of the UV radiation emitted by massive stars impinges
on the molecular gas associated with star formation. This results in a PDR
bounded on one side by an ionization front, and on the other side by cold
molecular gas which has not yet been appreciably affected by UV radiation.

As schematically illustrated in Figure 22, it is natural to expect that
interstellar PAHs are deuterated in such regions. Lets take the Orion Bar as
an example. The regions close to the ionization front are rich in UV photons
and deuteration could occur during photodissociation events. Upon absorption
of an energetic photon, one or more periphal H atoms could be stripped off and
replaced by D atoms. Due to the zero-point energy difference, C–D bonds are
more stable than C–H bonds and therefore PAHs of intermediate size are
expected to become deuterium enriched and obtain aromatic C–D units. In these
regions one expects to see the 4.4$\,{\rm\mu m}$ emission band arising from
aromatic C–D stretch.

When moving away from the ionization front, the regions intermediate beteen
the ionization front and the dissociation front, are rich in H and D atoms.
PAHs could be superhydrogenated (i.e., two H atoms share one C atom) and even
superdeuterated (i.e., one D atom and one H atom share one C atom). In these
regions one expects to see the aliphatic C–D stretching emission at
$\sim\,$4.6–4.8$\,{\rm\mu m}$ as well as the aliphatic C–H stretching emission
at 3.4$\,{\rm\mu m}$ (see Figure 22 for illustration).

Another way to deuterate PAHs and generate aliphatic C–D bonds is through the
addition of an D-substituted methyl group. This could occur in dense clouds
where PAHs and simple molecules (e.g., CH4) are frozen into the ice mantles
coated on dust grains. In D-enriched ices, exposed to UV irradiation, PAHs
could attain an D-substituted methyl group through chemical reactions with
deuterated methane (see Figure 23 for a schematic illustration). PAHs could
also become deuterated in D-enriched ices by exchanging H atoms with D atoms.
This H/D exchange would be enhanced for oxidized PAHs. Furthermore, UV-
stimulated addition of D atoms at aromatic edge sites could yield aliphatic
C–D bonds. Sandford et al. (2000) have experimentally demonstrated that, when
exposed to UV radiation, PAHs in D-enriched ices could be oxidized and easily
deuterated.

As deuterated PAHs made in UV-irradiated ices are likely also oxidized, we
also compute the vibrational spectra of deuterated, oxidized PAHs, taking
coronene as an example. We consider two coronene derivatives, with D-atom
exchange at oxidized edge sites or at aromatic edge sites. The structures of
these two molecules are also shown in Figure 2. We show in Figure 24 their
vibrational spectra, together with two “superdeuterated” species. The
aliphatic and aromatic C–D stretches as well as the O–D stretch are clearly
seen. We note that in Coro$\\_$O2D (in which two D atoms are attached to one C
atom) the aliphatic C–D stretch is significantly depressed.

##  5 Summary

We have explored whether interstellar PAHs could be a major reservoir of D to
account for its depletion from the gas phase, focusing on the aliphatic C–D
units of PAHs. The major results are as follows:

  1. 1.

We have employed the hybrid DFT method B3LYP in conjunction with the 6-311+G∗∗
basis set to calculate the vibrational spectra of deuterated PAHs with
aliphatic C–D units, including “superdeuterated” PAHs in which one D and one H
share an C atom and PAHs attached with an D-substituted methyl group.

  2. 2.

For all these aliphatically deuterated species and their isomers, we have
derived from the computed spectra the intrinsic band strengths of the
4.65$\,{\rm\mu m}$ aliphatic C–D stretch ($A_{4.65}$), and the 3.3$\,{\rm\mu
m}$ aromatic C–H stretch ($A_{3.3}$). Both the C–H and C–D stretches
predominantly arise from neutral PAHs. By averaging over all these molecules,
we have determined the mean band strengths to be $\langle
A_{4.65}\rangle\approx 15.3\pm 3.3\,{\rm km}\,{\rm mol}^{-1}$, $\langle
A_{3.3}\rangle\approx 14.8\pm 0.33\,{\rm km}\,{\rm mol}^{-1}$, and $\langle
A_{4.65}/A_{3.3}\rangle\approx 1.04\pm 0.21$.

  3. 3.

We have compiled all the observational data from ISO/SWS and AKARI/IRC and
derived $\left(I_{4.4}/I_{3.3}\right)_{\rm obs}$, the ratio of the power
emitted from the aliphatic C–D band at 4.6–4.8$\,{\rm\mu m}$ ($I_{4.65}$) to
that from the aromatic C–H band at 3.3$\,{\rm\mu m}$ ($I_{3.3}$). By comparing
the computationally-derived mean ratio of $\langle
A_{4.65}/A_{3.3}\rangle\approx 1.04$ for aliphatically deuterated PAHs, it is
found the degree of deuteration of some sources could be as high as
$\sim\,$20%. Combined with PAHs containing aromatic C–D units, the amount of
D/H tied up in deuterated PAHs could exceed that missing from the gas phase,
therefore consistent with the predictions from the BBN and GCE models and the
gas-phase D/H of the local ISM.

  4. 4.

It is not clear how the high deuteration of PAHs required to account for the
missing D is achieved. We have suggested that the deuteration of PAHs with
aromatic C–D units may mainly occur in UV-rich regions, while the deuteration
of PAHs with aliphatic C–D units may mainly occur in benign regions rich in H
and D atoms or in dense clouds.

We thank B.T. Draine, A.N. Witt and the anonymous referee for valuable
suggestions. XJY is supported in part by NSFC 12122302 and 11873041. AL is
supported in part by NASA grants 80NSSC19K0572 and 80NSSC19K0701.

## References

  * (1) Allamandola, L.J., Tielens, A.G.G.M., & Barker, J.R. 1985, ApJ, 290, L25
  * (2) Buragohain, M., Pathak, A., Sarre, P., Onaka, T., & Sakon, I. 2015, MNRAS, 454, 193
  * (3) Buragohain, M., Pathak, A., Sarre, P., Onaka, T., & Sakon, I. 2016, Planet. Space Sci., 133, 97
  * (4) Buragohain, M., Pathak, A., Sakon, I., & Onaka, T. 2020, ApJ, 892, 11
  * (5) Doney, K. D., Candian, A., Mori, T., Onaka, T., & Tielens, A. G. G. M. 2016, A&A, 586, 65
  * (6) Draine, B. T. 2006, in ASP Conf. Ser. 348, Astrophysics in the Far Ultraviolet: Five Years of Discovery with FUSE, ed. G. Sonneborn, H. Moos, & B. -G. Andersson (San Francisco, CA: ASP), 58
  * (7) Draine, B.T., & Li, A. 2001, ApJ, 551, 807
  * (8) Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2016, Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford CT
  * (9) Furton, D. G., Laiho, J. W., & Witt, A. N. 1999, ApJ, 526, 752
  * (10) Gao, J., Jiang, B.W., & Li, A. 2010, EPS, 62, 63
  * (11) Hudgins, D. M., Bauschlicher, C. W., Jr., & Sandford, S. A. 2004, ApJ, 614, 770
  * (12) Li, A., & Draine, B.T. 2001, ApJ, 554, 778
  * (13) McFadzean, A. D., Whittet, D.C.B., Bode, M.F., Adamson, A.J., & Longmore, A.J. 1989, MNRAS, 241, 873
  * (14) Mori, T., Onaka, T., Sakon, I., et al. 2022 ApJ, 933, 35
  * (15) Onaka, T., Mori, T. I., Sakon, I., et al. 2014, ApJ, 780, 114
  * (16) Onaka, T., Sakon, I., & Shimonishi, T., 2022, ApJ, 941, 190
  * (17) Pendleton, Y.J., Sandford, S.A., Allamandola, L.J., Tielens, A.G.G.M., & Sellgren, K. 1994, ApJ, 437, 683
  * (18) Peeters, E., Allamandola, L. J., Bauschlicher, C. W., Jr., et al. 2004, ApJ, 604, 252
  * (19) Sandford, S.A., Allamandola, L.J., Tielens, A.G.G.M., Sellgren, K., Tapia, M., & Pendleton, Y.J. 1991, ApJ, 371, 607
  * (20) Sandford, S. A., Bernstein, M. P., Allamandola, L. J., Gillette, J. S., & Zare, R. N. 2000, ApJ, 538, 691
  * (21) Tielens, A. G. G. M. 1997, Astrophysical Implications of the Laboratory Study of Presolar Materials, ed T. H. Bernatowicz & E. K. Zinner (New York: AIP) 523
  * (22) Verstraete, L., Puget, J. L., Falgarone, E., et al. 1996, A&A, 315, L337
  * (23) Yang, X. J., Glaser, R., Li, A., & Zhong, J. X. 2013, ApJ, 776, 110
  * (24) Yang, X. J., Glaser, R., Li, A., & Zhong, J. X. 2017, New Astron. Rev., 77, 1
  * (25) Yang, X. J., Li, A., & Glaser, R. 2020, ApJS, 251, 12
  * (26) Yang, X. J., Li, A., He, C. Y., & Glaser, R. 2021, ApJS, 255, 23

Figure 1:  Structures of methyl-deuterated PAHs. We refer to a methyl-
deuterated species by the abbreviation of the first four letters of the name
of its parent molecule followed by the position where the methyl group is
attached (e.g., Anth2 refers to anthracene with the deuterated methyl attached
at position 2). Since there are three H atoms in methyl and any of them could
be deuterated, therefore all three possibilities are taken into account.
Figure 2:  Structures of “superdeuterated” PAHs in which one D and one H atom
share an C atom. Since H is much more abundant than D in space, if a PAH
molecule is “superdeuterated”, it is very likely that it will also be
superhydrogenated. Also shown are doubly “superdeuterated” coronene (with two
pairs of H and D atoms), and deuterated, oxidized coronenes (see §4.3).
Figure 3:  Structures of “superdeuterated” PAH cations.  Figure 4:  Structures
of 23 isomers of deuterated pyrenes containing a methyl-deuterated sidegroup
as well as a periphal D atom. These molecules have both aliphatic and aromatic
C–D bonds.  Figure 5:  Structures of 24 isomers of “superdeuterated” pyrenes
containing an H+H pair, an H+D pair, as well as a periphal D atom. Similar to
that in Figure 4, these molecules also have both aliphatic and aromatic C–D
bonds.  Figure 6:  DFT-computed vibrational spectra of methyl-deuterated PAHs
(see Figure 1), obtained by averaging over all the isomers of each parent
molecule. Also shown is the overall mean spectrum obtained by averaging over
the mean spectra of all seven parent molecules (normalized by $N_{\rm C}$).
The aliphatic C–D, aliphatic C–H, and aromatic C–H stretches are labelled
respectively by “aliCD”, “aliCH” and “aroCH”.  Figure 7:  Band strengths of
the 4.65$\,{\rm\mu m}$ aliphatic C–D stretches ($A_{4.65}$; red), the
3.4$\,{\rm\mu m}$ aliphatic C–H stretches ($A_{3.4}$; green), and the
3.3$\,{\rm\mu m}$ aromatic C–H stretches ($A_{3.3}$; blue) computed at level
B3LYP/6-311+G∗∗ for methyl-deuterated PAHs.  Figure 8:  Same as Figure 6 but
for “superdeuterated” PAHs in which one H atom and one D atom share an C atom
(see Figure 2).  Figure 9:  Same as Figure 7 but for “superdeuterated” PAHs.
Figure 10:  Same as Figure 6 but for “superdeuterated” PAH cations (see Figure
3).  Figure 11:  Same as Figure 7 but for “superdeuterated” PAH cations.

Figure 12:  Same as Figure 6 but for pyrenes containing a methyl-deuterated
sidegroup as well as a periphal D atom (see Figure 4).  Figure 13:  Same as
Figure 7 but for pyrenes containing a methyl-deuterated sidegroup as well as a
periphal D atom.  Figure 14:  Same as Figure 6 but for “superdeuterated”
pyrenes containing an H+H pair, an H+D pair, as well as a periphal D atom (see
Figure 5).  Figure 15:  Same as Figure 7 but for “superdeuterated” pyrenes
containing an H+H pair, an H+D pair, as well as a periphal D atom.  Figure 16:
Band-strength ratios of $A_{4.65}/A_{3.3}$ (magenta) and $A_{4.65}/A_{3.4}$
(green) for methyl-deuterated PAHs (see Figure 1).  Figure 17:  Band-strength
ratios of $A_{4.65}/A_{3.3}$ (magenta) and $A_{4.65}/A_{3.4}$ (green) for
“superdeuterated” PAHs in which one H atom and one D atom share an C atom (see
Figure 2).  Figure 18:  Band-strength ratios of $A_{4.65}/A_{3.3}$ (magenta)
and $A_{4.65}/A_{3.4}$ (green) for “superdeuterated” PAH cations (see Figure
3).  Figure 19:  Band-strength ratios of $A_{4.65}/A_{3.3}$ (magenta)
$A_{4.65}/A_{3.4}$ (green), and $A_{4.4}/A_{3.3}$ (blue) for pyrenes
containing a methyl-deuterated sidegroup as well as a periphal D atom (see
Figure 4).  Figure 20:  Band-strength ratios of $A_{4.65}/A_{3.3}$ (magenta)
$A_{4.65}/A_{3.4}$ (green), and $A_{4.4}/A_{3.3}$ (blue) for “superdeuterated”
pyrenes containing an H+H pair, an H+D pair, as well as a periphal D atom (see
Figure 5).  Figure 21:  Upper left (a): Correlation of the observed
intensities of the 4.65$\,{\rm\mu m}$ aliphatic C–D band ($I_{4.65}$) with
that of the 4.4$\,{\rm\mu m}$ aromatic C–D band ($I_{4.4}$). Upper right (b):
Correlation of $I_{4.65}/I_{3.3}$ with $I_{4.4}/I_{3.3}$. Bottom left (c):
Correlation of $I_{3.4}/I_{3.3}$ with $I_{4.4}/I_{3.3}$. Bottom right (d):
Correlation of $I_{3.4}/I_{3.3}$ with $I_{4.65}/I_{4.4}$.  Figure 22:
Schematic illustration of the deuteration of PAHs in the Orion Bar PDR, from
regions near the ionization front where the C–D bonds are mostly aromatic, to
regions intermediate beteen the ionization front and the dissociation front
where PAHs could be superhydrogenated and even superdeuterated.  Figure 23:
Schematic illustration of the deuteration of PAHs in the M17 dense molecular
cloud. The deuteration occurs in the UV irradiated, D-enriched ice mantles
coated on dust grains and generates aliphatic C–D units.  Figure 24:
Vibrational spectra of “superdeuterated” coronene (a), doubly
“superdeuterated” coronene (b), and deuterated, oxidized coronenes (c, d).

Table 1: Wavelengths ($\lambda$) and Intensities ($A_{\lambda}$ in $\,{\rm km}\,{\rm mol}^{-1}$) of the Nominal “3.3$\,{\rm\mu m}$” Aromatic C–H Stretch, “3.4$\,{\rm\mu m}$” Aliphatic C–H Stretch and “4.65$\,{\rm\mu m}$” Aliphatic C–D Stretch Computed at the B3LYP/6-311+G∗∗ Level for All the Methyl-Deuterated PAHs Shown in Figure 1.  Compound | $\lambda_{3.3}$ | $A_{3.3}$ | $\lambda_{3.4}$ | $A_{3.4}$ | $\lambda_{4.65}$ | $A_{4.65}$ | $A_{4.65}/A_{3.3}$
---|---|---|---|---|---|---|---
|  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |
Benz$\\_$CH2D | 3.25 | 14.42 | 3.37 | 21.29 | 4.60 | 11.53 | 0.80
Naph$\\_$CH2D | 3.25 | 13.88 | 3.37 | 22.58 | 4.60 | 12.37 | 0.89
Anth$\\_$CH2D | 3.25 | 13.82 | 3.37 | 22.24 | 4.59 | 12.28 | 0.89
Phen$\\_$CH2D | 3.24 | 13.07 | 3.37 | 22.64 | 4.59 | 12.43 | 0.95
Pyre$\\_$CH2D | 3.25 | 14.66 | 3.38 | 23.93 | 4.60 | 13.62 | 0.93
Pery$\\_$CH2D | 3.25 | 13.73 | 3.38 | 24.38 | 4.60 | 13.73 | 1.00
Coro$\\_$CH2D | 3.25 | 16.58 | 3.38 | 24.83 | 4.60 | 14.24 | 0.86
Average | 3.25 | 14.31 | 3.37 | 23.13 | 4.60 | 12.89 | 0.90
stdev | 0.00 | 1.12 | 0.00 | 1.28 | 0.00 | 0.98 | 0.07
Table 2: Wavelengths ($\lambda$) and Intensities ($A_{\lambda}$ in $\,{\rm km}\,{\rm mol}^{-1}$) of the Nominal “3.3$\,{\rm\mu m}$” Aromatic C–H Stretch, “3.4$\,{\rm\mu m}$” Aliphatic C–H Stretch and “4.65$\,{\rm\mu m}$” Aliphatic C–D Stretch Computed at the B3LYP/6-311+G∗∗ Level for All the “Supedeuterated” PAHs Shown in Figure 2.  Compound | $\lambda_{3.3}$ | $A_{3.3}$ | $\lambda_{3.4}$ | $A_{3.4}$ | $\lambda_{4.65}$ | $A_{4.65}$ | $A_{4.65}/A_{3.3}$
---|---|---|---|---|---|---|---
|  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |
Benz$\\_$HD | 3.26 | 16.47 | 3.43 | 32.84 | 4.67 | 17.82 | 1.08
Naph$\\_$HD | 3.25 | 14.69 | 3.42 | 32.39 | 4.70 | 19.32 | 1.31
Anth$\\_$HD | 3.25 | 14.47 | 3.43 | 28.88 | 4.66 | 15.98 | 1.10
Phen$\\_$HD | 3.25 | 13.59 | 3.42 | 32.29 | 4.66 | 17.82 | 1.31
Pyre$\\_$HD | 3.25 | 15.07 | 3.42 | 32.18 | 4.62 | 17.69 | 1.17
Pery$\\_$HD | 3.24 | 13.71 | 3.43 | 33.05 | 4.65 | 17.82 | 1.30
Coro$\\_$HD | 3.27 | 17.24 | 3.40 | 37.89 | 4.71 | 6.97 | 0.40
Average | 3.25 | 15.03 | 3.42 | 32.79 | 4.66 | 16.20 | 1.10
Stdev | 0.01 | 1.36 | 0.01 | 2.65 | 0.03 | 4.18 | 0.32
Table 3: Wavelengths ($\lambda$) and Intensities ($A_{\lambda}$ in $\,{\rm km}\,{\rm mol}^{-1}$) of the Nominal “3.3$\,{\rm\mu m}$” Aromatic C–H Stretch, “3.4$\,{\rm\mu m}$” Aliphatic C–H Stretch “4.65$\,{\rm\mu m}$” Aliphatic C–D Stretch Computed at the B3LYP/6-311+G∗∗ Level for All the “Superdeuterated” PAH Cations Shown in Figure 3.  Compound | $\lambda_{3.3}$ | $A_{3.3}$ | $\lambda_{3.4}$ | $A_{3.4}$ | $\lambda_{4.65}$ | $A_{4.65}$ | $A_{4.65}/A_{3.3}$
---|---|---|---|---|---|---|---
|  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |
Benz$\\_$HD+ | 3.24 | 9.01 | 3.31 | 0.13 | 4.81 | 1.61 | 0.18
Naph$\\_$HD+ | 3.22 | 0.84 | 3.51 | 0.22 | 4.66 | 0.64 | 0.76
Anth$\\_$HD+ | 3.24 | 1.32 | 3.51 | 2.28 | 4.79 | 2.23 | 1.68
Phen$\\_$HD+ | 3.23 | 0.33 | 3.49 | 2.15 | 4.75 | 1.18 | 3.56
Pyre$\\_$HD+ | 3.24 | 0.46 | 3.50 | 1.39 | 4.78 | 0.95 | 2.09
Pery$\\_$HD+ | 3.22 | 1.10 | 3.51 | 2.74 | 4.78 | 2.27 | 2.07
Coro$\\_$HD+ | 3.23 | 1.92 | 3.48 | 5.20 | 4.75 | 2.62 | 1.36
Average | 3.23 | 2.14 | 3.47 | 2.02 | 4.76 | 1.64 | 1.67
Stdev | 0.01 | 3.08 | 0.07 | 1.73 | 0.05 | 0.75 | 1.00
Table 4: Wavelengths ($\lambda$) and Intensities ($A_{\lambda}$ in $\,{\rm km}\,{\rm mol}^{-1}$) of the Nominal “3.3$\,{\rm\mu m}$” Aromatic C–H Stretch, “3.4$\,{\rm\mu m}$” Aliphatic C–H Stretch, “4.4$\,{\rm\mu m}$” Aromatic C–D Stretch and “4.65$\,{\rm\mu m}$” Aliphatic C–D Stretch Computed at the B3LYP/6-311+G∗∗ Level for Deuterated Pyrenes Containing a Methyl-Deuterated Sidegroup as well as a Periphal D Atom Shown in Figure 4.  Compound | $\lambda_{3.3}$ | $A_{3.3}$ | $\lambda_{3.4}$ | $A_{3.4}$ | $\lambda_{4.4}$ | $A_{4.4}$ | $\lambda_{4.65}$ | $A_{4.65}$ | $A_{4.65}/A_{3.3}$
---|---|---|---|---|---|---|---|---|---
|  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |
Pyre1$\\_$D1 | 3.25 | 15.30 | 3.37 | 26.61 | 4.41 | 6.52 | 4.62 | 11.63 | 0.76
Pyre1$\\_$D2 | 3.25 | 15.79 | 3.37 | 26.75 | 4.38 | 4.15 | 4.62 | 11.51 | 0.73
Pyre1$\\_$D3 | 3.25 | 15.16 | 3.37 | 26.57 | 4.41 | 7.46 | 4.62 | 11.64 | 0.77
Pyre1$\\_$D4 | 3.25 | 15.10 | 3.37 | 26.61 | 4.41 | 7.67 | 4.62 | 11.63 | 0.77
Pyre1$\\_$D5 | 3.25 | 13.54 | 3.37 | 27.97 | 4.42 | 13.25 | 4.62 | 11.65 | 0.86
Pyre1$\\_$D6 | 3.25 | 15.25 | 3.37 | 26.56 | 4.41 | 7.11 | 4.62 | 11.62 | 0.76
Pyre1$\\_$D7 | 3.25 | 15.21 | 3.37 | 26.61 | 4.41 | 6.71 | 4.62 | 11.63 | 0.76
Pyre1$\\_$D8 | 3.25 | 14.25 | 3.37 | 26.59 | 4.40 | 12.23 | 4.62 | 11.63 | 0.82
Pyre1$\\_$D9 | 3.25 | 15.34 | 3.37 | 26.61 | 4.41 | 6.31 | 4.62 | 11.63 | 0.76
Pyre2$\\_$D1 | 3.26 | 14.62 | 3.38 | 23.99 | 4.41 | 6.65 | 4.57 | 15.39 | 1.05
Pyre2$\\_$D4 | 3.26 | 14.56 | 3.38 | 24.00 | 4.41 | 6.96 | 4.57 | 15.35 | 1.05
Pyre2$\\_$D5 | 3.25 | 14.40 | 3.38 | 24.02 | 4.43 | 7.72 | 4.57 | 14.91 | 1.04
Pyre2$\\_$D6 | 3.26 | 14.51 | 3.38 | 23.99 | 4.41 | 7.22 | 4.57 | 15.39 | 1.06
Pyre2$\\_$D7 | 3.26 | 13.65 | 3.38 | 23.99 | 4.40 | 11.60 | 4.57 | 15.39 | 1.13
Pyre4$\\_$D1 | 3.25 | 15.05 | 3.37 | 22.84 | 4.41 | 6.26 | 4.62 | 11.67 | 0.78
Pyre4$\\_$D2 | 3.25 | 14.88 | 3.37 | 22.81 | 4.41 | 7.19 | 4.62 | 11.67 | 0.78
Pyre4$\\_$D3 | 3.25 | 14.87 | 3.37 | 22.82 | 4.41 | 7.32 | 4.62 | 11.68 | 0.79
Pyre4$\\_$D4 | 3.25 | 15.68 | 3.37 | 22.90 | 4.37 | 3.36 | 4.62 | 11.60 | 0.74
Pyre4$\\_$D5 | 3.25 | 14.10 | 3.37 | 24.13 | 4.43 | 9.22 | 4.62 | 11.67 | 0.83
Pyre4$\\_$D6 | 3.25 | 14.00 | 3.37 | 22.82 | 4.40 | 12.06 | 4.62 | 11.68 | 0.83
Pyre4$\\_$D7 | 3.25 | 15.05 | 3.37 | 22.83 | 4.41 | 6.09 | 4.62 | 11.67 | 0.78
Pyre4$\\_$D8 | 3.25 | 14.99 | 3.37 | 22.81 | 4.41 | 6.66 | 4.62 | 11.67 | 0.78
Pyre4$\\_$D9 | 3.25 | 14.07 | 3.37 | 22.78 | 4.40 | 11.93 | 4.62 | 11.69 | 0.83
average | 3.25 | 14.76 | 3.37 | 24.68 | 4.41 | 7.90 | 4.61 | 12.43 | 0.85
stdev | 0.00 | 0.61 | 0.00 | 1.80 | 0.01 | 2.60 | 0.02 | 1.54 | 0.12
Table 5: Wavelengths ($\lambda$) and Intensities ($A_{\lambda}$ in $\,{\rm km}\,{\rm mol}^{-1}$) of the Nominal “3.3$\,{\rm\mu m}$” Aromatic C–H Stretch, “3.4$\,{\rm\mu m}$” Aliphatic C–H Stretch, “4.4$\,{\rm\mu m}$” Aromatic C–D Stretch and “4.65$\,{\rm\mu m}$” Aliphatic C–D Stretch Computed at the B3LYP/6-311+G∗∗ Level for All the “Superdeuterated” Pyrenes Containing an H+H Pair, an H+D Pair, as well as a Periphal D Atom. Shown in Figure 5.  Compound | $\lambda_{3.3}$ | $A_{3.3}$ | $\lambda_{3.4}$ | $A_{3.4}$ | $\lambda_{4.4}$ | $A_{4.4}$ | $\lambda_{4.65}$ | $A_{4.65}$ | $A_{4.65}/A_{3.3}$
---|---|---|---|---|---|---|---|---|---
|  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |  ($\mu$m) |  (km$\,{\rm mol}^{-1}$) |
Pyre1D$\\_$D1 | 3.26 | 15.29 | 3.44 | 37.86 | 4.41 | 5.17 | 4.59 | 18.00 | 1.18
Pyre1D$\\_$D2 | 3.26 | 15.06 | 3.44 | 37.86 | 4.41 | 5.66 | 4.59 | 18.02 | 1.20
Pyre1D$\\_$D3 | 3.26 | 14.65 | 3.44 | 37.86 | 4.41 | 7.96 | 4.59 | 18.00 | 1.23
Pyre1D$\\_$D4 | 3.25 | 14.43 | 3.44 | 37.89 | 4.42 | 9.79 | 4.59 | 17.47 | 1.21
Pyre1D$\\_$D5 | 3.26 | 15.13 | 3.44 | 37.86 | 4.41 | 6.31 | 4.59 | 18.03 | 1.19
Pyre1D$\\_$D6 | 3.26 | 13.70 | 3.44 | 37.85 | 4.40 | 13.22 | 4.59 | 18.00 | 1.31
Pyre1D$\\_$D7 | 3.25 | 14.76 | 3.44 | 37.86 | 4.41 | 8.09 | 4.59 | 18.01 | 1.22
Pyre1D$\\_$D8 | 3.25 | 14.97 | 3.44 | 38.19 | 4.43 | 5.01 | 4.59 | 18.23 | 1.22
Pyre2D$\\_$D1 | 3.26 | 15.28 | 3.45 | 35.10 | 4.41 | 5.16 | 4.60 | 22.58 | 1.48
Pyre2D$\\_$D2 | 3.26 | 15.05 | 3.45 | 35.09 | 4.41 | 5.67 | 4.60 | 22.59 | 1.50
Pyre2D$\\_$D3 | 3.26 | 14.65 | 3.45 | 35.09 | 4.41 | 7.96 | 4.60 | 22.58 | 1.54
Pyre2D$\\_$D4 | 3.25 | 14.28 | 3.44 | 35.52 | 4.42 | 9.06 | 4.60 | 22.60 | 1.58
Pyre2D$\\_$D5 | 3.26 | 15.13 | 3.45 | 35.09 | 4.41 | 6.35 | 4.60 | 22.58 | 1.49
Pyre2D$\\_$D6 | 3.26 | 13.69 | 3.45 | 35.09 | 4.40 | 13.21 | 4.60 | 22.59 | 1.65
Pyre2D$\\_$D7 | 3.25 | 14.75 | 3.45 | 35.09 | 4.41 | 8.10 | 4.60 | 22.58 | 1.53
Pyre2D$\\_$D8 | 3.25 | 15.11 | 3.45 | 35.05 | 4.43 | 5.61 | 4.60 | 22.25 | 1.47
Pyre4D$\\_$D1 | 3.25 | 15.20 | 3.41 | 26.54 | 4.42 | 7.47 | 4.59 | 18.20 | 1.20
Pyre4D$\\_$D2 | 3.25 | 15.65 | 3.41 | 26.08 | 4.41 | 7.01 | 4.59 | 18.16 | 1.16
Pyre4D$\\_$D3 | 3.25 | 15.64 | 3.41 | 26.09 | 4.41 | 7.00 | 4.59 | 18.17 | 1.16
Pyre4D$\\_$D4 | 3.25 | 15.11 | 3.41 | 26.08 | 4.40 | 8.74 | 4.59 | 18.16 | 1.20
Pyre4D$\\_$D5 | 3.25 | 15.47 | 3.41 | 26.09 | 4.41 | 8.00 | 4.59 | 18.16 | 1.17
Pyre4D$\\_$D6 | 3.25 | 15.48 | 3.41 | 26.08 | 4.41 | 8.01 | 4.59 | 18.16 | 1.17
Pyre4D$\\_$D7 | 3.25 | 15.11 | 3.41 | 26.09 | 4.40 | 8.71 | 4.59 | 18.17 | 1.20
Pyre4D$\\_$D8 | 3.25 | 15.36 | 3.41 | 26.12 | 4.42 | 8.23 | 4.59 | 17.64 | 1.15
Average | 3.25 | 14.96 | 3.43 | 33.06 | 4.41 | 7.73 | 4.59 | 19.54 | 1.31
Stdev | 0.00 | 0.52 | 0.02 | 5.13 | 0.01 | 2.16 | 0.00 | 2.18 | 0.17
Table 6: Mean Intensities of the 3.3$\,{\rm\mu m}$ Aromatic C–H Stretch ($A_{3.3}$) and the 4.65$\,{\rm\mu m}$ Aliphatic C–D Stretch ($A_{4.65}$) Computed at the B3LYP/6-311+G∗∗ Level for All the Deuterated PAHs Considered in This Work. The $A_{3.3}$ and $A_{4.65}$ Band Strengths Are on a Per C–H or C–D Bond Basis.  Band Strengths | Neutral PAHs | Cationic PAHs
---|---|---
$A_{3.3}$ (km$\,{\rm mol}^{-1}$) |  14.76$\pm{0.33}$ |  2.14${}^{+3.08}_{-2.14}$
$A_{4.65}$ (km$\,{\rm mol}^{-1}$) |  15.26$\pm{3.31}$ |  1.64$\pm{0.75}$
$A_{4.65}/A_{3.3}$ |  1.04$\pm{0.21}$ |  1.67$\pm{1.00}$

Table 7: A Summary of the Power Emitted from the 3.3$\,{\rm\mu m}$ Aromatic
C–H Stretch $\left(I_{3.3}\right)_{\rm obs}$, the 3.4$\,{\rm\mu m}$ Aliphatic
C–H Stretch $\left(I_{3.4}\right)_{\rm obs}$, the 4.4$\,{\rm\mu m}$ Aromatic
C–D Stretch $\left(I_{4.4}\right)_{\rm obs}$, and the 4.65$\,{\rm\mu m}$
Aliphatic C–D Stretch $\left(I_{4.65}\right)_{\rm obs}$ Compiled from
Astronomical Observations. Also Shown Are the Degree of Deuteration Derived
for Each Source

Object | Type | aromatic C–H |  | aliphatic C–H |  | aromatic C–D |  | aliphatic C–D
---|---|---|---|---|---|---|---|---
|  | $\left(I_{3.3}\right)_{\rm obs}^{\dagger}$ |  | $\left(I_{3.4}\right)_{\rm obs}^{\dagger}$ | $\left(I_{3.4}/I_{3.3}\right)_{\rm obs}^{\dagger}$ | $f_{\rm aliCH}$ |  | $\left(I_{4.4}\right)_{\rm obs}^{\dagger}$ | $\left(I_{4.4}/I_{3.3}\right)_{\rm obs}^{\dagger}$ | $f_{\rm aroCD}$ |  | $\left(I_{4.65}\right)_{\rm obs}^{\dagger}$ | $\left(I_{4.65}/I_{3.3}\right)_{\rm obs}^{\dagger}$ | $f_{\rm aliCD}$
G75.78+0.34a |  Hii Region |  22.2 $\pm{2.0}$ |  |  13.3 $\pm{2.1}$ |  59.9 $\pm{4.1}$ |  34.0 $\pm{4.7}$ |  |  0.30 $\pm{0.19}$ |  1.35 $\pm{0.73}$ |  3.3 $\pm{1.3}$ |  |  2.25 $\pm{0.86}$ |  10.1 $\pm{2.96}$ |  7.1 $\pm{4.1}$
NGC3603a |  Hii Region |  25.1 $\pm{2.4}$ |  |  16.7 $\pm{3.1}$ |  66.5 $\pm{6.0}$ |  37.8 $\pm{5.7}$ |  |  0.62 $\pm{0.52}$ |  2.47 $\pm{1.84}$ |  5.9 $\pm{2.8}$ |  |  8.70 $\pm{4.41}$ |  34.7 $\pm{14.3}$ |  20.6 $\pm{16.0}$
W51 obs.1a |  Hii Region |  20.0 $\pm{1.8}$ |  |  11.2 $\pm{1.5}$ |  56.0 $\pm{2.5}$ |  31.8 $\pm{4.0}$ |  |  0.41 $\pm{0.14}$ |  2.05 $\pm{0.52}$ |  5.0 $\pm{1.5}$ |  |  0.88 $\pm{0.47}$ |  4.40 $\pm{1.95}$ |  3.2 $\pm{2.1}$
W51 obs.2a |  Hii Region |  20.6 $\pm{1.8}$ |  |  11.5 $\pm{1.4}$ |  55.8 $\pm{1.9}$ |  31.7 $\pm{3.9}$ |  |  0.16 $\pm{0.14}$ |  0.78 $\pm{0.61}$ |  1.9 $\pm{0.9}$ |  |  0.86 $\pm{0.47}$ |  4.17 $\pm{1.92}$ |  3.0 $\pm{2.0}$
M8a |  Hii Region |  99.1 $\pm{8.9}$ |  |  55.2 $\pm{6.8}$ |  55.7 $\pm{1.9}$ |  31.6 $\pm{3.9}$ |  |  1.03 $\pm{0.53}$ |  1.04 $\pm{0.44}$ |  2.6 $\pm{0.9}$ |  |  1.87 $\pm{1.63}$ |  1.89 $\pm{1.48}$ |  1.4 $\pm{1.3}$
IRAS112073 obs.1a |  Hii Region |  11.7 $\pm{1.0}$ |  |  6.7 $\pm{0.8}$ |  57.3 $\pm{1.9}$ |  32.5 $\pm{4.0}$ |  |  0.30 $\pm{0.08}$ |  2.56 $\pm{0.46}$ |  6.1 $\pm{1.8}$ |  |  4.40 $\pm{0.76}$ |  37.6 $\pm{3.28}$ |  22.0 $\pm{13.1}$
IRAS112073 obs.2a |  Hii Region |  12.2 $\pm{1.0}$ |  |  6.1 $\pm{0.7}$ |  50.0 $\pm{1.6}$ |  28.4 $\pm{3.5}$ |  |  0.37 $\pm{0.08}$ |  3.03 $\pm{0.41}$ |  7.2 $\pm{2.1}$ |  |  2.92 $\pm{0.76}$ |  23.9 $\pm{4.27}$ |  15.2 $\pm{8.8}$
M17ba | PDR |  12.3 $\pm{1.5}$ |  |  6.6 $\pm{1.8}$ |  53.7 $\pm{8.1}$ |  30.5 $\pm{5.9}$ |  | \- - | \- - | \- - |  |  3.32 $\pm{1.94}$ |  27.0 $\pm{12.5}$ |  16.8 $\pm{13.2}$
M17b | PDR |  14.9 $\pm{0.2}$ |  |  4.8 $\pm{0.1}$ |  32.2 $\pm{0.2}$ |  18.3 $\pm{2.2}$ |  |  0.14 $\pm{0.04}$ |  0.94 $\pm{0.26}$ |  2.3 $\pm{0.7}$ |  |  0.32 $\pm{0.06}$ |  2.15 $\pm{0.37}$ |  1.6 $\pm{0.8}$
Orion Barb | PDR |  53.1 $\pm{0.6}$ |  |  19.2 $\pm{0.3}$ |  36.2 $\pm{0.2}$ |  20.5 $\pm{2.5}$ |  |  1.36 $\pm{0.11}$ |  2.56 $\pm{0.18}$ |  6.1 $\pm{1.7}$ |  |  0.76 $\pm{0.13}$ |  1.43 $\pm{0.23}$ |  1.1 $\pm{0.5}$
G18.14.0b | Reflection Nebula |  4.78 $\pm{0.06}$ |  |  1.99 $\pm{0.04}$ |  41.6 $\pm{0.3}$ |  23.7 $\pm{2.8}$ |  |  0.12 $\pm{0.01}$ |  2.45 $\pm{0.18}$ |  5.9 $\pm{1.6}$ |  |  0.03 $\pm{0.01}$ |  0.56 $\pm{0.20}$ |  0.4 $\pm{0.2}$

† Flux densities in units of $10^{-17}$ W m-2 arcsec-2

a Doney et al. 2016

b Onaka et al. 2014