# Non-adiabatic decay of Rydberg-atom-ion molecules

A. Duspayev<EMAIL_ADDRESS>Department of Physics, University of Michigan,
Ann Arbor, MI 48109, USA  G. Raithel  Department of Physics, University of
Michigan, Ann Arbor, MI 48109, USA

###### Abstract

The decay of Rydberg-atom-ion molecules (RAIMs) due to non-adiabatic couplings
between electronic potential energy surfaces is investigated. We employ the
Born-Huang representation and perform numerical simulations using a Crank-
Nicolson algorithm. The non-adiabatic lifetimes of rubidium RAIMs for the
lowest ten vibrational states, $\nu$, are computed for selected Rydberg
principal quantum numbers, $n$. The non-adiabatic lifetimes are found to
generally exceed the radiative Rydberg-atom lifetimes. We observe and explain
a trend of the lifetimes as a function of $\nu$ and $n$, and attribute
irregularities to quantum interference arising from a shallow potential well
in an inner potential surface. Our results will be useful for future
spectroscopic studies of RAIMs.

##  I Introduction

Ultralong-range Rydberg molecules (ULRM) Shaffer  _et al._ (2018); Fey  _et
al._ (2020) are an active direction in Rydberg-atom research. ULRMs can be
distinguished based upon their formation mechanisms. For instance, in Rydberg-
ground molecules Greene  _et al._ (2000); Bendkowsky  _et al._ (2009) a
ground-state atom resides within the Rydberg-atom wave-function, and a
molecular bond is formed due to scattering of the Rydberg electron at the
perturber atom. The studies on this type of ULRMs include detailed analyses of
association Hamilton  _et al._ (2002); Bendkowsky  _et al._ (2010); Tallant
_et al._ (2012); Bellos  _et al._ (2013); Anderson  _et al._ (2014a); DeSalvo
_et al._ (2015); Peper and Deiglmayr (2021), electronic structure Anderson
_et al._ (2014b), spin-orbit coupling Kleinbach  _et al._ (2017); Eiles and
Greene (2017); Deiß  _et al._ (2020) and scattering processes Saßmannshausen
_et al._ (2015); Böttcher  _et al._ (2016); Engel  _et al._ (2019); MacLennan
_et al._ (2019), calculations and measurements of lifetimes Butscher  _et al._
(2010, 2011); Camargo  _et al._ (2016) and permanent electric dipole moments
Khuskivadze  _et al._ (2002) that can vary from a few Li  _et al._ (2011); Bai
_et al._ (2020) to thousands of Debyes Booth  _et al._ (2015); Niederprüm  _et
al._ (2016), and interactions with external fields Lesanovsky  _et al._
(2006); Hummel  _et al._ (2019, 2021). In another type of Rydberg molecules,
referred to as macrodimers Boisseau  _et al._ (2002); Overstreet  _et al._
(2009); Saßmannshausen and Deiglmayr (2016), two Boisseau  _et al._ (2002) or
more Samboy and Côté (2013) Rydberg atoms with non-overlapping wave-functions
(LeRoy radius condition Le Roy (1974)) become bounded via multipolar
interactions Singer  _et al._ (2005); Schwettmann  _et al._ (2006); Deiglmayr
_et al._ (2014); Marcassa and Shaffer (2014). Their formation Samboy  _et al._
(2011); Kiffner  _et al._ (2012); Han  _et al._ (2018, 2019), vibrational
structure Hollerith  _et al._ (2019), lifetimes Schwettmann  _et al._ (2007);
Overstreet  _et al._ (2009); Saßmannshausen and Deiglmayr (2016) and alignment
with external fields Hollerith  _et al._ (2021) have been studied.

Recently, a Rydberg-atom-ion molecule (RAIM) Duspayev  _et al._ (2021); Deiß
_et al._ (2021) has been proposed that opens new perspectives at the interface
between the fields of Rydberg molecules and atom-ion interactions Schmid  _et
al._ (2010); Secker  _et al._ (2016, 2017); Schmid  _et al._ (2018); Ewald
_et al._ (2019); Wang  _et al._ (2020); Hirzler and Pérez-Ríos (2021);
Dieterle  _et al._ (2021). In RAIMs, multipolar interaction between a Rydberg
atom and an ion outside of the atom leads to bound molecular states. The non-
adiabatic decay rate of RAIMs was predicted to be negligibly small Duspayev
_et al._ (2021); Deiß  _et al._ (2021), based upon Landau-Zener (LZ) tunneling
probabilities. Since the assumptions of LZ tunneling are not satisfied in
RAIMs, a quantum theory on the non-adiabatic decay of RAIMs is needed.

Here, we develop a quantum model for non-adiabatic RAIM decay utilizing the
Born-Huang representation (BHR) Born and Huang (1954), in which the
vibrational motion is treated fully quantum-mechanically and non-adiabatic
couplings are accurately described. Being a common method to study non-
adiabatic processes in conventional molecules Agostini and Curchod (2019), the
BHR has not yet been applied to ULRMs, to our knowledge. After reviewing the
basic theory of RAIMs in Sec. II, we describe the BHR model for their non-
adiabatic dynamics in Sec. III. Results for selected RAIMs are presented and
discussed in Sec. IV. The paper is concluded in Sec. V.

FIG. 1: (a) Sketch of Rydberg-atom-ion molecule (RAIM). (b) PECs of rubidium
RAIMs for $m_{J}=1/2$ and 3/2 as a function of internuclear distance $R$.
Energies are relative to the field-free atomic $45P_{3/2}$ state. The wells in
the regions A, B and C are expected to support RAIMs. (c) Magnified view of
the region A with the three $m_{J}=1/2$-PECs used in our calculation of non-
adiabatic decay, and wave-function densities of the five lowest RAIM
vibrational states in the PEC labeled $V_{3}$. The vertical offsets of the
wave-function densities are for clarity and are not related with vibrational
energy. The ”X” marks the most relevant anti-crossing. The imaginary
absorbing-wall potential $W_{abs}$ used in the computational approach is also
shown.

##  II Theory of Rydberg-atom-ion molecules

The theory of RAIMs has been developed in Duspayev  _et al._ (2021). RAIMs,
sketched in Fig. 1(a), are formed between an ion and a neutral Rydberg atom
via electric-multipole interaction. The internuclear distance $R$ is larger
than the radius of the Rydberg atom. Adopting a $z-$axis aligned with the
internuclear axis and assuming a point-like positive ion, the interaction is,
in atomic units Schwettmann  _et al._ (2006); Deiglmayr  _et al._ (2014); Han
_et al._ (2018),

$\displaystyle
V_{int,m_{J}}(\mathbf{\hat{r}}_{e};R)=-\sum_{l=1}^{l_{max}}\sqrt{\frac{4\pi}{2l+1}}\frac{\hat{r}_{e}^{l}}{R^{l+1}}Y_{l0}(\hat{\theta}_{e},\hat{\phi}_{e}).$
(1)

Here, $m_{J}$ is the conserved magnetic quantum number of the Rydberg atom,
$n$ the principal quantum number, $\ell$ the orbital quantum number,
$\mathbf{\hat{r}}_{e}=(\hat{r}_{e},\hat{\theta}_{e},\hat{\phi}_{e})$ the
Rydberg-electron position operator relative to the atom’s center, and
$Y_{l0}(\hat{\theta}_{e},\hat{\phi}_{e})$ a spherical harmonic in which $l$ is
the multipole order of the atomic charge distribution. Diagonalization of the
Hamiltonian with atom-ion interaction given in Eq. 1 yields the molecular
potential energy curves (PECs). Some PECs exhibit deep wells conducive to
bound vibrational states of RAIMs, such as PECs in the vicinity of $nP_{J}$
Rydberg states of cesium and rubidium. The characteristics of these RAIMs
depend on quantum defects and other atomic parameters Duspayev  _et al._
(2021); Deiß  _et al._ (2021).

In Fig. 1(b), we show a case in which RAIMs are formed below the Rb $45P$
asymptotes. The RAIM potential wells are several hundreds of MHz deep and on
the order of 100 nm wide, which leads to tens of bound vibrational states. In
Fig. 1(c), we show the lowest 5 RAIM states in the molecular PEC labeled
$V_{3}$. The stability of the vibrational states may be affected by the non-
adiabatic couplings between $V_{3}$ and neighboring PECs, labeled $V_{1}$ and
$V_{2}$ in Fig. 1(c).

##  III Non-adiabatic dynamics

###  III.1 Time-dependent Schrödinger equation

In the Born-Oppenheimer approximation (BOA) Born and Oppenheimer (1927),
electronic and nuclear wave-functions are adiabatically separated to
facilitate the calculation of vibrational molecular states. However, the BOA
can ”break down” when vibrational and electronic time scales approach each
other, as is the case when PECs exhibit narrow anti-crossings, leading to non-
adiabatic coupling of PECs and to molecular decay. In Fig. 1(c), non-adiabatic
coupling from $V_{3}$ to the unbound PECs may contribute significantly to RAIM
decay.

A common method to model non-adiabatic effects is to write the TDSE in BHR
Born and Huang (1954) (see Agostini and Curchod (2019) and references therein
for a recent overview of BHR theory and applications). Here, we consider the
vibrational degree of freedom of a diatomic RAIM along its internuclear axis,
$\hat{\mathbf{R}}$. The TDSE in BHR is written as:

$\displaystyle i\hbar\frac{\partial\psi_{i}(R,t)}{\partial
t}=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}\psi_{i}(R,t)}{\partial
R^{2}}+V_{i}(R)\psi_{i}(R,t)$
$\displaystyle+\sum_{j}\mathcal{F}_{ij}(R)\psi_{j}(R,t).$ (2)

Here, $\psi_{i}(R,t)$ is the adiabatic vibrational RAIM wave-function on
$V_{i}(R)$, $\mu$ is the reduced mass, and the $\mathcal{F}_{ij}(R)$ are the
non-adiabatic couplings between the adiabatic wave-functions on PECs $V_{i}$
and on $V_{j}$. The non-adiabatic terms of the BHR, $\mathcal{F}_{ij}(R)$,
couple nuclear and electronic motion. Explicitly,

$\displaystyle\mathcal{F}_{ij}(R)=\mathcal{A}_{ij}(R)\cdot\frac{\partial}{\partial
R}+\mathcal{B}_{ij}(R),$ (3)

where $\mathcal{A}_{ij}(R)$ is referred to as the first-order non-adiabatic
coupling and is defined as (in one dimension):

$\displaystyle\mathcal{A}_{ij}(R)=-\frac{\hbar^{2}}{\mu}\left\langle\phi_{i}\right|\frac{\partial}{\partial
R}\left|\phi_{j}\right\rangle\quad.$ (4)

There, $|\phi_{i}(R)\rangle$ are the $R$-dependent electronic states of the
Rydberg atom, with wave-functions
$\phi_{i}({\bf{r}}_{e};R)=\langle{\bf{r}}_{e}|\phi_{i}(R)\rangle$. The inner
product in Eq. 4 is evaluated in the Rydberg-electron state space, i. e. it
involves, in principle, an integral over ${\bf{r}}_{e}$. In practice, the
Rydberg state is represented in the “diabatic” Rydberg-state basis
$\\{|n,\ell,J,m_{J}\rangle=:|\alpha\rangle\\}$, with the shorthand $\alpha$
for all diabatic-state quantum numbers. The electronic wave-function on PEC
$i$ then reads

$\phi_{i}({\bf{r}}_{e};R)=\sum_{\alpha}\quad
c_{i,\alpha}(R)\langle{\bf{r}}_{e}|\alpha\rangle\quad,$

with coefficient functions $c_{i,\alpha}(R)$. Because of the $R$-independence
and the orthonormality of the $|\alpha\rangle$,

$\displaystyle\mathcal{A}_{ij}(R)=-\frac{\hbar^{2}}{\mu}\sum_{\alpha}\quad
c^{*}_{i,\alpha}(R)\left[\frac{\partial}{\partial
R}c_{j,\alpha}(R)\right]\quad.$ (5)

The $\mathcal{B}_{ij}(R)$ in Eq. 3 is referred to as the second-order non-
adiabatic coupling and is

$\displaystyle\mathcal{B}_{ij}(R)$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{2\mu}\left\langle\phi_{i}\right|\frac{\partial^{2}}{\partial
R^{2}}\left|\phi_{j}\right\rangle$ (6) $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{2\mu}\sum_{\alpha}\quad
c^{*}_{i,\alpha}(R)\left[\frac{\partial^{2}}{\partial
R^{2}}c_{j,\alpha}(R)\right]\quad.$

Since $\mathcal{A}_{ij}(R)=-\mathcal{A}_{ji}(R)$, it is
$\mathcal{A}_{ii}(R)=0$. The remaining (generally non-zero) diagonal non-
adiabatic couplings, $\mathcal{B}_{ii}(R)$, are often combined with the
corresponding $V_{i}(R)$ into

$\displaystyle\tilde{V}_{i}(R)=V_{i}(R)+\mathcal{B}_{ii}(R)\quad.$ (7)

We may then re-write Eqs. 2 and 3 as

$\displaystyle i\hbar\frac{\partial\psi_{i}(R,t)}{\partial
t}=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}\psi_{i}(R,t)}{\partial
R^{2}}+\tilde{V}_{i}(R)\psi_{i}(R,t)$ $\displaystyle+\sum_{j\neq
i}\Big{[}\mathcal{A}_{ij}(R)\frac{\partial\psi_{j}(R,t)}{\partial
R}+\mathcal{B}_{ij}(R)\psi_{j}(R,t)\Big{]},$ (8)

We refer to $\tilde{V}_{i}(R)$ as “adiabatic potentials”, i.e. potentials in
which all diagonal non-adiabatic energy shifts of the adiabatic states have
been added to the PECs, $V_{i}(R)$. We also note that for $\mathcal{B}$ it is

$\mathcal{B}_{ij}+\mathcal{B}_{ji}^{*}=-\frac{\hbar^{2}}{\mu}\big{[}\frac{\partial}{\partial
R}\langle\phi_{i}|\big{]}\big{[}\frac{\partial}{\partial
R}|\phi_{j}\rangle\big{]}\quad.$

We use this identity as a check for numerical errors caused by the step size
in $R$.

It is apparent from Eqs. 4-6 that the non-adiabatic couplings follow from the
$R$-dependencies of the adiabatic Rydberg states $\phi_{i}({\bf{r}}_{e};R)$,
which are critically affected by the avoided crossings between the PECs. As
the general shapes of the PECs are the same for all cases studied here, in
Sec. IV we find a general trend for the non-adiabatic RAIM decay times.
However, detailed differences in the non-adiabatic $\mathcal{A}$\- and
$\mathcal{B}$-terms and in the wave-function dynamics on the dissociative PECs
lead to peculiar quantum effects that are also discussed.

###  III.2 Simulation method

We first numerically calculate PECs for selected Rydberg states $nP_{J}$ of
87Rb (see, for instance, Fig. 1(b) for $45P_{J}$). We investigate bound RAIM
vibrational states on PECs that asymptotically connect with the $nP_{1/2}$
levels. The corresponding PEC for $n=45$ is labeled $V_{3}$ in Fig. 1(c). This
PEC, and equivalent PECs for other $n$-values, exhibit non-adiabatic couplings
mostly to a pair of lower, dissociating PECs labeled $V_{1}$ and $V_{2}$ in
Fig. 1(c). The PEC calculation yields the PECs, the associated adiabatic
Rydberg states $|\phi_{i}(R)\rangle$, and the non-adiabatic terms
$\mathcal{A}_{ij}(R)$ and $\mathcal{B}_{ij}(R)$ according to the equations in
Sec. III.1. In view of the structure of the PEC anti-crossings evident in Fig.
1, the effect of non-adiabatic couplings of $V_{3}$ to PECs other than $V_{1}$
and $V_{2}$ is deemed negligible. The Rydberg-state basis sets used in the PEC
calculations include all Rydberg levels with $m_{J}=1/2$ and effective
principal quantum number differing by less than 5 from that of the molecular
RAIM states of interest.

The computational method to solve the TDSE in Eq. 8 on the three relevant PECs
$V_{1}$, $V_{2}$ and $V_{3}$ is a Crank-Nicolson (CN) algorithm Koonin and
Meredith (1990). The simulation is initialized with a RAIM vibrational state
on the adiabatic potential $\tilde{V}_{3}$ (which differs slightly from the
PEC $V_{3}$, according to Eq. 7). Over the course of the subsequent simulated
evolution, the norm of the wave-function decays due to the non-adiabatic
couplings, allowing us to extract the molecular lifetimes. In the following,
we describe additional details of the method.

The potential $\tilde{V_{3}}(R)$, constructed according to Eq. 7, is used to
calculate the initial vibrational RAIM state, $\Psi_{\nu}(R)$, with
vibrational quantum number $\nu$. The initial state for the CN simulation then
is $\psi_{1}(R,t=0)=\psi_{2}(R,t=0)=0$ and $\psi_{3}(R,t=0)=\Psi_{\nu}(R)$.
The wave-function is propagated in time for a duration of
$t_{total}=50~{}\mu$s with a step size $\Delta t=20$ ps. In order to reduce
transients from sudden “turn-on” of the non-adiabatic terms, the non-adiabatic
terms $\mathcal{A}_{ij}(R)$ and $\mathcal{B}_{ij}(R)$ are slowly ramped up at
the beginning of the time-propagation. We still find minor initial transients
in the $\psi_{i}(R,t)$, which cease at times $t_{0}\sim 10$ ns.

The dissociating potentials, $V_{1}$ and $V_{2}$, are unbound. As our CN
simulation employs a spatial box with fixed boundary conditions
$\psi_{i}(R)=0$ on all boundaries, the potentials must be modified such that
wave-functions propagating outward on $V_{1}$ and $V_{2}$ are absorbed rather
than reflected. To terminate the outgoing wave-function, we add an imaginary
part, $Im[V_{i}]=W_{i,abs}(R)$, on the unbound potentials $V_{1}$ and $V_{2}$,
as depicted by the dashed line in Fig. 1(c). The domain over which
$W_{i,abs}(R)$ differs from zero is placed far enough out in $R$ that it does
not affect the non-adiabatic dynamics of interest, which is restricted to
regions within which the non-adiabatic couplings differ from zero. The
absorbing wall $W_{i,abs}(R)$ exhibits a smooth turn-on, so as to avoid
reflections. We have checked the effectiveness of the absorbing wall as well
as the absence of wall reflections by calculating the quantum flux as a
function of $R$ (outside and inside the wall), and by verifying the absence of
standing-wave patters on $\psi_{1}$ and $\psi_{2}$ near the absorbing wall.

The absorbed outgoing flux leads to a decay of the overall wave-function norm
Muga  _et al._ (2004), allowing us to extract the RAIM lifetime. The
population in $\psi_{3}$, $p_{3}(t)=\int|\psi_{3}(R,t)|^{2}dR$, is determined
as a function of propagation time and fitted to the function:

$\displaystyle p_{3}(t)=p_{3}(t_{0})e^{-(t-t_{0})/\tau_{nad}},$ (9)

with fitting parameters $p_{3}(t_{0})$ and $\tau_{nad}$. Here,
$p_{3}(t_{0})\lesssim 1$ reflects the population after ramping up the non-
adiabatic terms and after allowing transients to cease, and $\tau_{nad}$ is
the non-adiabatic RAIM lifetime for the given $n$ and $\nu$.

Although most $\tau_{nad}$-values are longer than $t_{total}$, as seen in Fig.
2, the decrease of $p_{3}(t)$ during the interval $t_{total}$ allows for an
accurate determination of $\tau_{nad}$ in all cases studied. We have checked
that lowering the computation time step $\Delta t$ does not significantly
alter the $\tau_{nad}$.

FIG. 2: (a) Calculated non-adiabatic lifetimes, $\tau_{nad}$, of rubidium
RAIMs below the $nP_{1/2}$ asymptotes vs vibrational and principal quantum
numbers, $\nu$ and $n$, displayed on the indicated logarithmic color scale.
The X-marks and dashed circles mark cases in which $\tau_{nad}$ is less than
ten times the radiative decay time at 1.2 K and 300 K, respectively. (b) The
lifetimes $\tau_{nad}$ for the ground vibrational state, $\nu=0$,
corresponding to the first column in (a). The insets show the computed
population decrease, $\Delta p_{3}(t)$, of the ground RAIM states for $n=35$
(inset (A)) and $n=65$ (inset (B)) as a function of time. The
$\tau_{nad}$-values follow from the slopes according to Eq. 9.

##  IV Results and Discussion

We obtain the non-adiabatic lifetimes, $\tau_{nad}$, of the lowest ten
vibrational states of RAIMs below the $nP_{1/2}$ Rydberg-state asymptotes for
$n\in\\{25,35,45,55,65\\}$. The results are listed in Table 1 and visualized
in Fig. 2(a). These lifetimes are much longer than the radiative lifetimes of
$nP$ Rydberg states, $\tau_{r}$, which have values between
$\tau_{r}=32~{}\mu$s for $n=25$ and 680 $\mu$s for $n=65$ in a 1.2-K black-
body radiation field. In 300 K radiation fields, the $\tau_{r}$-values are
between 17 $\mu$s and 160 $\mu$s for $n=25$ and 65, respectively, because
upward and downward bound-bound transitions as well as thermal ionization
reduce the lifetime Gallagher (2005). Our $\tau_{r}$-values, obtained in
context with work involving Rydberg-atom-state diffusion in thermal radiation
fields Traxler  _et al._ (2013); Anderson  _et al._ (2013), are roughly in-
line with values of $\sim 20\mu$s at $n=25$ Branden  _et al._ (2009) and $\sim
150\mu$s at $n=65$ Beterov  _et al._ (2009) reported elsewhere.

Under the absence of other decay channels, the net RAIM decay time
$\tau=(1/\tau_{nad}+1/\tau_{r})^{-1}$. The symbols on the color map in Fig.
2(a) mark cases in which $\tau_{nad}<10\tau_{r}$ in 1.2 K and 300 K black-body
fields. Figure 2(a) shows that non-adiabatic RAIM decay, while not being the
dominant decay mechanism, should have a noticeable effect at higher $n$-values
and vibrational quantum numbers $\nu$. We note that the values for $\tau_{r}$
assumed in Fig. 2(a) are for $nP$ Rydberg states, whereas the electronic
states of RAIMs states carry up to about 50% admixture of longer-lived
high-$\ell$ Rydberg levels. This means non-adiabatic decay might be slightly
more relevant, on a relative scale, than suggested in Fig. 2(a).

Table 1:  Calculated non-adiabatic lifetimes $\tau_{nad}$ in seconds. $\nu$ | $n=25$ | $n=35$ | $n=45$ | $n=55$ | $n=65$
---|---|---|---|---|---
0 | 4.99 | 1.91 | 0.129 | 2.180 | 0.0595
1 | 1.67 | 0.116 | 0.0627 | 0.0175 | 0.0724
2 | 0.123 | 0.0393 | 0.0105 | 0.0425 | 0.00570
3 | 0.238 | 0.0216 | 0.0167 | 0.00297 | 0.00659
4 | 0.0367 | 0.00773 | 0.00324 | 0.00855 | 0.00259
5 | 0.103 | 0.00685 | 0.00631 | 0.00122 | 0.00103
6 | 0.0151 | 0.00320 | 0.00219 | 0.00135 | 0.00426
7 | 0.0295 | 0.00221 | 0.00132 | 0.00162 | 0.000475
8 | 0.00920 | 0.00262 | 0.0143 | 0.000411 | 0.000534
9 | 0.00702 | 0.000851 | 0.000572 | 0.00120 | 0.00139

To illustrate the inadequacy of LZ estimates for non-adiabatic RAIM decay, we
consider the RAIM vibrational ground state for $45P_{J}$, for which we have
computed $\tau_{nad}\sim 0.13$ s (see Fig. 2(b) and Table 1). From the PECs in
Fig. 1 and the vibrational energy levels we estimate a vibration frequency of
$f_{0}=18.2$ MHz, corresponding to a LZ decay “attempt rate” of
$R=2f_{0}=36.4$ MHz. The main avoided crossing has a gap size of $G\approx
h\times$350 MHz (see gap marked “X” in Fig. 1(c)). For a LZ RAIM decay
estimate, we assume a fixed particle velocity given by the vibration velocity,
$v_{max}$, at the minimum of $V_{3}$. From $(\mu/2)v_{max}^{2}=hf_{0}/2$, with
effective mass $\mu=43.5$ amu for 87Rb, we obtain $v_{max}=0.41$ m/s. The
differential slope of the level crossing, estimated from Fig. 1(b), is
$s=h\times 49$ GHz/$\mu$m. The LZ tunneling probability then is
$P_{LZ}=\exp(-2\pi(G/2)^{2}/\hbar/s/v_{max})=7\times 10^{-27}$, and the LZ
RAIM lifetime $\tau=1/(RP_{LZ})=4\times 10^{18}$ s. This estimate is about 20
orders of magnitude too large.

FIG. 3: Wave-function densities $|\psi_{1}|^{2}$ and $|\psi_{3}|^{2}$ in the
respective potentials $V_{1}$ and $V_{3}$ after reaching a quasi steady-state
for $55P_{1/2}$. The inset shows an enlarged view of $|\psi_{1}|^{2}$.

From Fig. 2 and Table 1 it is seen that the non-adiabatic lifetimes follow a
downward trend with increasing $n$ and increasing $\nu$. This is expected
because the anti-crossing gaps decrease with increasing $n$, and because the
$\partial/\partial R$-operator, which occurs in combination with the
$\mathcal{A}$-terms, exacerbates the non-adiabatic coupling at higher
$\nu$-values, where the vibrational wave-function gradients become larger.

Inspecting Fig. 2 and Table 1, we further note substantial deviations of
$\tau_{nad}$ from the overall trend. E.g., $\tau_{nad}$ for $\nu=0$ has an
outlier at $n=55$, for which $\tau_{nad}$ is larger than even for $n=35$. We
believe that the irregularities originate in a quantum-interference effect in
the RAIM decay. To illustrate this, in Fig. 3 we show PECs $V_{1}$ and $V_{3}$
and their adiabatic wave-function densities. The non-adiabatic couplings set
up a quasi-stationary $\psi_{1}$ and an associated probability-density flow
that causes the non-adiabatic RAIM decay. The density $|\psi_{1}|^{2}$
relative to $|\psi_{3}|^{2}$ is small, but non-zero (even at the right margin
of the plot), and it exhibits a standing wave in a shallow well in $V_{1}$
centered around $\approx 2.20~{}\mu$m. The standing wave is due to a 100-%
reflection on the rising side of $V_{1}$ near 2.18 $\mu$m, and a partial
quantum reflection at the three-level crossing near 2.23 $\mu$m. The net
outward flow on $V_{1}$ to the right of the anti-crossing is a superposition
of a contribution due to direct non-adiabatic coupling from $V_{3}$ onto
$V_{1}$, and a contribution that proceeds via non-adiabatic coupling from
$V_{3}$ into the shallow potential well in $V_{1}$ and time-delayed escape, as
visualized by the delay loop in Fig. 3. The superposition amplitude depends on
the phase difference between the contributions, which varies as a function of
$n$ and $\nu$, causing the irregularities seen in Fig. 2. A few combinations
of $n$ and $\nu$ appear to exhibit substantial destructive interference,
leading to lifetimes that are much longer than the overall trend would
suggest.

The situation portrayed in Fig. 3 applies to all values of $n$ studied. Our
interpretation of the lifetime irregularities in Fig. 2 in terms of a quantum
interference effect has been supported in additional test calculations, not
shown, in which an absorbing potential has been placed within the shallow
potential well in $V_{1}$. The test calculations show a smooth dependence of
$\tau_{nad}$ without irregularities.

##  V Conclusion

In summary, we have presented results of calculations of non-adiabatic decay
of Rydberg-atom-ion molecules of Rb. The lifetimes, extracted for five
representative values of $n$ and for the lowest ten vibrational states, follow
an overall trend that is in-line with the behavior of avoided-crossing gap
sizes and the structure of the vibrational wave-functions. Deviations from the
trend were attributed to a quantum interference effect. Our results confirm
that the RAIM states are quite stable against non-adiabatic decay, and that
their lifetimes are mainly limited by radiative decay of the Rydberg valence
electron. Future experimental work may reveal the existence and the lifetimes
of RAIM molecules. Future computational work could be devoted to non-adiabatic
decay of RAIMs in other potential wells evident in Fig. 1(b), as well as a
more detailed study of quantum interference effects in RAIM decay. The
formalism discussed here could also be applied to study non-adiabatic
processes in other types of Rydberg molecules.

## ACKNOWLEDGMENTS

We thank Bineet Kumar Dash and Ansh Shah for useful discussions. This work was
supported by NSF grant No. PHY-1806809 and in part through computational
resources and services provided by Advanced Research Computing at the
University of Michigan, Ann Arbor.

## References

  * Shaffer  _et al._ (2018) J. P. Shaffer, S. T. Rittenhouse, and H. R. Sadeghpour, “Ultracold Rydberg molecules,” Nat. Comm. 9, 1965 (2018).
  * Fey  _et al._ (2020) C. Fey, F. Hummel, and P. Schmelcher, “Ultralong-range Rydberg molecules,” Mol. Phys. 118, e1679401 (2020).
  * Greene  _et al._ (2000) C. H. Greene, A. S. Dickinson, and H. R. Sadeghpour, “Creation of polar and nonpolar ultra-long-range Rydberg molecules,” Phys. Rev. Lett. 85, 2458–2461 (2000).
  * Bendkowsky  _et al._ (2009) V. Bendkowsky, B. Butscher, J. Nipper, J. P. Shaffer, R. Löw, and T. Pfau, “Observation of ultralong-range Rydberg molecules,” Nature 458, 1005–1008 (2009).
  * Hamilton  _et al._ (2002) E. L. Hamilton, C. H. Greene, and H. R. Sadeghpour, “Shape-resonance-induced long-range molecular Rydberg states,” J. Phys. B 35, L199–L206 (2002).
  * Bendkowsky  _et al._ (2010) V. Bendkowsky, B. Butscher, J. Nipper, J. B. Balewski, J. P. Shaffer, R. Löw, T. Pfau, W. Li, J. Stanojevic, T. Pohl, and J. M. Rost, “Rydberg trimers and excited dimers bound by internal quantum reflection,” Phys. Rev. Lett. 105, 163201 (2010).
  * Tallant  _et al._ (2012) J. Tallant, S. T. Rittenhouse, D. Booth, H. R. Sadeghpour, and J. P. Shaffer, “Observation of blueshifted ultralong-range Cs2 Rydberg molecules,” Phys. Rev. Lett. 109, 173202 (2012).
  * Bellos  _et al._ (2013) M. A. Bellos, R. Carollo, J. Banerjee, E. E. Eyler, P. L. Gould, and W. C. Stwalley, “Excitation of weakly bound molecules to trilobitelike Rydberg states,” Phys. Rev. Lett. 111, 053001 (2013).
  * Anderson  _et al._ (2014a) D. A. Anderson, S. A. Miller, and G. Raithel, “Photoassociation of long-range $n{D}$ Rydberg molecules,” Phys. Rev. Lett. 112, 163201 (2014a).
  * DeSalvo  _et al._ (2015) B. J. DeSalvo, J. A. Aman, F. B. Dunning, T. C. Killian, H. R. Sadeghpour, S. Yoshida, and J. Burgdörfer, “Ultra-long-range Rydberg molecules in a divalent atomic system,” Phys. Rev. A 92, 031403(R) (2015).
  * Peper and Deiglmayr (2021) M. Peper and J. Deiglmayr, “Heteronuclear long-range Rydberg molecules,” Phys. Rev. Lett. 126, 013001 (2021).
  * Anderson  _et al._ (2014b) D. A. Anderson, S. A. Miller, and G. Raithel, “Angular-momentum couplings in long-range Rb2 Rydberg molecules,” Phys. Rev. A 90, 062518 (2014b).
  * Kleinbach  _et al._ (2017) K. S. Kleinbach, F. Meinert, F. Engel, W. J. Kwon, R. Löw, T. Pfau, and G. Raithel, “Photoassociation of trilobite Rydberg molecules via resonant spin-orbit coupling,” Phys. Rev. Lett. 118, 223001 (2017).
  * Eiles and Greene (2017) M. T. Eiles and C. H. Greene, “Hamiltonian for the inclusion of spin effects in long-range Rydberg molecules,” Phys. Rev. A 95, 042515 (2017).
  * Deiß  _et al._ (2020) M. Deiß, S. Haze, J. Wolf, L. Wang, F. Meinert, C. Fey, F. Hummel, P. Schmelcher, and J. Hecker Denschlag, “Observation of spin-orbit-dependent electron scattering using long-range Rydberg molecules,” Phys. Rev. Research 2, 013047 (2020).
  * Saßmannshausen  _et al._ (2015) H. Saßmannshausen, F. Merkt, and J. Deiglmayr, “Experimental characterization of singlet scattering channels in long-range Rydberg molecules,” Phys. Rev. Lett. 114, 133201 (2015).
  * Böttcher  _et al._ (2016) F. Böttcher, A. Gaj, K. M. Westphal, M. Schlagmüller, K. S. Kleinbach, R. Löw, T. Cubel Liebisch, T. Pfau, and S. Hofferberth, “Observation of mixed singlet-triplet Rb2 Rydberg molecules,” Phys. Rev. A 93, 032512 (2016).
  * Engel  _et al._ (2019) F. Engel, T. Dieterle, F. Hummel, C. Fey, P. Schmelcher, R. Löw, T. Pfau, and F. Meinert, “Precision spectroscopy of negative-ion resonances in ultralong-range Rydberg molecules,” Phys. Rev. Lett. 123, 073003 (2019).
  * MacLennan  _et al._ (2019) J. L. MacLennan, Y.-J. Chen, and G. Raithel, “Deeply bound ($24{D}_{J}+5{S}_{1/2}$) ${}^{87}\mathrm{Rb}$ and ${}^{85}\mathrm{Rb}$ molecules for eight spin couplings,” Phys. Rev. A 99, 033407 (2019).
  * Butscher  _et al._ (2010) B. Butscher, J. Nipper, J. B. Balewski, L. Kukota, V. Bendkowsky, R. Löw, and T. Pfau, “Atom–molecule coherence for ultralong-range Rydberg dimers,” Nat. Phys. 6, 97–974 (2010).
  * Butscher  _et al._ (2011) B. Butscher, V. Bendkowsky, J. Nipper, J. B. Balewski, L. Kukota, R Löw, T. Pfau, W. Li, T. Pohl, and J. M. Rost, “Lifetimes of ultralong-range Rydberg molecules in vibrational ground and excited states,” J. Phys. B 44, 184004 (2011).
  * Camargo  _et al._ (2016) F. Camargo, J. D. Whalen, R. Ding, H. R. Sadeghpour, S. Yoshida, J. Burgdörfer, F. B. Dunning, and T. C. Killian, “Lifetimes of ultra-long-range strontium Rydberg molecules,” Phys. Rev. A 93, 022702 (2016).
  * Khuskivadze  _et al._ (2002) A. A. Khuskivadze, M. I. Chibisov, and I. I. Fabrikant, “Adiabatic energy levels and electric dipole moments of Rydberg states of Rb2 and Cs2 dimers,” Phys. Rev. A 66, 042709 (2002).
  * Li  _et al._ (2011) W. Li, T. Pohl, J. M. Rost, Seth T. Rittenhouse, H. R. Sadeghpour, J. Nipper, B. Butscher, J. B. Balewski, V. Bendkowsky, R. Löw, and T. Pfau, “A homonuclear molecule with a permanent electric dipole moment,” Science 334, 1110–1114 (2011).
  * Bai  _et al._ (2020) S. Bai, X. Han, J. Bai, Y. Jiao, J. Zhao, S. Jia, and G. Raithel, “Cesium $n{D}_{J}+6{S}_{1/2}$ Rydberg molecules and their permanent electric dipole moments,” Phys. Rev. Research 2, 033525 (2020).
  * Booth  _et al._ (2015) D. Booth, S. T. Rittenhouse, J. Yang, H. R. Sadeghpour, and J. P. Shaffer, “Production of trilobite Rydberg molecule dimers with kilo-Debye permanent electric dipole moments,” Science 348, 99–102 (2015).
  * Niederprüm  _et al._ (2016) T. Niederprüm, O. Thomas, T. Eichert, C. Lippe, J. Pérez-Ríos, C. H. Greene, and H. Ott, “Observation of pendular butterfly Rydberg molecules,” Nat. Comm. 7, 12820 (2016).
  * Lesanovsky  _et al._ (2006) I. Lesanovsky, P. Schmelcher, and H. R. Sadeghpour, “Ultra-long-range Rydberg molecules exposed to a magnetic field,” J. Phys. B 39, L69–L76 (2006).
  * Hummel  _et al._ (2019) F. Hummel, C. Fey, and P. Schmelcher, “Alignment of $s$-state Rydberg molecules in magnetic fields,” Phys. Rev. A 99, 023401 (2019).
  * Hummel  _et al._ (2021) F. Hummel, K. Keiler, and P. Schmelcher, “Electric-field-induced wave-packet dynamics and geometrical rearrangement of trilobite Rydberg molecules,” Phys. Rev. A 103, 022827 (2021).
  * Boisseau  _et al._ (2002) C. Boisseau, I. Simbotin, and R. Côté, “Macrodimers: Ultralong range Rydberg molecules,” Phys. Rev. Lett. 88, 133004 (2002).
  * Overstreet  _et al._ (2009) K. R. Overstreet, A. Schwettmann, J. Tallant, D. Booth, and J. P. Shaffer, “Observation of electric-field-induced Cs Rydberg atom macrodimers,” Nat. Phys. 5, 581–585 (2009).
  * Saßmannshausen and Deiglmayr (2016) H. Saßmannshausen and J. Deiglmayr, “Observation of Rydberg-atom macrodimers: Micrometer-sized diatomic molecules,” Phys. Rev. Lett. 117, 083401 (2016).
  * Samboy and Côté (2013) N. Samboy and R. Côté, “Rubidium Rydberg linear macrotrimers,” Phys. Rev. A 87, 032512 (2013).
  * Le Roy (1974) R. J. Le Roy, “Long-range potential coefficients from RKR turning points: ${C}_{6}$ and ${C}_{8}$ for $\mathit{B}(3{\Pi}_{Ou}^{+})$-state ${C}l_{2}$, ${B}r_{2}$, and ${I}_{2}$,” Can. J. Phys. 52, 246–256 (1974).
  * Singer  _et al._ (2005) K. Singer, J. Stanojevic, M. Weidemüller, and R. Côté, “Long-range interactions between alkali Rydberg atom pairs correlated to the ns–ns, np–np and nd–nd asymptotes,” J. Phys. B 38, S295–S307 (2005).
  * Schwettmann  _et al._ (2006) A. Schwettmann, J. Crawford, K. R. Overstreet, and J. P. Shaffer, “Cold Cs Rydberg-gas interactions,” Phys. Rev. A 74, 020701(R) (2006).
  * Deiglmayr  _et al._ (2014) J. Deiglmayr, H. Saßmannshausen, P. Pillet, and F. Merkt, “Observation of dipole-quadrupole interaction in an ultracold gas of Rydberg atoms,” Phys. Rev. Lett. 113, 193001 (2014).
  * Marcassa and Shaffer (2014) L. G. Marcassa and J. P. Shaffer, “Interactions in ultracold Rydberg gases,” Adv. At. Mol. Opt. Phys. 63, 47–133 (2014).
  * Samboy  _et al._ (2011) N. Samboy, J. Stanojevic, and R. Côté, “Formation and properties of Rydberg macrodimers,” Phys. Rev. A 83, 050501(R) (2011).
  * Kiffner  _et al._ (2012) M. Kiffner, H. Park, W. Li, and T. F. Gallagher, “Dipole-dipole-coupled double-Rydberg molecules,” Phys. Rev. A 86, 031401(R) (2012).
  * Han  _et al._ (2018) X. Han, S. Bai, Y. Jiao, L. Hao, Y. Xue, J. Zhao, S. Jia, and G. Raithel, “Cs $62{D}_{J}$ Rydberg-atom macrodimers formed by long-range multipole interaction,” Phys. Rev. A 97, 031403(R) (2018).
  * Han  _et al._ (2019) X. Han, S. Bai, Y. Jiao, G. Raithel, J. Zhao, and S. Jia, “Adiabatic potentials of cesium $n{D}_{2}$ Rydberg-Rydberg macrodimers,” J. Phys. B 52, 135102 (2019).
  * Hollerith  _et al._ (2019) S. Hollerith, J. Zeiher, J. Rui, A. Rubio-Abadal, V. Walther, T. Pohl, D. M. Stamper-Kurn, I. Bloch, and C. Gross, “Quantum gas microscopy of Rydberg macrodimers,” Science 364, 664–667 (2019).
  * Schwettmann  _et al._ (2007) A. Schwettmann, K. R. Overstreet, J. Tallant, and J. P. Shaffer, “Analysis of long-range Cs Rydberg potential wells,” J. Mod. Opt. 54, 2551–2562 (2007).
  * Hollerith  _et al._ (2021) S. Hollerith, J. Rui, A. Rubio-Abadal, K. Srakaew, D. Wei, J. Zeiher, C. Gross, and I. Bloch, “Microscopic electronic structure tomography of Rydberg macrodimers,” Phys. Rev. Research 3, 013252 (2021).
  * Duspayev  _et al._ (2021) A. Duspayev, X. Han, M. A. Viray, L. Ma, J. Zhao, and G. Raithel, “Long-range Rydberg-atom–ion molecules of Rb and Cs,” Phys. Rev. Research 3, 023114 (2021).
  * Deiß  _et al._ (2021) M. Deiß, S. Haze, and J. Hecker Denschlag, “Long-range atom–ion Rydberg molecule: A novel molecular binding mechanism,” Atoms 9 (2021), 10.3390/atoms9020034.
  * Schmid  _et al._ (2010) S. Schmid, A. Härter, and J. H. Denschlag, “Dynamics of a cold trapped ion in a Bose-Einstein condensate,” Phys. Rev. Lett. 105, 133202 (2010).
  * Secker  _et al._ (2016) T. Secker, R. Gerritsma, A. W. Glaetzle, and A. Negretti, “Controlled long-range interactions between Rydberg atoms and ions,” Phys. Rev. A 94, 013420 (2016).
  * Secker  _et al._ (2017) T. Secker, N. Ewald, J. Joger, H. Fürst, T. Feldker, and R. Gerritsma, “Trapped ions in Rydberg-dressed atomic gases,” Phys. Rev. Lett. 118, 263201 (2017).
  * Schmid  _et al._ (2018) T. Schmid, C. Veit, N. Zuber, R. Löw, T. Pfau, M. Tarana, and M. Tomza, “Rydberg molecules for ion-atom scattering in the ultracold regime,” Phys. Rev. Lett. 120, 153401 (2018).
  * Ewald  _et al._ (2019) N. V. Ewald, T. Feldker, H. Hirzler, H. A. Fürst, and R. Gerritsma, “Observation of interactions between trapped ions and ultracold Rydberg atoms,” Phys. Rev. Lett. 122, 253401 (2019).
  * Wang  _et al._ (2020) L. Wang, M. Deiß, G. Raithel, and J. H. Denschlag, “Optical control of atom-ion collisions using a Rydberg state,” J. Phys. B 53, 134005 (2020).
  * Hirzler and Pérez-Ríos (2021) H. Hirzler and J. Pérez-Ríos, “Rydberg atom-ion collisions in cold environments,” Phys. Rev. A 103, 043323 (2021).
  * Dieterle  _et al._ (2021) T. Dieterle, M. Berngruber, C. Hölzl, R. Löw, K. Jachymski, T. Pfau, and F. Meinert, “Transport of a single cold ion immersed in a Bose-Einstein condensate,” Phys. Rev. Lett. 126, 033401 (2021).
  * Born and Huang (1954) M. Born and K. Huang,  _Dynamical theory of crystal lattices_ (Oxford University Press, 1954).
  * Agostini and Curchod (2019) F. Agostini and B. F. E. Curchod, “Different flavors of nonadiabatic molecular dynamics,” Wiley Interdiscip. Rev. Comput. Mol. Sci. 9, e1417 (2019).
  * Born and Oppenheimer (1927) M. Born and R. Oppenheimer, “Zur quantentheorie der molekeln,” Ann. Phys. 389, 457–484 (1927).
  * Koonin and Meredith (1990) S.E. Koonin and D.C. Meredith, “Computational physics Fortran version,” (Addison-Wesley, Reading, MA, 1990) pp. 169–180.
  * Muga  _et al._ (2004) J.G. Muga, J.P. Palao, B. Navarro, and I.L. Egusquiza, “Complex absorbing potentials,” Phys. Rep. 395, 357–426 (2004).
  * Gallagher (2005) T. F. Gallagher,  _Rydberg Atoms_ , Vol. 3 (Cambridge University Press, 2005).
  * Traxler  _et al._ (2013) M. Traxler, R. E. Sapiro, K. Lundquist, E. P. Power, and G. Raithel, “Coupled internal-state and center-of-mass dynamics of Rydberg atoms in a magnetic guide,” Phys. Rev. A 87, 053418 (2013).
  * Anderson  _et al._ (2013) D. A. Anderson, A. Schwarzkopf, R. E. Sapiro, and G. Raithel, “Production and trapping of cold circular Rydberg atoms,” Phys. Rev. A 88, 031401 (2013).
  * Branden  _et al._ (2009) D. B. Branden, T. Juhasz, T. Mahlokozera, C. Vesa, R. O. Wilson, M. Zheng, A. Kortyna, and D. A. Tate, “Radiative lifetime measurements of rubidium Rydberg states,” J. Phys. B 43, 015002 (2009).
  * Beterov  _et al._ (2009) I. I. Beterov, I. I. Ryabtsev, D. B. Tretyakov, and V. M. Entin, “Quasiclassical calculations of blackbody-radiation-induced depopulation rates and effective lifetimes of Rydberg $ns$, $np$, and $nd$ alkali-metal atoms with $n\leq 80$,” Phys. Rev. A 79, 052504 (2009).


Vibrational Dynamics and Spectroscopy of Water at Porous g-C3N4 and C2N
Materials

Deepak Ojha, Christopher Penschke and Peter Saalfrank∗

Theoretische Chemie, Institut für Chemie, Universität Potsdam, Karl-
Liebknecht-Strasse 24-25, D-14476 Potsdam-Golm, Germany

Abstract

Porous graphitic materials containing nitrogen are promising catalysts for
photo(electro)chemical reactions, notably water splitting, but can also serve
as “molecular sieves”. Hearby nitrogen, among other effects, increases the
hydrophilicity of the graphite parent material. A deeper understanding of how
water interacts with C and N containing layered materials, if and which
differences exist between materials with different N content and pore size,
and what the role of water dynamics are – a prerequsite for catalysis and
sieving, is largely absent, however. Vibrational spectroscopy can answer some
of these questions. In this work, the vibrational dynamics and spectroscopy of
deuterated water molecules (D2O) mimicking dense water layers at room
temperature on the surfaces of two different C/N based materials with
different N content and pore size, namely graphitic C3N4 (g-C3N4) and C2N are
studied using Ab Initio Molecular Dynamics (AIMD). In particular, Time-
Dependent vibrational Sum-Frequency Generation spectra (TD-vSFG) of the OD
modes and also time-averaged vSFG spectra and OD frequency distributions are
computed. This allows us to distinguish “free” (dangling) OD bonds from OD
bonds which are bound in a H-bonded water network or at the surface – with
subtle differences between the two surfaces and also to a pure water / air
interface. It is found that the temporal decay of OD modes is very similar on
both surfaces with a correlation time near 4 ps. In contrast, TD-vSFG spectra
reveal that the inter-conversion time from “bonded” to “free” OD bonds being
about 8 ps for water on C2N and thus twice as long as for g-C3N4,
demonstrating a propensity of the former material to stabilize bonded OD
bonds.

∗E-mail<EMAIL_ADDRESS>
##  1 Introduction

Water at interfaces and on surfaces is of immense significance for it plays a
role as a solvent and as a catalyst in many chemical processes1, 2, 3. In
biology, interfacial water at lipid bilayers and proteins contribute to the
overall biological activity and selectivity4. Water on semiconductor surfaces
offers an alternative to conventional energy sources by means of electro- or
photo-catalytic water splitting. Many photo(electro)catalysts contain metals
which are expensive and sometimes hazardous5, 6, 7. In contrast, recently
synthesized organic materials like graphitic carbon nitride (g-C3N4)8, 9, 10
and C2N11, 12 not only have demonstrated their propensity to facilitate water-
splitting, they are also easy to synthesize and non-hazardous to the
environment. These quasi-two dimensional, porous, layered materials offer a
large surface area, and their properties can be tuned by the C/N ratio and the
size of the pores. Nitrogen is particularly important in these materials,
notably for uptake of water, as it considerably increases the hydrophilicity
of the material w.r.t. graphite or graphene12, 13. This makes these materials
also interesting in a non-catalytic context, e.g., as sieves to separate
liquids when acting as water-selective membranes.

Both for catalytic and non-catalytic purposes like selective adsorption and
separation, a precise understanding of the interaction of water with
C/N-containing porous materials is desirable. Several theoretical studies have
been performed which explored the geometric and electronic structure of these
carbon-based materials and how the latter is modified by water adsorption15,
16, 17, 18, 12, 13. Often (but not always) single- or few-layer models were
studied, and low coverages using stationary models. For example, in a first-
principles Density Functional Theory (DFT) study including van-der Waals
corrections and employing a (3$\times$3) periodic supercell model it was
found, that a single water molecule preferentially adsorbs in a triangular
pore of the g-C3N4 surface, with an adsorption energy between -0.55 and -0.62
eV depending on functional18. The water molecule lies almost flat in the pore
of the buckled g-C3N4 surface, forming two hydrogen bonds with accepting N
atoms of the cavity. Also cluster models using highly correlated wavefunction
methods support these results13, 14. (Note that with smaller (1$\times$1) unit
cells the surface is falsely predicted to be flat and other adsorption
geometries are more stable16.) In Ref.18, by a Potential Energy Surface (PES)
scan also other (metastable) adsorption sites for water were found, and
barriers for lateral diffusion of water to these sites were determined – in
the order of 0.2 eV. In Ref.17, in an attempt to go to higher-coverage
situations, multi-layer systems and dynamical behaviour relevant for water-
selective membranes and “molecular sieves”, the diffusion of water (and
ethanol) through pristine and modified g-C3N4 sheets was studied using
classical molecular dynamics simulations based on empirical forcefields. An
inverse correlation between the diffusion coefficient of water and the
lifetime of intermolecular hydrogen bonds (HBs) was found, i.e. short
lifetimes caused larger diffusion coefficients (and vice versa). Diffusion
coefficients (for water) on / between g-C3N4 were also found to be smaller
than for bulk water, and the HB-lifetime was correspondingly larger for water
at / in g-C3N4 compared to bulk water17.

In case of C2N, which consists of flat layers under ideal conditions and with
hexagonal-shaped pores slightly larger than those of g-C3N4, both periodic DFT
and also cluster calculations predict a similar adsorption energy for a single
H2O on a C2N monolayer as for g-C3N4. In this case, however, the water resides
above the (center of the) pore, but still forming two H-bonds to two adjacent,
basic N atoms12, 13. Beyond single-layer models for C2N and after loading with
larger amounts of water, according to the combined experimental-theoretical
work in Ref.12, water molecules reside between the layers, forming H-bonds
with N atoms of C2N which are stronger than the H-bonds between different
water molecules in bulk water. The adsorption energy decreases with increasing
water load.

In order to understand the interaction with and adsorption of liquid water on
2D C/N-containing materials better, it is necessary to consider high
coverages, finite temperatures, and dynamical behaviour which will naturally
arise at finite temperature. (In the studies mentioned above, this was only
done within the classical MD work reported in Ref.17.) Our goal here is to
unravel differences of the C/N material / water systems compared to bulk water
or water / air interfaces, and also differences between various C/N materials.
This is, in the end, a prerequisite to fully understand the performance of
N-containing graphitic materials for water splitting or as sieves, and their
rational design. Another motivation for our work is the mentioned observation
that the H-network but also transport or chemical behaviour of water in
confined spaces / at surfaces can be very different from what is observed in
the bulk (see also 19, 20).

Appropriate “fingerprints” are needed by which differences between various
systems and general effects of a confining surface can be accessed. In this
context, we note that the spectral signature peaks for the concerted
intermolecular motion of hydrogen bonded water molecules as well as the
intramolecular fluctuations of the OH (or OD) modes of water molecules can be
seen in vibrational spectra, such as the infrared (IR) spectrum. In
particular, the IR spectrum provides a fingerprint region to decipher the
hydrogen-bond dynamics in aqueous solutions and at interfaces3. Nevertheless,
conventional vibrational spectroscopy cannot disseminate the contribution of
interfacial molecules from the bulk. Vibrational Sum Frequency Generation
Spectroscopy (vSFG) enables to selectively study the structure, hydrogen-bond
network and orientational profile of interfaces and surfaces21, 22, 23, 24,
25, 26, 27, 28. Moreover, the temporal evolution of the interfacial molecules
can be explored using Time-Dependent vSFG (TD-vSFG)29, 30, 31.

Existing theoretical studies have not explored the vibrational dynamics of
water molecules on two-dimensional, N-containing graphitic materials. In the
present work, we have investigated the vibrational dynamics and vSFG spectra
of interfacial water on g-C3N4 and C2N surfaces using Ab Initio Molecular
Dynamics. We concentrate on surfaces in this work, using single-layer models
for the 2D material, and leaving investigations of the effects of other layers
and the bulk to future work. We will consider thin layers of D2O molecules
(rather than H2O, for practical reasons) adsorbed on one side of the 2D layer,
and analyze OD vibrations in detail. Fluctuations in the vibrational frequency
of OD modes were obtained using the wavelet transform of the time-series
analysis and the surface-specific Velocity Velocity Autocorrelation Function
(ssVVAF) approach32 was used to calculate the vSFG and TD-vSFG spectra of
interfacial water molecules.

The paper is organized as follows. In the following section 2, we briefly
outline theoretical methods and computational models. In Sec.3, we study the
vibrational dynamics, time-averaged vibrational spectroscopy and time-resolved
vibrational spectroscopy, respectively. Sec.4 summarizes and concludes our
work.

##  2 Computational Details and Models

###  2.1 Models

We performed AIMD simulations for D2O molecules adsorbed either on C2N or
g-C3N4 monolayers, using periodic supercell models. For g-C3N4, a supercell of
dimension (20.4 $\times$20.4$\times$24.0) Å with 126 atoms (C54N72) was
employed representing a (3$\times 3$) 2D trigonal elementary cell to model the
wave-like reconstruction found elsewhere 18, 13 (see Fig.1(a) in Ref.13 for a
pictorial representation). We then adsorbed 60 D2O molecules on one side (the
surface) of the monolayer. A snapshot resulting from an AIMD simulation (see
below), is shown in Fig.1(a). As a result, about 14 Å of vacuum separate
different periodic images along $z$ (perpendicular to the surface).

C2N was modelled using a primitive unit cell of dimensions
$(8.32\times{8.32}\times 20.0)$ Å with 18 atoms (composition C12N6),
corresponding to a 2D trigonal elementary cell, see Fig.1(c) in Ref.13). On
one side of this layer, 12 D2O molecules were added, as shown in the form of
an AIMD snapshot in Fig.1(b). This way, a vacuum gap of around 12 Å remained
between water-covered layers.

(a) | (b)
---|---
|
Figure 1: Snapshots of water (D2O) molecules simulated on (a) a C3N4
(3$\times$3) monolayer and (b) a C2N monolayer. The employed unit cells are
shown in side-views, with C: green, N: blue, O: red, D: white.

###  2.2 AIMD simulations

First-principles simulations of water molecules on the g-C3N4 and C2N surfaces
were performed using the Vienna Ab initio Simulation Package33, 34 (VASP,
version 5.4). The inner core electrons were represented by projector
augmented-wave (PAW) pseudopotentials36, 35 and valence electrons using the
Perdew Becke Ernzerhof (PBE) exchange correlation functional37 together with
Grimme’s D3 dispersion correction38, 39. (For a general assessment of DFT for
water, see Ref. 40.) The plane-wave kinetic energy cutoff was set to 400 eV
for both systems. For C2N/water, a 3$\times$3$\times$1 k-point grid was used
for Brillouin zone sampling; for g-C3N4/water, only the $\Gamma$-point was
included.

Simulations were performed in the NVT ensemble using the Nosé-Hoover
thermostat41 at 300 K with deuterium masses for hydrogen and a time step of 1
fs for integrating the equations of motion. The vibrational spectra were
sampled by using single, long NVT trajectories of 45 ps length (C3N4/D2O) and
55 ps (C2N/D2O), respectively. D2O instead of H2O was used in order to allow
for longer timesteps and total propagation times in AIMD, hence better
statistics, and also to inspire experiments which are often done with
deuterated water 42.

###  2.3 Frequency distributions

A (time-dependent) distribution of vibrational frequencies of OD modes of
water molecules on g-C3N4 and C2N surfaces was determined using the wavelet
transform of a time-series analysis. The method is well documented43 and is
based on the principle that a time-dependent function ($f(t)$) can be
expressed in terms of basis functions obtained by the translations and
dilations of a mother wavelet

$\psi_{a,b}(t)=a^{-\frac{1}{2}}\psi\Bigg{(}\frac{t-b}{a}\Bigg{)}\quad,$ (1)

which is represented in the so-called Morlet-Grossman form in our present
study and is mathematically given as,

$\psi(t)=\frac{1}{\sigma\sqrt{2\pi}}e^{2\pi i\lambda
t}e^{-\frac{t^{2}}{2\sigma^{2}}}\quad.$ (2)

The parameters $\lambda$ and $\sigma$ are assigned the values of 1 and 2 sec-1
in our present study. The coefficients of the wavelet expansion are given by
the wavelet transform of $f(t)$, i.e.

$L_{\psi}f(a,b)=a^{-\frac{1}{2}}\int_{-\infty}^{+\infty}f(t){\psi}\Bigg{(}\frac{t-b}{a}\Bigg{)}dt\quad,$
(3)

where $a$ and $b$ are both real quantities. Here, $a>0$ is a scale parameter
which is directly related to the instantaneous frequency content of the system
over a small time-window centered around $b$. Accordingly, based on the
fluctuations in the time-series $f(t)$, the wavelet transform $L_{\psi}f(a,b)$
provides the vibrational frequency for the given small time-window around
$t=b$. Since the frequency value is proportional to $\frac{1}{a}$, the time-
window narrows for high frequency (small $a$) and widens for low frequency
(large $a$). The value of parameter $a$ which maximizes the modulus of the
wavelet transform of time series $f$ at time $t=b$ is used to calculate the
most important frequency component for the given interval.

The time-series of interest in this work, $f(t)$, is constructed as a complex
function with its real and imaginary parts corresponding to fluctuations of
the bond length and momentum of an OD mode projected along the OD bond, namely
(possible unit factors omitted)

$f(t)=\delta r^{\mathrm{OD}}(t)+i\delta p^{\mathrm{OD}}(t)\quad.$ (4)

Here, the fluctuations are defined as

$\delta q(t)=q(t)-\left\langle q\right\rangle$ (5)

where $q(t)$ denotes the instantaneous property ($=r^{\mathrm{OD}}(t)$ or
$p^{\mathrm{OD}}(t)$), and $\left\langle q\right\rangle$ the corresponding
average over the trajectory. Here,
$r^{\mathrm{OD}}(t)=|\underline{r}^{\mathrm{OD}}(t)|$ where
$\underline{r}^{\mathrm{OD}}(t)=\underline{r}_{\mathrm{O}}-\underline{r}_{\mathrm{D}}$
is the OD bond vector for a given OD bond, obtained from atomic position
vectors. Further, the momentum projected along the OD mode is given as

$p^{\mathrm{OD}}(t)=\left(\frac{m_{\mathrm{D}}\cdot
m_{\mathrm{O}}}{m_{\mathrm{D}}+m_{\mathrm{O}}}\right)\cdot(\underline{v}_{\mathrm{O}}-\underline{v}_{\mathrm{D}})\cdot\underline{\hat{r}}^{\mathrm{OD}}\quad,$
(6)

where $m_{\mathrm{D}},m_{\mathrm{O}}$ are atomic masses of deuterium and
oxygen, ${\underline{v}}_{\mathrm{O}}$, ${\underline{v}}_{\mathrm{D}}$ are
atom velocities, and
${\underline{\hat{r}}}^{\mathrm{OD}}={\underline{r}}^{\mathrm{OD}}/r^{\mathrm{OD}}$
the unit vector along the OD bond. This method is then applied to all the OD
modes present in a given system. For further details and recent applications
of the method within AIMD, see Refs.43, 44, 45, 46.

##  3 Results

###  3.1 Vibrational frequency distribution and dynamics

####  3.1.1 Time-averaged frequency distribution

We first look at time-averaged vibrational frequency distributions $P(\omega)$
of OD modes, obtained from time-averaging the time-dependent vibrational
frequencies of OD modes resulting from the wavelet transform of $f(t)$.
Practically, an averaged vibrational frequency distribution is obtained by
dividing the frequency range of 2000-3000 cm-1 in equally spaced bins. If at a
given instant the frequency of an OD mode falls within a given bin, the bin
height is incremented by unity. This process is repeated for each OD mode and
along the entire trajectory.

Figure 2:  Time-averaged frequency distribution $P(\omega)$ (in arbitrary
units) of the OD modes of water molecules on the (a) g-C3N4 and (b) C2N
surface.

The time-averaged vibrational frequency distribution of OD modes of the
interfacial water on the g-C3N4 surface was found to be a distribution ranging
from about 2000 cm-1 to about 2900 cm-1 with a mean at 2496 cm-1 as shown in
Fig.2(a). The distribution is similar to a distribution found for bulk water
at 300 K using a similar methodology (see Fig.2(a) of Ref.46, for example),
however, with two differences. First, the distribution in our Fig.2(a) is
slightly broader, and, more importantly, a clear shoulder peak now appears in
the high frequency region around 2750 cm-1 which is absent in bulk water.
Closer analysis shows that this shoulder corresponds to the free or dangling
OD modes, i.e., OD bonds which are not part of a H-bond network and not bound
to the surface and stick out into vacuum instead. Similarly, the average
vibrational frequency of OD modes on the C2N surface is 2487 cm-1 (Fig.2(b)),
and the whole time-averaged distribution is similar to the other surface. Also
for C2N, a shoulder peak in the high-frequency region around 2750 cm-1 is
found, again mainly due to the free/dangling OD modes of the interfacial
water. Due to the non-centrosymmetric environment at the surface, the
tetrahedral hydrogen-bond network is distorted for both surfaces and thus the
population of free OD modes shows an overall increment. The high-frequency
shoulder is more prominently observable for g-C3N4 which implies comparatively
higher population of free/dangling OD modes. Interestingly, high-frequency
shoulders are also a characteristic feature of vSFG (vibrational Sum Frequency
Generation) spectra of other systems with interfacial water24, 25, 26, 47, 42.

####  3.1.2 Time-resolved frequency distribution and dynamics

Further, the temporal dynamics of the OD modes of water molecules on the
surface of g-C3N4 and C2N were studied using time-dependent joint probability
distributions. The time-dependent joint probability distribution used can be
expressed mathematically as

$P(\omega_{3},t_{2},\omega_{1})=\left\langle\delta\left(\omega(t_{2})-\omega_{3}\right)\cdot\delta\left(\omega(0)-\omega_{1}\right)\right\rangle\quad.$
(7)

Here, $\delta(a(t)-b)=a(t)-b$ in analogy to Eq.(5), and the average is taken
over all OD bonds. The joint probability distribution gives the probability
that a given OD mode which was oscillating with the vibrational frequency
$\omega_{1}$ evolves to a frequency $\omega_{3}$ within a time interval
$t_{2}$. (Note that for centrosymmetric systems other than those studied here,
$P(\omega_{3},t_{2},\omega_{1})$ is similar to a two-dimensional IR (2D-IR)
spectrum at long waiting times $t_{2}$ 46.)

Figure 3:  Joint probability frequency distribution
$P(\omega_{3},t_{2},\omega_{1})$ of OD modes of water molecules on the g-C3N4
surface for waiting times $t_{2}=10,100,750$ and $3500$ fs.

In Fig.3, we have plotted the distributions for the waiting times
$t_{2}=10,100,750$ and $3500$ fs for OD modes of water molecules on the g-C3N4
surface respectively. For small waiting times, i.e., $t_{2}=10$ fs, the
distribution is predominantly a straight line aligned along the diagonal
$\omega_{1}=\omega_{3}$, indicating that the OD modes keep close to their
initial vibrational frequency. However, for the longer waiting times, the
distribution eventually gets elongated and evolves into a spherical
distribution which implies that the OD modes undergo vibrational spectral
diffusion due to hydrogen bond rearrangement. Moreover, within 3.5 ps the OD
modes have sampled all accessible vibrational states / frequencies (cf.
Fig.2(a)). A key observation in the time-resolved frequency distributions is
that the high frequency shoulder is evidently observable at low as well as
long waiting times. Further, the vibrational dephasing of the modes
corresponding to the shoulder peak is seen progressing as vertical slab along
2750 cm-1.

A similar analysis for the water molecules on the C2N surface corresponding to
the waiting times of $t_{2}=10,100,750$ and $3500$ fs is shown in Fig.4. The
vibrational dephasing follows a similar pattern, and overall differences
between both surfaces seem small.

Figure 4: Joint probability frequency distribution
$P(\omega_{3},t_{2},\omega_{1})$ of OD modes of water molecules on the C2N
surface for waiting times $t_{2}=0,100,750$ and $3500$ fs.

To more quantitatively address this point, the temporal dynamics of the OD
modes of water molecules on the g-C3N4 and C2N surfaces is also analyzed using
a frequency time correlation function (FTCF) which is given as,

$C_{\omega\omega}(t)=\frac{\left\langle\delta\omega(0)\cdot\delta\omega(t)\right\rangle}{\left<\delta\omega\left(0\right)^{2}\right>}\quad.$
(8)

The time-dependent decay of the correlation function of OD modes of water
molecules on the g-C3N4 and C2N surfaces is shown in Fig.5(a).

Figure 5:  (a) Frequency correlation function $C_{\omega\omega}(t)$ and (b)
hydrogen bond correlation function $S_{HB}(t)$ for OD modes of water molecules
on the g-C3N4 and C2N surfaces.

Evidently, the decay is biphasic with short time regime which extends up to
100 fs and is governed by the ultrafast librational dynamics or concerted
collective motion of water molecules48 whereas the long time decay which
extends up to a few picoseconds and is mainly attributed to the hydrogen bond
rearrangement dynamics. The timescale of vibrational dephasing was determined
by using a bi-expontential fitting function,

$C_{\omega\omega}(t)=a_{0}\exp\left(-\frac{t}{\tau_{0}}\right)+(1-a_{0})\exp\left(-\frac{t}{\tau_{1}}\right)\quad.$
(9)

The short time constant $\tau_{0}$ was found to be 76 fs and the long time
decay constant $\tau_{1}$ was 5.7 ps for water on g-C3N4. For C2N, we find
$\tau_{0}=104$ fs and $\tau_{1}=4.4$ ps. We note that for both systems the
long time constant, $\tau_{1}$, is significantly longer than for bulk liquid
water49, 46. According to Ref. 46, a computed $\tau_{1}$ for bulk D2O is 1.6
ps. Thus, vibrational dynamics of OD modes are slowed down in proximity of a
the g-C3N4 surface, and slightly less so near C2N.

Since frequency fluctuations of the OD modes are associated with the
rearrangement of the local hydrogen bond network, we have also calculated a
“hydrogen bond lifetime” using a continuous HB correlation function given as,

$S_{HB}(t)=\frac{\left<h(0)H(t)\right>}{\left<h(0)^{2}\right>}\quad.$ (10)

Here, $h(t)$ and $H(t)$ are hydrogen bond population labels, defined as
follows. For a given pair of hydrogen-bonded, different water molecules at an
instant $t$, $h(t)$ is assigned a value of 1 and otherwise 0. (Time zero
refers to the beginning of the segment of the trajectory into which the
trajectory was divided.) Further, if the hydrogen bond between a pair of water
molecules remains intact from time $t=0$ till time $t$, $H(t)$ is assigned a
value of 1 or else 0. The water molecules are considered to be hydrogen-
bonded, if the intermolecular O $\dots$ D distance is less than 2.45 Å
corresponding to the first minimum in the O $\cdots$ D pair correlation
function. The temporal decay of the hydrogen bond correlation function
indicates the timescale of rearrangement of the H-bond network and is shown
for water (D2O) on g-C3N4 and C2N surfaces in Fig.5(b). Evidently, also within
this measure, the temporal decay of the $S_{HB}(t)$ correlation function is
somewhat faster for C2N/water in comparison to g-C3N4/water, as can be seen by
visual inspection of Fig.5(b). More quantitatively, the long hydrogen bond
lifetime $\tau_{1}$, obtained from bi-exponential fits of $S_{HB}(t)$, are
found to be $\tau_{1}=1.8$ ps for g-C3N4 and 1.3 ps for water on the C2N
surface, respectively. Calculations of $S_{HB}(t)$ have been done also for
bulk deuterated water and for the water-vacuum interface in Ref.50. The H-bond
lifetimes computed there are of similar length, but hard to compare
quantitatively to our numbers since they were obtained with another DFT
functional and with Car-Parrinello MD rather than Born-Oppenheimer MD. Still,
in Ref. 50 it was found that the bulk $S_{HB}$ function decays faster than the
interface one. Here, we have the additional effect of the porous surface and
we see differences between the two surfaces, g-C3N4 and C2N.

###  3.2 Vibrational Sum Frequency Generation (vSFG) spectra

Further insight into the water dynamics on the studied 2D materials, more
directly accessible by experiment, comes from vibrational spectroscopy. Here
we concentrate on vibrational SFG as a particularly surface-sensitive method.

####  3.2.1 Time-averaged vSFG spectra

We start with time-averaged vSFG spectra. A time-averaged vSFG spectrum
provides relevant frequencies and gives information about the orientation
profile of molecular dipoles at the surface/interface. The time-averaged vSFG
spectrum is determined from the second-order optical susceptibility tensor of
the system, with tensor elements27, 28

${\chi^{2}_{abc}}\propto\int e^{-i\omega
t}\left\langle{{{\alpha}_{ab}}(t){{{\mu}_{c}}}(0)}\right\rangle dt\quad.$ (11)

Here, $a$, $b$, and $c$ are Cartesian coordinates, $\alpha_{ab}(t)$ is the
($a$,$b$)-component of the polarizability tensor at time $t$, $\mu_{c}(0)$ the
$c$-component of the dipole vector at time $t=0$, and $\omega$ the frequency.
Squared, summed (and weighted) components of ${\chi^{2}_{abc}}$ give the vSFG
spectrum. Here, we concentrate on so-called ssp-polarization spectra, which we
approximate by only considering a single ${\chi^{2}_{abc}}$ element and
setting $a=b=x$ and $c=z$, respectively, where $x$ is a direction parallel to
the surface and $z$ perpendicular to it. We show real and imaginary parts of
the second-order susceptibility, ${\chi^{2}_{xxz}}(\omega)$, whereby the
imaginary part can be attributed to the orientation of the interfacial
molecules with respect to the surface21, 22, 23. In practice, we calculate
here the ssVVAF method mentioned earlier32, 47. Then, Eq.(11) can be
rewritten, for our D2O / surface vSFG problem (with only O-D stretch
vibrations considered), as32, 47

${\chi^{2}_{abc}}(\omega)\propto\begin{dcases}\dfrac{\mu^{\prime}_{\mathrm{OD}}\alpha^{\prime}_{\mathrm{OD}}}{i\omega^{2}}\int^{\infty}_{0}dt\
e^{-i\omega
t}\left\langle\sum_{i,j}^{M}{{v}}^{\text{OD}}_{c,i}(0)\dfrac{{\underline{v}}^{{\text{OD}}}_{j}(t)\cdot{{\underline{r}}}^{{\text{OD}}}_{j}(t)}{|{\underline{r}}^{\text{OD}}_{j}(t)|}\right\rangle,&a=b\\\
0,&a\neq b\quad.\end{dcases}$ (12)

Here, $\underline{r}^{{\text{OD}}}_{i}$ is the coordinate vector of an OD bond
($i$) as defined earlier (there are $M=2\cdot N$ O-D bonds for $N$ water
molecules),
$\underline{v}^{{\text{OD}}}_{i}=d\underline{r}^{{\text{OD}}}_{i}/dt$ the
corresponding velocity, ${{v}}^{\text{OD}}_{c,i}$ its $c$-component, and
$\mu^{\prime}_{\mathrm{OD}}$ and $\alpha^{\prime}_{\mathrm{OD}}$ are
derivatives of the dipole moments and polarizabilities, respectively. A few
things are worth mentioning: The ssVVAF expression (12) is approximate in the
sense that only O-D stretch vibrations are considered. Further, we also use a
parameterized (linear) version of dipole moments and polarizabilities without
explicitly calculating them – in fact we also do also not include the
frequency dependence of $\mu^{\prime}_{\mathrm{OD}}$ and
$\alpha^{\prime}_{\mathrm{OD}}$ and set them as constants. Finally, a
sometimes used “quantum correction factor” $Q(\omega)$32 is also neglected.
Still, the approach followed here will give insight into the O-D dynamics and
in particular into the dynamics of the interfacial H-bond network. The method
has been shown to give for the O-H stretch region of a water / air interface,
quite similar vSFG spectra compared to a non-parameterized procedure with
calculated dipole moments and polarizabilities28.

Figure 6:  Real and imaginary parts of $\chi^{2}_{xxz}$ for D2O at (a) g-C3N4
and (b) C2N surfaces. The y-axis is in arbitrary units, with the largest
positive intensity of imaginary or real parts normalized to 1.

In Fig.6, we show the real and imaginary parts of $\chi^{2}_{xxz}$ for D2O at
g-C3N4 (a) and C2N, respectively, again in the frequency range of the OD
stretch vibrations. To gain information on the orientation of OD bonds, we
focus on the imaginary part of $\chi^{2}_{xxz}$. (The real part of the vSFG
spectrum is similar to a conventional infrared spectrum and provides
information about the frequency required to excite vibrations of the adsorbed
molecules.) For D2O/g-C3N4, the imaginary part of the susceptibility is
characterized by a broad, (mostly) negative peak for the frequency domain of
2400-2700 cm-1 which refers to the OD modes which are hydrogen-bonded with the
other water molecules or with the g-C3N4 surface (Fig.6(a)). We also note a
sharp positive peak for the frequency domain of 2700-2800 cm-1, referring to
free, dangling, or weakly bound OD units (with D2O not having four D-bonded
neighbours), which are orientated away from the surface. For example, O-D
bonds pointing into vacuum are of this type. Note that positive intensities
refer to dominantly “upward” and negative intensities dominantly to “downward”
orientations of OD bonds. We also note that the intensity of this peak remains
unchanged relative to the (unsupported) pure water/air interface computed in
Ref.30 (Fig.2 there), while the intensity of bonded, red-shifted OD modes is
relatively reduced as compared to this reference. This may imply an overall
reduction of contributing bonded OD bonds, relative to the free/non-bonded
ones, in case of water on C/N surfaces compared to the water/air interface.
This is in agreement with the frequency distribution of water molecules on
g-C3N4 as shown in Fig.1(a), where the population of free/dangling OD modes is
clearly visible as a high-frequency shoulder.

Further, in Fig.6(b) the second-order susceptibility of the water molecules on
the C2N surface is shown. The imaginary part of the second order
susceptibility for this system is more similar to the pure water-air
interface30. Namely, similar to the latter, the negative-intensity peak
covering the spectral domain of 2200-2700 cm-1 for the water/C2N system is
broader and of larger intensity (comparable to the high-frequency peak)
compared to the water/g-C3N4 surface. (Note that in Refs.28, 30 non-deuterated
water was considered so absolute oxygen-hydrogen frequencies are different.)

To summarize, the vSFG spectra of OD modes on the g-C3N4 and C2N surfaces have
similar features for the high frequency region around 2700-2800 cm-1. But for
the frequency region of 2200-2700 cm-1, the spectral intensity is depleted in
the case of g-C3N4, while it is not for water/C2N which behaves similar to the
water-air interface28. We attribute this behaviour to a propensity of the C2N
material to stabilize bonded OD bonds.

####  3.2.2 Time-resolved vSFG spectra

Time-averaged vSFG spectra provide information about the orientation of
molecular dipoles at the interface. Nevertheless, dynamics of these molecular
dipoles cannot be studied with the conventional vSFG technique. Here, time-
resolved vSFG spectra offer more insight. A time-dependent vibrational SFG
spectrum can be obtained by modifying Eq.(11) as

${\chi^{2}_{abc}}\lvert_{\omega({t}^{\prime})={\omega}^{\prime}}\hskip
11.38092pt=\int e^{-i\omega
t}\left\langle{{\alpha}_{ab}}({t}^{\prime}+T_{w})\cdot{{\mu}_{c}}({t}^{\prime}+T_{w}+t)\right\rangle
dt\quad.$ (13)

Here, we compute the second-order susceptibility of a given vibrational
oscillator provided that the oscillator was vibrating at a frequency
${\omega}^{\prime}$ at a time instant ${t}^{\prime}$. By systematically
varying the value of the waiting time, $T_{w}$ in Eq.(13), the delay time
between a pump pulse and the overlapping IR/vis pulses which generate an SFG
signal, we can obtain frequency-resolved, time-dependent susceptibilities and
from there, TD-vSFG spectra30. It is important to note that in the given
computational framework, a priori information of the instantaneous frequency
of the oscillators/OD modes is needed to obtain the TD-vSFG spectrum. This
information is obtained from the wavelet transform as described above.
Further, Eq.(13) can easily be rewritten in the ss-VVAF form (12), which we
use in the following to compute Im $\chi^{(2)}(\omega;T_{w})$ (actually, Im
$\chi^{(2)}(\omega_{IR};T_{w},\omega_{pump})$).

Figure 7:  Imaginary part of $\chi^{(2)}(\omega)$ at different delay times
$T_{w}$, computed for water on g-C3N4 for (a) the H-bonded OD modes and (b)
the free/dangling OD modes. For C2N, the same information is given in panels
(c) for bonded OD modes and (d) for free/dangling OD modes respectively. Now,
for y-axes arbitrary but non-normalized units are used to make the time-
evolution of signal intensities visible.

To calculate the TD-vSFG spectra of the hydrogen bonded OD modes, we have
applied a broadband IR pump pulse corresponding to the frequency domain of
2200-2700 cm-1. Technically, this was done by computing the susceptibilities
of the water molecules based on their vibrational excitation frequency within
this frequency range. The imaginary parts of the second-order susceptibility
of the hydrogen-bonded OD modes of the interfacial water molecules on g-C3N4
corresponding to the waiting times $T_{w}=0,400,2000,3000$ and $4000$ fs are
shown in Fig.7(a). The spectrum for the waiting time $T_{w}=0$ fs, the
“impulsive limit”, is dominated by a peak with negative intensity centered
around 2600 cm-1. With increasing waiting times, e.g. for $T_{w}=$ 400 fs, we
see a formerly shallow shoulder in the high frequency region around 2750 cm-1
growing, corresponding to the free / dangling OD modes. Further, within 4.0
ps, the peak has attained the maximum intensity and becomes saturated. The
figure demonstrates the inter-conversion from hydrogen/surface-bonded
(“bound”) OD modes to free OD modes, taking place on this timescale.

Similarly, we have calculated spectra of free / dangling OD modes of g-C3N4
for waiting times $T_{w}=0,200,1000$ and 2000 fs in Fig.7(b). Here we applied
a narrowband IR pump pulse in the range of 2700-2900 cm-1. The spectrum for
$T_{w}=0$ fs is a single, positive peak centered around 2750 cm-1. With
increasing waiting times $T_{w}=200,1000,2000$ fs, the frequency domain
between 2200 and 2700 cm-1 corresponding to the H-bonded region shows finite-
valued, negative intensity, while at the same time the high-frequency,
positive peak looses intensity. This indicates interconversion of free OD
modes to (H-)bound OD modes. Closer inspection shows, as can also be seen from
the shorter timescale shown in Fig.7(b), that the free $\rightarrow$ bound
interconversion proceeds faster (within about 2 ps) in (b) than the bound
$\rightarrow$ free interconversion (about 4 ps) in (a). This is due to the
free OD modes having a strong propensity to form hydrogen bonds and stabilize.
When compared to the pure water / air interface30, 31, it is found that the
hydrogen-bond inter-conversion dynamics for the bonded OD modes are overall
slowed down whereas for the free OD, the timescale remains relatively
unchanged in the presence of the g-C3N4 surface.

Along the same lines, we have also obtained the spectra of the bonded and
dangling/free OD bonds of water on the C2N surface as shown in Fig.7(c) and
(d). The spectra of the bonded OD modes are shown for waiting times of
$T_{w}=$0, 400, 2000, 5000 and 8000 fs, respectively. While the spectral
features as seen in the TD-vSFG spectra are similar to previous calculations
for the water / air interface30, a key difference is a very slow temporal
evolution of the peak for the free OD modes which extends up to nearly 8 ps.

The spectra of dangling/free OD modes are shown for the waiting times $T_{w}=$
0, 200, 1000, 2000 and 5000 fs respectively. Again, the free OD modes have a
strong propensity to form hydrogen bonds due to which the spectral domain
corresponding to the bonded OD modes proceeds mostly within 2 ps. The spectra
of free OD modes of interfacial water molecules are interestingly independent
of the surface, i.e., g-C3N4 or C2N.

To summarize, we note that although the overall vibrational dynamics of OD
modes on g-C3N4 and C2N are similar according to the frequency correlation
function, Fig.5(a) (with a tendency of faster H-bond dynamics for C2N), the
inter-conversion rates of free and bonded OD modes can be still significantly
different (and of other order). This refers in particular to the bound
$\rightarrow$ free interconversion of OD bonds, where D2O/C2N exhibits a
significantly longer ($\sim 8$ ps) timescale than D2O/g-C3N4 ($\sim 4$ ps).
Here we want to emphasize that the time-resolved vSFG spectra of hydrogen-
bonded OD modes on the g-C3N4/C2N gives a comparative perspective of
stabilization of the OD modes on two surfaces. The observed differences in the
inter-conversion rates for the two systems imply that g-C3N4 could be more
suitable for catalyzing interfacial reactions involving water, e.g., water
dissociation because the OD modes have a relatively higher propensity to break
the hydrogen bonds ($\sim 4$ ps). In contrast, C2N may be more suited for
processes like desalination or sieving which require stabilization of water by
the surface with an average inter-conversion rate of nearly 8 ps .

The key differences between both systems is that the g-C3N4 has a buckled
surface whereas C2N is flat. Moreover, from the stoichiometric perspective
g-C3N4 has a higher percentage of nitrogen as compared to C2N. These factors
obviously affect the HB dynamics of adsorbed water. Finally we note that the
present inferences are drawn based on the TD-vSFG spectra of water molecules
on a single monolayer of g-C3N4 or C2N. The role of a multi-layered surface as
well as density of water molecules corresponding to liquid phase in the
simulation, may bring more insight into the catalytic properties of the two
materials. Work along these lines is in progress in our laboratory.

##  4 Summary

To summarize, we have studied the vibrational dynamics and (time-dependent)
vSFG spectra of deuterated water molecules at g-C3N4 and C2N surfaces
respectively, using DFT based AIMD simulations. The most important results are
as follows.

  * •

The time-averaged frequency distribution for the OD modes of water molecules
on both surfaces shows a high-frequency peak corresponding to free OD modes.
The “bound” OD bonds, either engaged in the intermolecular H-bond network or
interacting with the C/N surface, are usually weakened and therefore red-
shifted w.r.t. the high-frequency peak.

  * •

The time-dependent decay of the OD frequency (measured by a frequency
correlation function) is quite similar for both surfaces with a timescale of
about 4 ps which is considerably slower than that of the bulk liquid
(deuterated) water. Similar trends for the two studied C/N surfaces arise when
analyzing the H-bond network via a H-bond measure $S_{HB}(t)$, which implies a
slightly longer timescale for g-C3N4 than for C2N.

  * •

Time-independent vSFG spectra (measured by $\chi^{2}_{xxz}$), indicate that
the dangling OD bonds, pointing upward away from the water / vacuum interface
are dominant features in the spectra and are characterized by positive
imaginary parts of $\chi^{2}_{xxz}$ at high frequencies. “Bound” OD modes are
usually of lower intensity, have negative imaginary parts of $\chi^{2}$ and
are red-shifted. As a difference between the two surfaces, the spectral domain
around 2000-2600 cm-1 of the vSFG spectrum for D2O/C2N is marked by a stronger
negative intensity peak compared to g-C3N4.

  * •

Detailed insight into OD bond dynamics is provided by TD-vSFG spectra, with
interesting trends for the two surfaces. Although the vibrational frequency
correlation functions are similar, the time-resolved vSFG of bonded OD modes
saturates within about 4 ps for g-C3N4 whereas the same process takes nearly 8
ps for the C2N surface. The “free” $\rightarrow$ “bound” interconversion is
faster ($\sim 2$ ps) and more similar for both surfaces. Thus we can infer on
the basis of slower interconversion rates for the bonded to free OD state that
the C2N surface has a propensity to stabilize bonded OD moieties.

All of these observations imply different performance of porous C/N-containing
materials in either catalysis or sieving. That key factors like pore size, C/N
ratio, or the three-dimensional structure of the materials (e.g., flat or
buckled) influence the vibrational dynamics of adsorbed water has been
demonstrated, however, how these factors act in detail is a matter of further
research.

## Acknowledgment

P.S. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG),
through project Sa 547-18, and Ch.P. by Deutsche Forschungsgemeinschaft within
Germany’s Excellence Strategy – EXC 2008/1-390540038, UniSysCat. The authors
thank the Paderborn Center for Parallel Computing (PC2) for computation time
on supercomputer “Noctua 2”.

## References

  * 1 O. Björneholm, M.H. Hansen, A. Hodgson, L.-M. Liu, D.T. Limmer, A. Michaelides, P. Pedevilla, J. Rossmeisl, H. Shen, G. Tocci, E. Tyrode, M.-M. Walz, J. Werner, and H. Bluhm,  Chem. Rev. 2016, 116, 7698-7726.
  * 2 Y. Jung and R. A. Marcus  J. Amer. Chem. Soc. 2007, 129, 5492–5502.
  * 3 J. L. Skinner, P. A. Pieniazek, and S. M. Gruenbaum,  Acc. Chem. Res. 2012, 45, 93-100.
  * 4 P. Jungwirth,  J. Chem. Phys. Lett. 2015, 6, 2449–2451.
  * 5 G. Gonella et al.,  Nat. Rev. Chem. 2021, 5, 466–485.
  * 6 J. Zhu et al.  Angew. Chem., Int. Ed. 2015, 54, 9111-9114.
  * 7 N. G. Hörrmann, Z. Guo, F. Ambrosio, O. Andreussi, A. Pasquarello, and N. Marzari,  npj Comput Mater 2019, 5, 100-105.
  * 8 X. Wang, K. Maeda, A. Thomas, K. Takanabe, G. Xin, J. M. Carlsson, K. Domen, and M. Antonietti ,  Nat. Mater. 2009, 8 , 76-80.
  * 9 Y. Wang, X. Wang, and M. Antonietti,  Angew. Chem. Int. Ed 2009 , 51, 68-89.
  * 10 X. Wang , S. Blechert and M. Antonietti,  ACS Catal. 2012, 2, 1596-1606, .
  * 11 R. Walczak, B. Kurpil, A. Savateev, T. Heil, J. Schmidt, Q. Qin, M. Antonietti, and M. Oschatz,  Angew. Chem. Int. Ed. 2018, 57, 10765-10770.
  * 12 J. Heske, R. Walczak, J.D. Epping, S. Youk, S.K. Sahoo, M. Antonietti, T.D. Kühne, M. Oschatz,  J. Mat. Chem. A 2021, 9, 22563-22572.
  * 13 Ch. Penschke, A. Zehle, N. Jahn, R.E. von Zander, R. Neumann, A. Beqiraj, and P. Saalfrank,  Phys. Chem. Chem. Phys. 2022, 24, 14709-14726.
  * 14 T. Schäfer, A. Gallo, A. Irmler, F. Hummel, and A. Grüneis,  J. Chem. Phys. 2021, 155, 244103.
  * 15 H.-Z. Wu, L.-M. Liu, S.-J. Zhao  Phys. Chem. Chem. Phys., 2014, 16, 3299-3304.
  * 16 S. M. Aspera1, M. David, and H. Kasai  Jpn. J. Appl. Phys 2010, 49, 115703.
  * 17 X. Zou, M. Li, S. Zhou, Ch. Chen, J. Zhong, A. Xue, Y. Zhang, and Y. Zhao,  J. Membr. Sci. 2019, 585, 81-89.
  * 18 J. Wirth, R. Neumann, M. Antonietti and P. Saalfrank,  Phys. Chem. Chem. Phys 2014, 16, 15917-15926.
  * 19 I.-C. Yeh and G. Hummer,  J. Phys. Chem. B 2004, 108, 15973-15879.
  * 20 N. Kavokine, R. R. Netz, and L. Bocquet, Annu. Rev. Fluid Mech. 2021, 53, 377–410.
  * 21 Q. Du, R. Superfine, E. Freysz, and Y. R. Shen,  Phys. Rev. Lett., 1993, 70, 2313–2316.
  * 22 Y. R. Shen, and V. Ostroverkhov,  Chem. Rev., 2006, 106, 1140–1154.
  * 23 V. Ostroverkhov, G. A. Waychunas, Y. R. Shen,  Phys. Rev. Lett., 2005, 94, 046102.
  * 24 A. Morita and J. T. Hynes,  Chem. Phys. 2000 258, 371–390.
  * 25 A. Morita and J. T. Hynes,  J. Phys. Chem. B, 2002, 106, 673–685.
  * 26 T. Ishiyama, T. Imamura and A. Morita,  Chem. Rev., 2014, 114, 8447–8470.
  * 27 A. Morita, Theory of Sum Frequency Generation Spectroscopy (Springer, Singapore 2018)
  * 28 D. Ojha and T. D. Kühne,  Molecules, 2020, 25, 3939.
  * 29 E. H. G. Backus, J. D. Cyran, M. Grechko, Y. Nagata, and M. Bonn,  J. Phys. Chem. A 2018, 122, 2401–2410.
  * 30 D. Ojha, N. K. Kaliannan, T. D. Kühne  Commun. Chem., 2019, 2, 116.
  * 31 D. Ojha and T. D. Kühne,  Sci. Rep., 2021, 11, 2456.
  * 32 T. Ohto, K. Usui, T. Hasegawa, M. Bonn, and Y. Nagata,  J. Chem. Phys., 2015, 143, 124702.
  * 33 G. Kresse and J. Furthmüller,  Comput. Mater. Sci. 1996, 6, 15–50.
  * 34 G. Kresse and J. Furthmüller,  Phys. Rev. B 1996, 54, 11169–11186.
  * 35 P. E. Blöchl,  Phys. Rev. B 1994, 50, 17953–17979.
  * 36 G. Kresse and D. Joubert,  Phys. Rev. B 1999, 59, 1758–1775.
  * 37 J. P. Perdew, K. Burke and M. Ernzerhof,  Phys. Rev. Lett. 1996, 77, 3865–3868.
  * 38 S. Grimme, J. Antony, S. Ehrlich and H. Krieg,  J. Chem. Phys. 2010, 132, 154104.
  * 39 S. Grimme, S. Ehrlich and L. Goerigk,  J. Comput. Chem. 2011, 32, 1456–1465.
  * 40 M.J. Gillan, D. Alfé, and A. Michaelides,  J. Chem. Phys. 2016, 144, 13091.
  * 41 S. Nosé,  J. Chem. Phys. 1984, 81, 511–519.
  * 42 G. Melani, Y. Nagata, R.K. Campen, and P. Saalfrank,  J. Chem. Phys., 2019, 150, 244701.
  * 43 A. Semparithi and S. Keshavamurthy,  Phys. Chem. Chem. Phys., 2003, 5, 5051–5062.
  * 44 D. Ojha, A. Henao, and T. D. Kühne,  J. Chem. Phys., 2018, 148, 102328.
  * 45 D. Ojha, K. Karhan, and T. D. Kühne,  Sci. Rep., 2018, 8, 16888.
  * 46 D. Ojha and T.D. Kühne,  Phys. Chem. Chem. Phys., 2023, 25, 13442.
  * 47 G. Melani, Y. Nagata, J. Wirth, and P. Saalfrank,  J. Chem. Phys., 2018, 149, 014707.
  * 48 M. Heyden, J. Sun, S. Funkner, G. Mathias, H. Forbert, M. Havenith, and D. Marx Proc. Nat. Acad. Sci., 2010, 107, 12068–12073.
  * 49 D. Ojha and A. Chandra, Chem. Phys. Lett. 2020, 751, 137493.
  * 50 J.R. Choudhuri and A. Chandra, J. Chem. Phys. 2014, 141, 194705.