The nuclear-spin-forbidden rovibrational transitions of water from first principles
Andrey Yachmenev
Guang Yang
Emil Zak
Sergei Yurchenko
Jochen Küpper

The water molecule occurs in two nuclear-spin isomers that differ by the value of the total
nuclear spin of the hydrogen atoms, , $I=0$ for para-and $I=1$ for
ortho-. Spectroscopic transitions between rovibrational states of ortho and
para water are extremely weak due to the tiny hyperfine nuclear-spin-rotation interaction
of only $\ordsim30$ kHz and so far were not observed. We report the first comprehensive
theoretical investigation of the hyperfine effects and transitions in
H$_2$$^{16}$O due to nuclear-spin-rotation and spin-spin interactions. We also present the
details of our newly developed general variational approach to the simulation of hyperfine
effects in polyatomic molecules. Our results for water suggest that the strongest transitions with room-temperature intensities on the order of
$10^{-31}$ cm/molecule are about an order of magnitude larger than previously predicted values
and should be detectable in the mid-infrared $\nu_2$ and near-infrared
$2\nu_1+\nu_2$ and $\nu_1+\nu_2+\nu_3$ bands by current spectroscopy experiments.

§ INTRODUCTION

The water molecule is abundant in nature. It has two nuclear spin isomers, ortho, with a
total nuclear spin of hydrogen atoms $I=1$, and para, with a total nuclear spin of hydrogens
$I=0$. In isolated-molecule conditions the ortho and para nuclear spin isomers show
tremendously long-lasting stability to inter-conversion <cit.>, can be spatially separated <cit.>, and exhibit distinct physical and chemical
properties <cit.>. Thus the nuclear spin isomers of water
are frequently treated as distinct molecular species.

The concept of stable nuclear spin isomers is appealing to astrophysicists, as it allows to deduce
temperatures, below 50 K, in cometary comae, star- and planet-forming regions from the observations
of relative abundance of ortho and para species <cit.>. Some astronomical observations however reported anomalous
ortho-para ratios (OPR), corresponding to spin temperatures that are much lower than
gas kinetic temperatures in the same region <cit.>. These observations pose the intriguing question if
the OPR values could be altered as a result of internal conversion, which can possibly be
enhanced by natural factors, such as molecular collisions <cit.>, interaction with catalytic
surfaces <cit.>, external fields <cit.> and
radiation <cit.>. Low nuclear-spin temperatures have been attributed to
the photodesorption of water from colder icy grains <cit.>. However, this
theory was benchmarked and disputed in a number of recent laboratory
experiments <cit.>.
Arguably there could be another yet unknown mechanism of spin-non-destructive desorption of water
molecules from ice.

The OPR values can change as a result of the interaction between the nuclear spins and an induced
internal magnetic field of the rotating molecule, which is called the hyperfine spin-rotation
interaction. For the main water isotopologue H$_2^{16}$O, considered here, the $^{16}$O has zero
nuclear spin, and the hyperfine coupling between the spins of the protons is very weak, providing a
fundamental rationale for neglecting the conversion in practical applications. However,
it can be significantly enhanced by the accidental resonances between the ortho and
para states, produced by molecular collisions and interactions with strong external fields
and field gradients. The accurate modeling of these processes may unravel previously unknown
mechanisms contributing to the observed anomalous OPR of water in space. Precise knowledge of the
molecular hyperfine states and corresponding transitions is mandatory for the understanding of such
conversion mechanisms. This information can also be important for cold-molecule precision
spectroscopy relying on controlled hyperfine transitions and hyperfine-state changing
collisions <cit.>.

Here, we report a complete linelist of rovibrational hyperfine transitions in
H$_2$$^{16}$O at room-temperature that we computed using an accurate variational
approach <cit.> with an empirically refined potential energy surface
(PES) <cit.> and a high-level ab initio spin-rotation tensor surface.
The spin-spin coupling was modelled as the magnetic dipole-dipole interaction between the two
hydrogen nuclei. We show that the strongest forbidden transitions are on the order of
$10^{-31}$ cm/molecule, which is about ten times stronger than previously reported calculations for
the same lines <cit.>. We also present the details of our variational approach for
computing hyperfine effects, which is general and not restricted by the numbers and specific
magnitudes of the molecules' nuclear spins.

§ THEORETICAL DETAILS

§.§ Spin-rotation and spin-spin coupling

In this section we describe the implementation of the hyperfine spin-rotation and spin-spin coupling
terms within the general variational framework of the nuclear motion approach
TROVE <cit.>. Implementation details of the hyperfine nuclear quadrupole coupling can be
found in our previous works <cit.>.

The spin-rotation coupling is the interaction between the rotational angular momentum $\mathbf{J}$
of the molecule and the nuclear spins $\mathbf{I}_n$ of different
nuclei <cit.>

\begin{align}
\label{hsr_cart}
H_\text{sr} = \sum_{n}^{N_I} \mathbf{I}_n \cdot \mathbf{M}_n \cdot \mathbf{J},
\end{align}

where $\mathbf{M}_n$ is the second-rank spin-rotation tensor relative to the nucleus $n$ and the sum
runs over all nuclei $N_I$ with non-zero spin. The interaction between the nuclear spins
$\mathbf{I}_n$ of different nuclei is given by the spin-spin coupling as

\begin{align}
\label{hss_cart}
H_\text{ss} = \sum_{n>n'}^{N_I} \mathbf{I}_n \cdot \mathbf{D}_{n,n'} \cdot \mathbf{I}_{n'},
\end{align}

where $\mathbf{D}_{n,n'}$ is the second-rank spin-spin tensor, which is traceless and symmetric.
Using the spherical-tensor representation <cit.>, the spin-rotation and
spin-spin Hamiltonians can be expressed as

\begin{align}
\label{hsr_spher}
H_\text{sr} = &\frac{1}{2}\sum_{n}^{N_I} \sum_{\omega=0}^2 \sqrt{2\omega+1}
\left(-\frac{1}{\sqrt{3}}\right) \mathbf{I}_n^{(1)}\cdot \\
\nonumber
&\cdot\left( (-1)^\omega \left[ \mathbf{M}_n^{(\omega)} \otimes \mathbf{J}^{(1)} \right]^{(1)} +
\left[ \mathbf{J}^{(1)} \otimes \mathbf{M}_n^{(\omega)} \right]^{(1)} \right)
\end{align}

\begin{equation}
\label{hss_spher}
H_\text{ss} = \sum_{n>n'}^{N_I} \mathbf{D}_{n,n'}^{(2)} \cdot \left[ \mathbf{I}_n^{(1)}
\otimes \mathbf{I}_{n'}^{(1)} \right]^{(2)},
\end{equation}

where $\mathbf{M}_n^{(\omega)}$, $\mathbf{D}_{n,n'}^{(2)}$, $\mathbf{J}^{(1)}$, and
$\mathbf{I}_n^{(1)}$ denote the spherical-tensor representations of operators in (<ref>)
and (<ref>) and the square brackets are used to indicate the tensor product of two
spherical-tensor operators. Because the spin-rotation tensor is generally not symmetric, the second
term in the sum (<ref>) is added to ensure that the Hamiltonian is Hermitian.

The nuclear-spin operator $\mathbf{I}_n$ and the rotational-angular-momentum operator $\mathbf{J}$
are coupled using a nearly-equal coupling scheme, ,
$\mathbf{I}_{1,2}=\mathbf{I}_1+\mathbf{I}_2$, $\mathbf{I}_{1,3}=\mathbf{I}_{1,2}+\mathbf{I}_3$,
…, $\mathbf{I}\equiv\mathbf{I}_{1,N}=\mathbf{I}_{1,N-1}+\mathbf{I}_{N}$, and
$\mathbf{F}=\mathbf{J}+\mathbf{I}$. The nuclear-spin functions $\ket{I,m_I,\mathcal{I}}$ depend on
the quantum numbers $I$ and $m_I$ of the collective nuclear-spin operator $\mathbf{I}$ and its
projection onto the laboratory $Z$ axis, respectively. The set of auxiliary quantum numbers
$\mathcal{I}=\{I_1,I_{1,2},I_{1,3},\ldots,I_{1,N-1}\}$ for the intermediate spin angular momentum
operators provide a unique assignment of each nuclear-spin state. The total spin-rovibrational wave
functions $\ket{F,m_F,u}$ are built as symmetry-adapted linear combinations of the coupled products
of the rovibrational wave functions $\ket{J,m_J,l}$ and the nuclear-spin functions
$\ket{I,m_I,\mathcal{I}}$. Here, $J$ and $F$ are the quantum numbers of $\mathbf{J}$ and
$\mathbf{F}$ operators with $m_J$ and $m_F$ of their $Z$-axis projections. $l$ and $u$ denote the
rovibrational and hyperfine state indices, respectively, and embrace all quantum numbers, ,
rotational $k$ and vibrational quantum numbers $v_1,v_2,\ldots$, that are necessary to characterize
a nuclear spin-rovibrational state.

The symmetrization postulate requires the total wavefunction of the molecule to change sign
upon exchange of the protons, , to transform as one of the irreducible representations $B_1,B_2$
of its $\textbf{C}_\text{2v}$(M) symmetry group. Accordingly, the ortho spin state
$\ket{I=1}$ of $A_1$ symmetry can be coupled with the rovibrational states of $B_1$ and $B_2$
symmetries and the para state $\ket{I=0}$ of $B_2$ symmetry can be coupled with the
rovibrational states of $A_1$ and $A_2$ symmetries.

The matrix representations of the spin-rotation and spin-spin Hamiltonians in the basis of the
$\ket{F,m_F,u}$ functions are diagonal in $F$ and $m_F$, with the explicit expressions given by

\begin{align}\label{hsr_me}
\langle &F,m_F,u'|H_\text{sr}|F,m_F,u\rangle = \\ \nonumber
&=\frac{1}{2}(-1)^{I+F}\sqrt{(2J+1)(2J'+1)} \left\{ \begin{array}{ccc}I'&J'&F\\J&I&1\end{array}
\right\} \\ \nonumber
&\times\sum_{n}^{N_I}\sum_{\omega=0}^2 N_\omega
\left[
(-1)^\omega J \left\{ \begin{array}{ccc}\omega&1&1\\J&J'&J\end{array}
\right\}\left( \begin{array}{ccc}J&1&J\\-J&0&J\end{array} \right)^{-1}
\right. \\ \nonumber
& +\left.
J' \left\{ \begin{array}{ccc}1&\omega&1\\J&J'&J'\end{array}\right\}\left(
\begin{array}{ccc}J'&1&J'\\-J'&0&J'\end{array} \right)^{-1}
\right] \\ \nonumber
&\times\mathcal{M}_{\omega,n}^{(J'l',Jl)} \langle I'||\mathbf{I}_n^{(1)}||I\rangle
\end{align}

\begin{align}\label{hss_me}
\langle &F,m_F,u'|H_\text{ss}|F,m_F,u\rangle =  \\ \nonumber
&=(-1)^{I+J'+J+F}\sqrt{(2J+1)(2J'+1)} \left\{ \begin{array}{ccc}I'&J'&F\\J&I&2\end{array}
\right\} \\ \nonumber
\langle I'|| [ \mathbf{I}_n^{(1)} \otimes \mathbf{I}_{n'}^{(1)} ]^{(2)} ||I\rangle,
\end{align}

with the normalization constant $N_\omega=1$, $-\sqrt{3}$, and $\sqrt{5}$ for $\omega=0$, 1, and 2,
respectively. The expressions for the reduced matrix elements of the nuclear-spin operators
$\langle I'||\mathbf{I}_n^{(1)}||I\rangle$ and
$\langle I'|| [ \mathbf{I}_n^{(1)} \otimes \mathbf{I}_{n'}^{(1)} ]^{(2)} ||I\rangle$ depend on the
total number of coupled spins and can be computed using a general recursive procedure as described,
for example, in ref. Yachmenev:JCP147:141101. Here, for the two equivalent hydrogen
spins $I_1=I_2=1/2$, the reduced matrix elements are

\begin{align}
& \langle I'||\mathbf{I}_n^{(1)}||I\rangle = (-1)^{I\delta_{n,1}+I'\delta_{n,2}} I_1  \\
\nonumber
& \times \sqrt{(2I+1)(2I'+1)}
\left\{ \begin{array}{ccc}I_1&I'&I_1\\I&I_1&1\end{array} \right\}
\left( \begin{array}{ccc}I_1&1&I_1\\-I_1&0&I_1\end{array} \right)^{-1},
\end{align}

with the explicit values $\langle 0||\mathbf{I}_n^{(1)}||0\rangle=0$,
$\langle 1||\mathbf{I}_n^{(1)}||1\rangle=\sqrt{3/2}$,
$\langle 0||\mathbf{I}_n^{(1)}||1\rangle=\pm\sqrt{3}/2$ for $n=1$ and $2$, respectively, and
$\langle 1||\mathbf{I}_n^{(1)}||0\rangle=\mp\sqrt{3}/2$.

The expressions for the $\mathcal{M}_{\omega,n}^{(J'l',Jl)}$ and $\mathcal{D}_{n,n'}^{(J'l',Jl)}$
tensors in Eqs. (<ref>) and (<ref>) depend on the chosen rovibrational wave functions
$\ket{J,m_J,l}$, which are represented by the molecular rovibrational eigenfunctions calculated with
the variational approach TROVE. The functions $\ket{J,m_J,l}$ are linear combinations of products of
vibrational wave functions $\ket{\nu}=\ket{v_1,v_2,\ldots,v_{M}}$ ($M$ is the number of vibrational
modes) and symmetric-top rotational functions

\begin{align}
|J,m_J,l\rangle = \sum_{\nu,k} c_{\nu,k}^{(J,l)}
\ket{\nu}\ket{J,k,m_J}.
\end{align}

In this basis, the $\mathcal{M}_{\omega,n}^{(J'l',Jl)}$ and $\mathcal{D}_{n,n'}^{(J'l',Jl)}$ tensors are

\begin{align}\label{m_tens}
& \mathcal{M}_{\omega,n}^{(J'l',Jl)} = \sum_{\nu' k'}\sum_{\nu k}
\left[c_{\nu' k'}^{(J',l')}\right]^* \, c_{\nu k}^{(J,l)} \, (-1)^{k'} \\ \nonumber
& \times \sum_{\sigma=-\omega}^\omega\sum_{\alpha,\beta=x,y,z}
\left(\begin{array}{ccc} J & \omega & J' \\ k & \sigma & -k'\end{array}\right)
U_{\omega\sigma,\alpha\beta}^{(2)} \langle \nu'|\bar{M}_{\alpha\beta,n}|\nu\rangle
\end{align}

\begin{align}\label{d_tens}
& \mathcal{D}_{n,n'}^{(J'l',Jl)} = \sum_{\nu' k'}\sum_{\nu k}
\left[c_{\nu' k'}^{(J',l')}\right]^* \, c_{\nu k}^{(J,l)} \, (-1)^{k'} \\ \nonumber
& \times \sum_{\sigma=-2}^2\sum_{\alpha,\beta=x,y,z}
\left(\begin{array}{ccc} J & 2 & J' \\ k & \sigma & -k'\end{array}\right)
U_{2\sigma,\alpha\beta}^{(2)} \langle \nu'|\bar{D}_{\alpha\beta,nn'}|\nu\rangle
\end{align}

where $\bar{M}_{\alpha\beta,n}$ and $\bar{D}_{\alpha\beta,nn'}$ ($\alpha,\beta=x,y,z$) are
spin-rotation and spin-spin interaction tensors in the molecule-fixed frame and the $9\times9$
constant matrix $U_{\omega\sigma,\alpha\beta}^{(2)}$ ($\omega=0,\dots,2$,
$\sigma=-\omega,\ldots,\omega$) defines the transformation of a general second-rank
Cartesian tensor operator into its spherical-tensor representation, see,
, (5.41)–(5.44) in ref. Zare:AngularMomentum.

The total Hamiltonian $H$ is composed of a sum of the pure rovibrational Hamiltonian $H_\text{rv}$
and hyperfine terms $H_\text{sr}$ and $H_\text{ss}$. In the basis of TROVE wave functions, the
rovibrational Hamiltonian $H_\text{rv}$ is diagonal, its elements are given by the rovibrational

\begin{multline}
\label{tot_ham}
\langle F,m_F,u'|H|F,m_F,u\rangle \\
= E_u\delta_{u,u'} + \langle F,m_F,u'|H_\text{sr}|F,m_F,u\rangle \\
+ \langle F,m_F,u'|H_\text{ss}|F,m_F,u\rangle,
\end{multline}

where $\delta_{u,u'}=\delta_{J,J'}\delta_{l,l'}\delta_{I,I'}\delta_{\mathcal{I},\mathcal{I'}}$.

The above equations were implemented in the  module of the Python software package
Richmol <cit.>, which uses rovibrational molecular states
calculated in TROVE as a variational basis. Alternative approaches using Watson-type effective
Hamiltonians <cit.> are also implemented in the Richmol package.

The hyperfine energies and wave functions are computed in a three step procedure. First, we solve
the full rovibrational problem using TROVE and obtain the rovibrational energies and wave functions
for all states with energies below a selected threshold. In the next step, the rovibrational matrix
elements of the spin-rotation and spin-spin tensors are computed in the form given by
Eqs. (<ref>) and (<ref>). These matrix elements are later used to build the
spin-rotation and spin-spin interaction Hamiltonians using Eqs. (<ref>) and (<ref>).
The total Hamiltonian is composed of the sum of a purely rovibrational part, which is diagonal and
given by the rovibrational state energies, and non-diagonal spin-rotation and spin-spin parts. In
the final step, the hyperfine energies and wave functions are obtained by diagonalizing the total

The computation of the dipole transition intensities also proceeds in two steps. First, the
rovibrational matrix elements of the dipole moment surface are computed and cast into a tensor form
similar to (<ref>),

\begin{align}\label{k_tens}
& \mathcal{K}_{\omega}^{(J'l',Jl)} = \sum_{\nu' k'}\sum_{\nu k}
\left[c_{\nu' k'}^{(J',l')}\right]^* \, c_{\nu k}^{(J,l)} \, (-1)^{k'} \\ \nonumber
& \times \sum_{\sigma=-\omega}^\omega\sum_{\alpha,\beta=x,y,z}
\left(\begin{array}{ccc} J & \omega & J' \\ k & \sigma & -k'\end{array}\right)
U_{\omega\sigma,\alpha}^{(1)} \langle \nu'|\bar{\mu}_{\alpha}|\nu\rangle,
\end{align}

where $\bar{\mu}_{\alpha}$ ($\alpha=x,y,z$) is the permanent dipole moment in the molecule-fixed
frame and the $3\times3$ constant matrix $U_{\omega\sigma,\alpha}^{(1)}$ ($\omega=1$,
$\sigma=-\omega,\ldots,\omega$) defines the transformation of a general first-rank Cartesian tensor
operator into its spherical-tensor representation, see, ,
(5.4) in ref. Zare:AngularMomentum. In the second step, the dipole matrix elements are
transformed into the basis of hyperfine wave functions, ,

\begin{multline} {\mathcal{K}_{\omega}^{(F',u',F,u)} = \sum_{I',{\mathcal
I}',J',l'}\sum_{I,{\mathcal I},J,l} \left[ c_{I',{\mathcal
I}',J',l'}^{(F',u')} \right]^* c_{I,\mathcal{I},J,l}^{(F,u)} (-1)^{I}} \\
{\times\sqrt{(2J'+1)(2J+1)}\left\{\begin{array}{ccc}J'&F'&I \\ F&J&\omega\end{array}\right\}
\mathcal{K}_{\omega}^{(J',l',J,l)}\delta_{I',I}\delta_{\mathcal{I}',\mathcal{I}}},
\end{multline}

where $c_{I,\mathcal{I},J,l}^{(F,u)}$ are hyperfine wave function coefficients obtained by
diagonalization of the total Hamiltonian. Finally, the line strengths for transitions between
hyperfine states $\ket{f}=\ket{F',u'}$ and $\ket{i}=\ket{F,u}$ are computed
as <cit.>

\begin{align}
S(f\leftarrow i) = (2F'+1)(2F+1)\left|\mathcal{K}_{1}^{(F'u',Fu)} \right|^2,
\end{align}

where we sum over all degenerate $m_F$ and $m_F'$ components. The expression for the integrated
absorption coefficient of the dipole transition in units of cm/molecule reads

\begin{align}
I(f\leftarrow i) = \frac{8\pi^3 \nu_{if}
e^{-hcE_{i}/kT}\left(1-e^{-hc\nu_{if}/kT}\right)}{3hcZ(T)}S(f\leftarrow i),
\label{eq:intensity}
\end{align}

where $\nu_{if}=|E_i - E_f|$ is the transition wavenumber, $E_i$ and $E_f$ are energy term values of
the initial and final states in , $Z(T)$ is the temperature dependent partition function, $h$
(erg$\cdot$s) is the Planck constant, $c$ (cm/s) is the speed of light and $k$ (erg/K)
is the Boltzmann constant.

§.§ Electronic structure calculations

The molecule-fixed frame spin-rotation tensors $\bar{M}_{\alpha\beta,n}$ ($\alpha,\beta=x,y,z$,
$n=1,2$) were calculated ab initio on a grid of 2000 different molecular geometries with
electronic energies ranging up to 30 000 cm$^{-1}$ above the equilibrium energy. We used the
all-electron CCSD(T) (coupled-cluster singles, doubles, and perturbative triples) method with the
augmented core-valence correlation-consistent basis set aug-cc-pwCVTZ <cit.>
and aug-cc-pVTZ <cit.> for the oxygen and hydrogen atoms,
respectively. The basis sets were downloaded from the Basis Set Exchange
library <cit.>. The
calculations employed second-order analytical derivatives <cit.> together with
the rotational London orbitals <cit.>, as implemented in the
quantum chemistry package CFOUR <cit.>.

The electronic structure calculations used the principal axes of inertia coordinate frame. For
variational calculations another frame was employed, defined such that the $x$ axis is parallel to
the bisector of the valence bond angle with the molecule lying in the $xz$ plane at all instantaneous
molecular geometries. In this frame, the $z$ axis coincides with the molecular axis at the linear
geometry. The computed spin-rotation tensors were rotated from the principal axis of inertia to the
new frame. The permutation symmetry is such, that exchange of the two hydrogen atoms transforms
$\bar{M}_{\alpha\beta,1}$ into $\bar{M}_{\alpha\beta,2}$ followed by a sign change for non-diagonal
elements ($\alpha\neq\beta$).

The expression for the spin-rotation tensor, as computed in CFOUR, contains multiplication by the
inverse of the tensor of inertia, see (3) and (7) in ref. Gauss:MolPhys91:449. For
linear and closely linear geometries of the molecule, the inertial tensor becomes singular, which
creates a discontinuity in the dependence of $xz$ and $zz$ elements of spin-rotation tensor on the
bending angle. To circumvent this problem, we have multiplied the computed spin-rotation tensors on
the right side by the corresponding inertial tensors. The resulting data for the inertia-scaled
spin-rotation tensor was parameterized through least-squares fitting, using a power series
expansions to fourth order in terms of valence bond coordinates, with $\sigma_\text{rms}\leq0.3$ kHz
for all tensor components. Later, when computing the rovibrational matrix elements of the
spin-rotation tensor, we have multiplied the inertia-scaled tensor with the inverse moment of
inertia. The divergence of the spin-rotation tensor in the vicinity of linear geometries is exactly
canceled by the basis functions chosen to satisfy the kinetic cusp condition at the linear
geometry <cit.>.

The spin-spin tensor elements were computed as magnetic dipole-dipole interaction between two
hydrogen nuclei H$_1$ and H$_2$,

\begin{align}
D_{\alpha\beta,12} = \frac{\mu_0}{4\pi}\frac{\mu_1\mu_2}{I_1 I_2
r_{12}^3}\left(\mathbf{I}-3\mathbf{n} \otimes \mathbf{n}\right)_{\alpha\beta},
\end{align}

where $\mu_1=\mu_2=2.79284734$ are the magnetic dipole moments of H$_1$ and H$_2$ in units of the
nuclear magneton, $I_1=I_2=1/2$ are the corresponding hydrogen nuclear spins, $r_{12}$ is the
distance between the hydrogen nuclei, and $\mathbf{n}$ is the unit vector directed from one hydrogen
to another. The indirect spin-spin coupling constants mediated by the electronic motions were not
considered here, as they are typically two orders of magnitude smaller than the direct
constants <cit.>.

§.§ Nuclear motion calculations

We employed TROVE to calculate the rovibrational states using the exact kinetic-energy operator
formalism recently developed for triatomic molecules <cit.>. This
formalism is based on the use of associated Laguerre polynomials $L_{n}^{l}(x)$ as bending basis
functions, which ensures a correct behavior of the rovibrational wave functions at linear molecular
geometry <cit.>. The bisecting frame embedding was selected as a non-rigid
reference frame, with the $x$ axis oriented parallel to the bisector of the valence bond angle and
the molecule placed in the $xz$ plane. In this frame, the $z$ axis coincides with the linearity axis
at linear molecular geometry. Accurate empirically refined PES of H$_2^{16}$O was
employed <cit.>.

The primitive-stretching vibrational basis functions were generated by numerically solving the
corresponding one-dimensional Schrödinger equations on a grid of 2000 points using the
Numerov-Cooley approach <cit.>. The primitive basis
functions were then symmetry-adapted to the irreducible representations of the
$\textbf{C}_\text{2v}$(M) molecular symmetry group using an automated numerical procedure
<cit.>. The total vibrational basis set was formed as a direct product of the
symmetry-adapted stretching and bending basis functions, contracted to include states up to a polyad
48. It was used to solve the $J=0$ eigenvalue problem for the complete vibrational Hamiltonian of
H$_2$O. A product of the $J=0$ eigenfunctions and symmetry-adapted rigid rotor wavefunctions formed
the final rovibrational basis set. The rovibrational wavefunctions $\ket{J, m_J, l}$ for rotational
excitations up to $J=40$ and four irreducible representations $A_1$, $A_2$, $B_1$ and $B_2$ were
computed by diagonalizing the matrix representation of the total rovibrational Hamiltonian
$H_\text{rv}$ in the rovibrational basis set. More details about the variational approach and the
basis-symmetrization procedure for the case of triatomic molecules can be found in
ref. Yurchenko:JCP153:154106.

Strongest predicted transitions in H$_2$$^{16}$O at
$T=296$ K with the $10^{-31}$ cm/molecule intensity cut-off.

$\nu_1'$     $\nu_2'$     $\nu_3'$     $F'$     $J'$     $k_a'$     $k_c'$     $I'$     $E'$ ()
$\nu_1$     $\nu_2$     $\nu_3$     $F$     $J$     $k_a$     $k_c$     $I$     $E$
Freq. ()     Int. (cm/molec.)

0       1       0       3       4      2       3     $o$      1908.016319         0       0       0       4       4
4       0
$p$       488.134170      1419.882149     $2.26\times 10^{-31}$

0       1       0       3       3      3       1     $p$      1907.450231         0       0       0       3       4
3       2

$o$       382.516901      1524.933330     $1.36\times 10^{-31}$

0       1       0       3       3      3       1     $p$      1907.450231         0       0       0       3       4
1       4

$o$       224.838381      1682.611850     $1.12\times 10^{-31}$

0       1       0       3       4      2       3     $o$      1908.016319         0       0       0       3       3
2       2

$p$       206.301430      1701.714889     $1.02\times 10^{-31}$

0       1       0       3       3      3       1     $p$      1907.450231         0       0       0       2       3
1       2

$o$       173.365811      1734.084420     $2.05\times 10^{-31}$

0       1       0       3       4      2       3     $o$      1908.016319         0       0       0       2       2
2       0

$p$       136.163927      1771.852392     $3.28\times 10^{-31}$

2       1       0       3       4      1       4     $o$      8979.657423         0       0       0       4       4
1       3

$p$       275.497051      8704.160372     $3.36\times 10^{-31}$

2       1       0       3       4      1       4     $o$      8979.657423         0       0       0       3       3
1       3

$p$       142.278493      8837.378930     $1.01\times 10^{-31}$

2       1       0       3       4      1       4     $o$      8979.657423         0       0       0       2       2
1       1

$p$        95.175936      8884.481487     $6.41\times 10^{-31}$

1       1       1      15      14      3      11     $o$     11067.083574         0       0       0      14      14
0      14

$p$      2073.514207      8993.569367     $1.92\times 10^{-31}$

1       1       1      15      15      2      13     $p$     11067.089122         0       0       0      14      13
1      12

$o$      2042.309821      9024.779300     $2.04\times 10^{-31}$

§.§ Linelist simulations

The linelist of hyperfine rovibrational transitions for H$_2$$^{16}$O was computed with an energy
cutoff at $15\,000$ and includes transitions up to $F=39$ ($J=40$). To further improve the
accuracy of the linelist, after solving the pure rovibrational problem and before entering the
hyperfine calculations, the rovibrational energies $E_u$ in (<ref>) were replaced with the
high-resolution experimental IUPAC values from ref. Tennyson:JQSRT117:29, where
available. Such empirical adjustment of the rovibrational energies have been adopted and tested,
, for the production of molecular linelists as part of the ExoMol
project <cit.>. Recently, this approach was proven accurate for computing
the ultra-weak quadrupole transitions in water <cit.>
and carbon dioxide <cit.>, which enabled their
first laboratory (H$_2$O and CO$_2$) and astrophysical (CO$_2$) detection.

The final linelist has been calculated at room temperature
($T=296$ K) with the corresponding partition function
$Z=174.5813$ <cit.>, and a threshold of
$10^{-36}$ cm/molecule for the absorption intensity based on (<ref>). The linelist
stored in the ExoMol <cit.>

format is provided in the supplementary information.

§ RESULTS AND DISCUSSION

Overview of the H$_2^{16}$O dipole absorption spectrum at $T=296$ K. The
ortho-ortho and para-para transitions are marked with blue circles, whereas the
transitions are given by red circles.

Comparison of calculated hyperfine transitions (red stems) with experimental data
(dashed lines) from (a) ref. Bluyssen:PLA9:482 and (b-h)
ref. Cazzoli:CPL473:21. Different panels show hyperfine transitions for different
rotational bands $J'_{k_a',k_c'}\leftarrow J_{k_a,k_c}$. The measured (calculated)
zero-crossing frequencies, in MHz, are 22235.0447 (22235.0322), 321225.6363 (321225.6311),
380197.3303 (380197.3361), 439150.7746 (439150.7857), 443018.3358 (443018.4016), 448001.0538
(448001.0359), 556935.9776 (556935.9849), 620700.9334 (620700.8889) for panels (a)–(h),

An overview of the calculated H$_2^{16}$O dipole absorption stick spectrum at $T=296$ K is shown in
<ref>. The forbidden ortho-para transitions are plotted as red circles. Despite
being, at least, 10 orders of magnitude weaker than the corresponding allowed transitions, for some
of the strongest ortho-para transitions the predicted absorption intensities are close to the
sensitivity threshold of modern cavity ring-down spectroscopic techniques <cit.>. All predicted ortho-para transitions
with line intensity larger than $10^{-31}$ cm/molecule are listed in <ref>.
These transitions all occur in the
fundamental $\nu_2$ bending and the overtone $2\nu_1 + \nu_2$ and $\nu_1+\nu_2+\nu_3$ bands. The
off-diagonal elements of molecular-frame spin-rotation tensor $\bar{M}_{\alpha\beta, n}$, which lead
to interaction, are highly dependent on the bending vibrational coordinate, indicating
significance of the $\nu_2$ band in transitions. The size of the off-diagonal
spin-rotation matrix elements increases for bending angles close to 180, , the linear
geometry of the molecule. This leads to an increase in the interaction for rovibrational
energies close to the linearity barrier at $\ordsim8254~\invcm$ above the zero-point energy. The
spin-rotation coupling in these vibrationally excited states is responsible for the transitions. For example, the final transition state $F=3$, $J_{k_a,k_c}=4_{2,3}$ (ortho) with
energy $E=1908.016319$ is mixed with the state $F=3$, $J_{k_a,k_c}=3_{3,1}$ (para) with
energy $E=1907.450231$ . The size of the rovibrational matrix element of spin-rotation tensor,
$\mathcal{M}_{\omega,n}^{(J'l',Jl)}$ in (<ref>) for this transitions is $\pm0.95$ kHz and
$\pm6.3$ kHz ($\pm$ for $n=1,2$) for $\omega=1$ and 2, respectively. Note that following
(<ref>) only the spin-rotation tensor with $\omega=1$ contributes to the coupling.
Allowed transitions into these states from the ground state are quite strong, $2.07\times10^{-20}$
and $3.52\times10^{-20}$ cm/molecule, respectively. Accordingly, intensity borrowing as a result of
the spin-rotation interaction of excited states leads to non-zero intensities of the two
corresponding forbidden transitions on the order of $10^{-31}$ molecule/cm. Similarly for other of
the strongest forbidden transitions listed in <ref>, the enhancement occurs due to
intensity borrowing effect from strongly allowed transitions with coincident near resonance between
the excited states, accompanied by a relatively large value of the spin-rotation matrix element

Absorption cross sections computed at $T=296$ K for selected rotational bands, using
Gaussian lineshapes with HWHMs of 0.01 (solid lines), 0.005 (dashed lines),
and 0.001 (dotted lines). The cross sections for allowed ortho-ortho and
para-para transitions are plotted with blue colour lines and cross sections for
forbidden ortho-para transitions are plotted with red colour lines.

Though ortho-para transitions are yet to be observed in , there are several spectroscopic
studies of the allowed hyperfine transitions in the pure rotational spectrum of
$\text{H}_2^{16}\text{O}$ <cit.>. We used these data to validate the accuracy of our predictions.
In <ref> the calculated transitions (stems) are compared with the available experimental
data (dashed lines), demonstrating an excellent agreement, within 1–4 kHz, for the hyperfine
splittings. For example, the root-mean square (rms) deviation of the predicted hyperfine splittings
from experiment is 2.1 kHz in <ref>a]fig2, while for the absolute line positions it is
12.3 kHz. The latter can be explained by the discrepancies in predictions of the pure rotational
transitions. The errors in predictions of the hyperfine splittings can be attributed to the level of
electronic structure theory, in particular the basis set, employed in the calculations of
spin-rotation tensor surface. The basis set convergence of the equilibrium spin-rotation constants
of was investigated elsewhere <cit.>. According to the results, the employed
aug-cc-pwCVTZ basis set produces an average error of 1.3 kHz with a maximum of 1.8 kHz for one of
the off-diagonal elements, when compared with the results obtained with the aug-cc-pwCV6Z basis set.
There are several predicted splittings in <ref>d–h]fig2 that are less than 12 kHz and were not
resolved in the experiment <cit.>. Indeed, by visual inspection of the Lamb-dip
spectrum plotted in Fig. 1 of ref. Cazzoli:CPL473:21, which was provided as an example
of the experimental resolution achieved in that work, the transition profiles' full width at a half
maximum is about 13 kHz.

The sensitivity and resolution required to observe the transitions in a prospective
experiment can be estimated from the simulated absorption spectrum, shown <ref> for
selected wavenumber ranges with strong transitions. Since the Doppler linewidth would be
around 0.01 at room temperature and even much higher-resolution spectroscopy was
demonstrated <cit.>, we used simple Gaussian line profiles with half-width at
half-maximum (HWHM) fixed at 0.01, 0.005, and 0.001 and computed absorption cross sections at
$T=296$ K using  <cit.> to predict the experimental spectra.
The transitions In <ref>a,c]fig3 (red) show considerable overlap with the allowed
transitions (blue) for purely rotational transitions and in the fundamental $\nu_2$ excitation band
and could only be detected with an experimental HWHM below 0.005 at an experimental
sensitivity of $10^{-30}$ and $10^{-29}$ cm$^2$/molecule, respectively. In <ref>b,d]fig3,
showing parts of the $\nu_2$ and $\nu_2+\nu_3$ bands, the predicted transitions are
better separated from the allowed transitions and should already be detectable at lower resolution,
, at HWHM of 0.01 , but demand a greater sensitivity of $10^{-30}$ and
$10^{-31}$ cm$^2$/molecule, respectively. Such high-sensitivity measurements of intensities on the
scale of $10^{-30}$ cm$^2$/molecule are currently within reach, for example, using continuous wave
laser cavity ring down spectroscopy <cit.>.

§ CONCLUSIONS

We developed and performed comprehensive variational calculations of the room temperature linelist
of with hyperfine resolution, including forbidden transitions. The calculations were
based on accurate rovibrational energy levels and wavefunctions produced using the variational
approach TROVE. The nuclear hyperfine effects were modeled as spin-rotation and direct spin-spin
interactions, with the spin-rotation coupling surface calculated at a high level of the
electronic-structure theory. We found excellent agreement between the calculated transition
frequencies and available hyperfine-resolved spectroscopic data of allowed transitions.

The predicted transitions are useful for guiding future experimental spectroscopic
studies in search of these forbidden transitions in the laboratory as well as in astrophysical
environments. Our accurate predictions of hyperfine effects complement the spectroscopic data for

The variational approach we developed for computing these hyperfine effects is general. It includes
nuclear quadrupole <cit.>, spin-rotation, and
spin-spin interactions, and can be applied to other molecular systems without restrictions on the
number and values of nuclear spins.

§ SUPPLEMENTARY MATERIAL

The computed hyperfine-linelist data for are available at

§ AUTHOR DECLARATIONS

§.§ Conflict of interests

The authors have no conflicts to disclose.

§ DATA AVAILABILITY

The computer codes used in this work are available from git repositories at
<https://github.com/Trovemaster/TROVE> and <https://github.com/CFEL-CMI/richmol>.

§ ACKNOWLEDGMENTS

We acknowledge support by Deutsches Elektronen-Synchrotron DESY, a member of the Helmholtz
Association (HGF). This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the
priority program “Quantum Dynamics in Tailored Intense Fields” (QUTIF, SPP1840, YA 610/1) and the
cluster of excellence “Advanced Imaging of Matter” (AIM, EXC 2056, ID 390715994) and through the
Maxwell computational resources operated at Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany.
S.Y. acknowledges support from the UK Science and Technology Research Council (STFC,
No. ST/R000476/1) and the European Research Council under the European Union’s Horizon 2020 research
and innovation programme through an Advanced Grant (883830). The authors acknowledge the use of the
Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of
Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC
component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and
ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.
G.Y. gratefully acknowledges the financial support by the China Scholarship Council (CSC).


H$\beta$ red peak | $51.6^{+7.0}_{-7.9}$ | $54.5^{+4.5}_{-11.0}$ | $0.80\pm 0.19$ | $53.1^{+7.1}_{-9.4}$ | $54.5^{+6.0}_{-10.0}$ | $0.86\pm 0.74$
H$\beta$ red wing | $54.4^{+6.0}_{-10.3}$ | $56.0^{+4.5}_{-11.5}$ | $0.81\pm 0.35$ | $56.3^{+5.9}_{-9.6}$ | $57.0^{+5.5}_{-8.0}$ | $0.91\pm 0.71$
H$\beta$ blue vs. red wing | $2.9^{+13.4}_{-18.5}$ | $3.0^{+18.0}_{-22.0}$ | $0.47\pm 0.10$ |  |  |
H$\beta$ blue vs. red peak | $1.0^{+8.8}_{-9.8}$ | $0.0^{+10.0}_{-10.0}$ | $0.50\pm 0.11$ |  |  |
Table 10: Estimates of the black hole virial product M${}_{bh}^{vir}$, based
on the Balmer and helium line profile properties ($\sigma_{line}$ and
separation of the blue and red peak in the double-peaked profiles) and the
time delays $\tau_{cent}$ (80 % threshold) of the broad emission-line flux
variations. Virial Black Hole Product $M_{bh}^{vir}$ in $10^{8}M_{\odot}$
---
using $\sigma_{line}$ of the line profiles
feature | Fλ(5100Å) |  g-band
| M${}_{bh}^{vir}$(mean) | M${}_{bh}^{vir}$(rms) | M${}_{bh}^{vir}$(mean) | M${}_{bh}^{vir}$(rms)
(1) | (2) | (3) | (4) | (5)
H$\alpha\,\lambda$ 6563 | $2.33^{+0.17}_{-0.16}$ | $2.57^{+0.25}_{-0.24}$ | $2.33^{+0.10}_{-0.27}$ | $2.57^{+0.20}_{-0.34}$
H$\beta\,\lambda$ 4861 | $2.62^{+0.20}_{-0.18}$ | $2.69^{+0.29}_{-0.27}$ | $2.50^{+0.17}_{-0.19}$ | $2.57^{+0.25}_{-0.27}$
H$\gamma\,\lambda$ 4340 | $0.99^{+0.37}_{-0.33}$ | $1.67^{+0.62}_{-0.55}$ | $1.81^{+0.14}_{-0.19}$ | $3.05^{+0.23}_{-0.33}$
He ii $\lambda 4686$ | $2.08^{+0.43}_{-0.41}$ | $1.26^{+0.36}_{-0.35}$ | $1.89^{+0.56}_{-0.33}$ | $1.14^{+0.40}_{-0.30}$
average | $2.01^{+0.67}_{-0.48}$ | $2.05^{+0.38}_{-0.36}$ | $2.13^{+0.24}_{-0.21}$ | $2.33^{+0.35}_{-0.37}$
using the separation of the blue and red peak
feature | Fλ(5100Å) |  g-band
| M${}_{bh}^{vir}$(mean) | M${}_{bh}^{vir}$(rms) | M${}_{bh}^{vir}$(mean) | M${}_{bh}^{vir}$(rms)
(1) | (2) | (3) | (4) | (5)
H$\alpha\,\lambda$ 6563 | $1.74^{+0.14}_{-0.13}$ | $1.91^{+0.21}_{-0.21}$ | $1.74^{+0.09}_{-0.21}$ | $1.91^{+0.18}_{-0.27}$
H$\beta\,\lambda$ 4861 | $1.41^{+0.12}_{-0.10}$ | $1.87^{+0.23}_{-0.22}$ | $1.34^{+0.10}_{-0.11}$ | $1.78^{+0.21}_{-0.22}$
H$\gamma\,\lambda$ 4340 | $1.10^{+0.41}_{-0.37}$ | $0.88^{+0.42}_{-0.39}$ | $2.02^{+0.19}_{-0.23}$ | $1.61^{+0.50}_{-0.51}$
average | $1.42^{+0.35}_{-0.32}$ | $1.55^{+0.48}_{-0.46}$ | $1.70^{+0.26}_{-0.33}$ | $1.77^{+0.29}_{-0.31}$


# Designing a Pseudo-Random Bit Generator with a Novel 5D-Hyperchaotic System

Ngoc T. Nguyen, Toan Q. Bui, Ghyslain Gagnon, _Member, IEEE_ , Pascal Giard,
_Senior Member, IEEE_ ,
and Georges Kaddoum, _Member, IEEE_ The authors are with the LaCIME,
Department of Electrical Engineering, École de technologie supérieure,
Montreal, QC, Canada, (e-mails<EMAIL_ADDRESS>quang-
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>georges.kaddoum@etsmtl.ca).

###### Abstract

Dynamic and non-linear systems are emerging as potential candidates for random
bit generation. In this context, chaotic systems, which are both dynamic and
stochastic, are particularly suitable. This paper introduces a new continuous
chaotic system along with its corresponding implementation, which targets
field-programmable gate arrays (FPGAs). This chaotic system has five
dimensions, which exhibit complex chaotic dynamics, thus enabling the
utilization of chaotic signals in cryptography. A mathematical analysis is
presented to demonstrate the dynamic characteristics of the proposed
hyperchaotic system. A novel digital implementation of the proposed system is
presented. Moreover, a data scrambling circuit is implemented to eliminate the
bias effect and increase the randomness of the bitstream generated from the
chaotic signals. We show that the proposed random bit generator has high
randomness. The generated bits successfully pass well-known statistical
randomness test-suites, i.e., NIST SP800-22, Diehard and TestU01. The ready-
to-use random bit generator is deployed on a Xilinx Zynq-7000 SoC ZC702
Evaluation Kit. Experimental results show that the proposed random bit
generator can achieve a maximum throughput of 6.78 Gbps, which is over 3.6
times greater than state-of-the-art designs, while requiring under 4% of the
resources available on the targeted FPGA.

###### Index Terms:

Chaos, hyperchaos, pseudo-random number generator, NIST, TestU01,
cryptography, security, FPGA, System Generator.

##  I Introduction

Random number generators are critical components that are responsible for
generating public keys, private keys, and other kinds of random numbers that
are utilized in cryptographic applications and security [1, 2, 3].
Differential chaotic systems, which present high dynamic characteristics and
multi-dimensional signals, are superior in terms of generating random bits due
to their ability to achieve a high level of randomness [4]. Therefore, in this
work, we focus on the design of such systems, targeted at field-programmable
gate arrays (FPGAs).

A hyperchaotic system exhibits rich dynamics since the system states it hosts
are expanded exponentially in several directions simultaneously. This property
makes the hyperchaotic system an interesting candidate for the generation of
random keys used in miscellaneous applications in engineering, such as secure
communications, cryptosystems, and encryptions [5, 6]. Therefore, in the
present work, we develop a 5D hyperchaotic system with three positive Lyapunov
exponents to provide better dynamic characteristics than the state of the art.
Moreover, high-dimensional chaotic systems provide multiple outputs which
improve the throughput of the overall random bit generator.

The implementation methodology adopted has a great impact on the
digitalization of differential chaotic systems. The research in [7] presented
an implementation of numerical techniques, including the Euler, mid-point, and
fourth-order Runge-Kutta (RK) methods on FPGAs. Among these, the Euler method
has the shortest data path, but the least accuracy. The fourth-order RK method
shows the highest accuracy, but also has the longest time frame. In state-of-
the-art RK algorithm implementations, mapping functions are used four times
and implemented separately [8, 9]. We propose the fourth-folding Runge-Kutta
(FFRK) method to allow reusing the mapping functions (the latter are
implemented only once and reused four times), and as a result, the device
resources required to implement the mapping functions are reduced by 75%. The
iterations are controlled by adding MUXs and control signals, which use
insignificant amounts of device resources. Moreover, taking advantage of the
multiple dimensional signals in hyperchaotic systems, we implement an
effective post-processing, in which five chaotic outputs are used to generate
ready-to-use random bitstreams. In summary, the contributions of this work
include i) a new 5D hyperchaotic system to generate stochastic signals, ii) a
novel implementation of the fourth-order RK algorithm, and iii) a simple and
effective data post-processing circuit.

The remainder of this paper is organized as follows. Related works are
reviewed in Section II. The mathematical model of the hyperchaotic system is
clarified in Section III. The implementation of the proposed hyperchaotic
system in a FPGA hardware device using the new FFRK method is detailed in
Section IV. The experimental resuts are presented and discussed in Section V.
Then, the randomness evaluation of the random bit generator is presented in
Section VI. Finally, Section VII summarizes the work in this paper.

##  II Related Works

In this section, we review the differential chaotic systems and their
implementations in the digital world to emphasize our research contribution.
There are several commonly used differential chaotic systems in secure
communications, including the Lorenz system, Chua’s system, Liu’s system and
Lu’s system [10, 11]. The Lorenz system has been used in different
applications in many research works [10]. The research in [10] presents the
implementation of a Lorenz system in FPGA hardware devices and co-simulation
with Matlab. Liu’s system is implemented in [11] using the Grunward-Letniknov
algorithm. However, the above chaotic systems utilize multiple multiplexers,
which are resource-demanding in hardware implementations. Moreover, the
authors in [12] provide both analog and digital implementations of the Bruke-
Shaw chaotic system in a Virtex-6 FPGA chip, which consumes a lot of device
resources. The research in [13] presents an FPGA realization of Chua’s system
with multiple scrolls by using different saturated functions. The work in [8]
implements a 3D chaotic system with no equilibrium points. A high-speed FPGA
implementation of a 3D continuous chaotic system, which achieves a maximum
operating frequency of 293 MHz, is presented in [9].

The authors in [14] present a new hybrid PRNG using Sprott 94 G chaotic system
which improves security of the advanced encryption standard (AES) in
cryptographic systems. Although instructive, these chaotic systems are
3-dimensional systems which has only one positive Lyapunov exponent. This
means that thay have limited dynamic characteristics as the signal expands
exponentially in only one direction. Speed optimization is applied to a 3D
chaotic system and then expanded to a 4D chaotic system in [15], achieving a
maximum throughput of 1882 Mbps for the random bit generator.

Nowadays, differential chaotic systems achieve higher complexity by conveying
integer-order systems into the fractional-order domain [16]. However,
fractional-order systems are complicated to implement in hardware design due
to their memory dependency. The hardware implementation of fractional-order
differentiators and integrators requires careful consideration [11].

Although chaotic systems are unpredictable, and have random-like state
trajectories, they can be studied and recovered by using computational tools.
A cyber-attack has a high possibility of success if the target system uses a
well-known and self-excited oscillator for its random bit generator. After the
transient process, a trajectory, starting from the point of an unstable
manifold in a small neighborhood of unstable equilibrium, can be revealed.
Therefore, the system’s parameters can be computed to recover the target
system. From this point of perspective, the development of modern computers
enables the numerical simulation of complex nonlinear dynamical systems, and
therefore the structure of their trajectories can be deduced. However, this
approach shows very limited success when it comes to hidden attractor chaotic
systems [17, 18]. The hyperchaotic systems with hidden attractors in [19, 20]
were proposed to overcome such attacks. As compared to previous systems, a
mathematical analysis of our hyperchaotic system with hidden attractors shows
rich dynamic characteristics that are suitable for use in security and
cryptographic applications.

##  III 5D Differential Chaotic System Model

The following section presents the proposed hidden attractor hyperchaotic
system, which is expressed by five differential equations as given in (1). The
theoretical analysis is divided into three parts. The first part presents the
proposed chaotic system and a theoretical study of chaotic characteristics,
while the second part addresses the stability of the equilibrium. The
simulation of chaotic transitions is investigated in the third part.

###  III-A Mathematical Analysis

The proposed 5D dynamic system, which is developed from the 4D hyperjerk
chaotic system presented in [21], is expressed as $S^{\prime}=F(S)$,
$S=(x,y,z,u,v)\in\mathbb{R}^{5}$, with

$F(S)=\begin{cases}x^{\prime}=y,\\\ y^{\prime}=z,\\\ z^{\prime}=u,\\\
u^{\prime}=-z-0.5\times u+(x-1)\times y\\\ v^{\prime}=-u-0.5\times
v+(x-1)\times z.\ \end{cases}$ (1)

The ordinary differential equations (ODEs) are solved and simulated in MATLAB,
based on the fourth-order RK integration algorithm with a step size of
$10^{-2}$. Here, a small step size is chosen to provide a higher resolution
and better accuracy. The equilibrium points of the proposed system are located
on the line $E(x,0,0,0,0)$. Therefore, the proposed system is a dynamical
system with hidden attractors. According to [20], it is difficult to expose a
chaotic attractor, and then reveal the chaotic system architecture by choosing
an arbitrary initial condition. In other words, the chaotic attractors are
invisible to attackers and their basin localization does not dissolve the
chaotic oscillator. In what follows, we present a mathematical analysis of the
novel chaotic system in terms of Lyapunov exponents and bifurcation diagrams.

  1. Lyapunov exponents: The Lyapunov exponents of the system are defined as: $L_{1}=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{\|\partial x(t)\|}{\|\partial x(0)\|}$, $L_{2}=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{\|\partial y(t)\|}{\|\partial y(0)\|}$, $L_{3}=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{\|\partial z(t)\|}{\|\partial z(0)\|}$, $L_{4}=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{\|\partial u(t)\|}{\|\partial u(0)\|}$, and $L_{5}=\lim_{t\rightarrow\infty}\frac{1}{t}\log\frac{\|\partial v(t)\|}{\|\partial v(0)\|}$ The Lyapunov exponents of the novel 5D chaotic system are: $L_{1}=0.093790,L_{2}=0.001101,L_{3}=0.000107,L_{4}=-0.500100$, and $L_{5}=-0.594898$. The initial conditions of the system are chosen as $x_{0}=0.0002,y_{0}=0.0005,z_{0}=0.00005,v=0.001,\text{and }v_{0}=0$. The divergence of (1) is evaluated based on the following conditions:

$\begin{split}\frac{\partial\dot{x}}{\partial
x}+\frac{\partial\dot{y}}{\partial y}+\frac{\partial\dot{z}}{\partial
z}+\frac{\partial\dot{u}}{\partial u}+\frac{\partial\dot{v}}{\partial
v}=-1<0,\\\ L=L_{1}+L_{2}+L_{3}+L_{4}+L_{5}<0.\end{split}$ (2)

  2. Lyapunov dimension: The Lyapunov dimension ($D_{L}$), which is closely related to the correlation dimension, is commonly used to evaluate the chaotic complexity [22]. A higher Lyapunov dimension value is associated with a chaotic system with a higher level of complexity.

$D_{L}=j+\frac{1}{|L_{j+1}|}\sum_{i=1}^{j}L_{i},$ (3)

where $j$ is the largest index of the positive Lyapunov exponent. In the
proposed system, $j=3$, and the Lyapunov dimension is therefore:

$D_{L}=3+\frac{L_{1}+L_{2}+L_{3}}{|L_{4}|}=3.1899.$ (4)

Figure 1: Poincaré section of the phase-space in the plain ($x=1$).

Fig. 1 shows the Poincaré section of the phase-space in the plane $x=1$.
Distinct set of points in the Poincaré section indicates the chaotic region of
dynamic system.

###  III-B Stability of Equilibria

The stability evaluation of the equilibria is an important step to find the
chaotic region of a dynamic system. To evaluate the stability of equilibria,
the Jacobian matrix at the equilibria $E=[c,0,0,0]$ of the proposed
hyperchaotic system is calculated. The eigenvalues of the Jacobian matrix
satisfy the condition $J_{E}-\lambda I=0$.

$J_{E}=\begin{bmatrix}0&1&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\
0&(c-1)&-1&-0.5&0\\\ 0&0&(c-1)&-1&-0.5\\\ \end{bmatrix}.$ (5)

Assuming a small perturbation from the fixed points $x(t)=x(0)+\Delta x$. If
$\Delta x\approx e^{\lambda}t$, the characteristic polynomial equation is
derived as:

$\lambda(\lambda+0.5)(\lambda^{2}+1.5\lambda+(c-0.5))=0.$ (6)

This equation indicates that the Jacobian matrix at equilibrium points has one
zero value and three non-zero values. According to the Routh-Hurwitz
criterion, apart from the zero eigenvalues, the real parts of the roots of (6)
are negative if and only if $(c-0.5)>0$. To make the equilibrium set $E$
unstable, thereby enabling the possibility of chaos occurrence, the initial
condition $x(0)=c$ must satisfy $c<0.5$. In other words, depending on the
initial value, the proposed system has stable or unstable saddle-focus points.
Thus, the dynamical behaviors of the line equilibrium chaotic system are
heavily dependent on the initial state of the variable $x(t)$.

###  III-C Transitions to Chaotic Region

Figure 2: Bifurcation diagram according the initial condition. Figure 3:
Lyapunov spectrum according to initial condition x(0)=c. Figure 4: Transition
to chaos regions depending to initial condition $x(0)=c$ : a) stable region
for $c=0.6$, b) limited dynamic region for $c=0.4$, c) chaotic region with
limited periods for $c=0.2$, and d) rich dynamic and chaotic region for
$c=0.05$.

The bifurcation diagram and the Lyapunov spectrum of the state variable $x(t)$
of the proposed system are shown in Fig. 2 and Fig. 3, respectively. Both
figures illustrate that the value of the initial condition $c$ has an impact
on the characteristics of the system. Moreover, a stable region is observed
for $0.5<c<0.92$ in Fig. 2, where the data space converges. This corresponds
to the negative Lyapunov exponents in Fig. 3. A dynamic and unbounded region
is observed for $c>0.92$, a region that corresponds to positive high-value
Lyapunov exponents in Fig. 3. Finally, in the $0<c<0.5$ range, the system is
dynamic and bounded in a limited data space, as illustrated in the bifurcation
diagram. This corresponds to the presence of positive Lyapunov exponents
$L_{1}$ and/or $L_{2}$, as shown by the Lyapunov spectrum. However, the system
only exhibits rich and chaotic dynamics for $c<0.05$, where a high positive
value of $L_{1}$ is presented, as shown in Fig. 3. The gradual transition from
periodic orbits to chaotic regions based on various values of the initial
condition $c$ is observed in Fig. 4. The chaotic region with rich dynamics
illustrated under Fig. 4-d) is obtained with $c=0.05$ which is compatible to
the stability analysis of equilibria in the previous section.

Figure 5: The proposed pseudo-random number generator scheme.

##  IV Pseudo-Random Number Generator Scheme

In this section, we present a pseudo-random number generator scheme based on
the chaotic output signals. We use fixed-point 32-bits to represent data, with
1-bit for the sign and 27 bits for the fractions. Fig. 5 shows the scheme of
the proposed random bit generator which includes the implementation of the 5D
chaotic system block using FFRK with 32-bit fixed-point data. The truncation
block that follows keeps the 12 least-significant bits of the stream. Then,
the data is up-sampled and converted from parallel to serial before entering
the post-processing block. In the post-processing block, a data scrambler is
proposed using a simple LFSR of one of the chaotic outputs, while the random
bits are obtained by XORing the output of the data scrambler with the other
chaotic signals.

Due to the long data path of the RK algorithm and low self-oscillator
frequency of the chaotic system, signal truncation is applied in order to have
better randomness qualification in the following block. Therefore, 32-bit
chaotic outputs are truncated, where only 12 least significant bits (LSB) are
used in the following post-processing. 12-LSBs are up-sampled and serialized
before the post-processing step. In what follows, we present a novel FPGA-
based implementation of the fourth-order RK algorithm in a new manner to solve
the differential chaotic equations. Moreover, we present the post-processing
process with a data scrambler to increase the randomness of the generated
bitstreams and eliminate bias effects.

Figure 6: Hyperchaos function $F(S_{n})$ implementation in Xilinx System
Generator. Figure 7: Block diagram of the proposed fourth-folding Runge-Kutta
(FFRK) algorithm. Clock signals are omitted for clarity. Figure 8:
Trajectories of the chaotic outputs.

###  IV-A Implementation of the 5D Hyperchaotic System

To solve differential equations, there are three well-known numerical methods,
with different complexity and accuracy. These are, the Euler method, the
Midpoint method, and the RK method. Compared to the other two candidates, the
RK method has the longest calculation path and requires larger areas [7].
However, it produces the most accurate results. The proposed FFRK
implementation optimizes the resources of the 5D hyperchaotic system. Compared
to previous RK implementations, the proposed FFRK provides the most accurate
solution for differential equations and consumes fewer device resources. The
system state at the $n$-th time step is denoted by the vector
$S_{n}=[x_{n},y_{n},z_{n},u_{n},v_{n}]$. The proposed hyperchaotic system
discrete-time implementation is calculated as follows:

$F(S_{n})=\begin{cases}F_{x}=y_{n}\\\ F_{y}=z_{n}\\\ F_{z}=u_{n}\\\
F_{u}=-z_{n}-0.5\times u_{n}+(x_{n}-1)\times y_{n}\\\ F_{v}=-u_{n}-0.5\times
v_{n}+(x_{n}-1)\times z_{n}.\end{cases}\\\ $ (7)

The fourth-order RK algorithm is obtained by defining the primitive parameters
$k_{1}$, $k_{2}$, $k_{3}$, and $k_{4}$ which are calculated as:

$\begin{split}&k_{1}=F(S_{n})\\\ &k_{2}=F(S_{n}+\frac{1}{2}h\times k_{1})\\\
&k_{3}=F(S_{n}+\frac{1}{2}h\times k_{2})\\\ &k_{4}=F(S_{n}+h\times
k_{3}),\end{split}$ (8)

where $h$ is the discrete step size. The next system state vector $S_{n+1}$ is
then evaluated using the previous system state $S_{n}$ and these above
mentioned primitive parameters. Fig. 6 presents the block diagram
implementation of $F(S_{n})$ in (7), which is reused to implement the
parameters $k_{i}$ with different input variables. The high-level tool (XSG)
automatically maps those functions to FPGA resources. The $z^{-i}$ notation
indicates the insertion of $i$ registers. We can see that the longest path
amounts to $5$ registers, which results in $5$ clock cycles of output latency
for the function $F(S_{n})$. The block diagram presented in Fig. 7 shows the
implementation of the FFRK algorithm. The discrete step size $h=0.01$ provides
enough resolution for chaotic signals. First, the parameter $k_{1}$ is
calculated while the register enable control signal $E_{1}$ is high.

$S_{n+1}=\begin{bmatrix}x_{n+1}\\\ y_{n+1}\\\ z_{n+1}\\\ u_{n+1}\\\
v_{n+1}\end{bmatrix}=S_{n}+\frac{1}{6}h\times(k_{1}+2k_{2}+2k_{3}+k_{4}).$ (9)

The output $k_{1}$ from the block $F(S_{n})$ will be temporary stored in
register $R_{1}$. Then, the parameters $k_{2}$, $k_{3}$, and $k_{4}$ are
produced at the high level of the control signals $E_{2}$, $E_{3}$ and
$E_{4}$, respectively. Therefore, four primitive parameters
$k_{1},k_{2},k_{3}$, and $k_{4}$ are stored in four registers. Finally, the
next state of the system is produced when all these parameters are ready at
the high period of the enable signal $E_{0}$ of register $R_{0}$. The enable
signals are generated by the control unit which uses delay blocks. The phase
spaces of the chaotic outputs of the proposed implementation shown in Fig. 8
are compatible with the simulation results and mathematical analysis.

###  IV-B Data Post-Processing Implementation

Before we address the post-processing process, we truncate 32-bit chaotic
signals as indicated in Fig. 5. In previous works, the truncation was
mentioned and applied [23, 24]. In [23], random bits are generated from the
$33^{rd}$ bit to the $40^{th}$ bit in each 52-bit data. The authors in [24]
used 16-bits out of a 32-bit signal for the post-processing. It is obvious
that the full data-width in chaotic signal is not used to generate random bits
due to the low self-oscillators. Nevertheless, to obtain sufficient randomness
quality, no clear reason is presented for choosing the truncation position.
The least significant bits (LSBs) are preferred due to their higher
fluctuation. Therefore, it is understood that the bits going from the LSB to
the truncation position will be used in the post-processing. Our research
shows that the signal entropy is strongly related to the truncation position.
To evaluate the truncated chaotic signal’s entropy, Shannon’s entropy of a
bitstream is calculated. For a given random variable $\tilde{X_{n}}$, with a
possible outcome $x_{i}$, $i\subset[0,N-1]$, each with probability $p_{i}$,
the signal entropy is:

$H(\tilde{X_{n}})=-\sum_{i=0}^{N-1}p_{i}\log_{2}p_{i},$ (10)

where $N$ is the number of symbols. The average entropy per bit is calculated
as:

$\bar{E}=\frac{H(\tilde{X_{n}})}{N_{b}}$ (11)

in which, $\tilde{X_{n}}$ denotes the truncated chaotic signal and
$N_{b}=\log_{2}N$ is the truncation position. Fig. 9 provides the average
entropy per bit with different lengths of truncated signals. The average
entropy per bit decreases significantly when the number of remaining bits is
larger than 12. Therefore, in this work, a 12-bit truncated chaotic vector
$\tilde{S_{n}}=[\tilde{X_{n}},\tilde{Y_{n}},\tilde{Z_{n}},\tilde{U_{n}},\tilde{V_{n}}]$
is used in the data post-processing.

Figure 9: Average entropy per bit with different length of truncated chaotic
signals. Figure 10: Post-processing design. Clock signals omitted for clarity.

The proposed data post-processing hardware implementation based on shift
registers and XOR operations is presented in Fig. 10. The data scrambler is
composed of four shift registers (m-shift registers) and XOR operations [25,
26]. The working principle is to evaluate the incoming bits from the truncated
chaotic signals and to then reuse them. One of the chaotic signal outputs of
the proposed 5D hyper-chaotic system ($\tilde{V_{n}}$) is used to generate the
seed for the data scrambler. Then, the other serialized truncated chaotic
signals are XOR-ed with the output bit from the data scrambler. Based on the
experimental result, $m=6$ shift registers, which are implemented by delay
units (flip-flops (FFs)) in the Xilinx System Generator, are utilized to
obtain high success rates in the statistical tests, . Five random bitstreams
$(B_{1}\rightarrow B_{5})$ are collected for the evaluation phase, which will
be presented below.

##  V Experimental Results

In this section, we present the deployment of the random bit generator on the
Xilinx Zynq-7000 SoC ZC702 Evaluation Kit. The implementation is performed
using the Xilinx System Generator tool. The random number generator is
generated as an IP-core and evaluated on Xilinx Vivado tools with the default
configuration.

###  V-A Comparison and Discussion

We start this section with a comparison with the state-of-the-art chaotic
systems, i.e., a Lorenz system as well as both Elwaki’s systems of [27]. To
illustrate the advantage of our FFRK implementation, we carried out our own
implementation of the systems, in which we produced versions that either
integrated the traditional RK implementation or our proposed FFRK method. The
implementation is produced on the evaluation kit mentioned above. Four
iterations of the mapping functions are used to generate $k_{i}$ in the
traditional fourth-order RK method.

Fig. 11 shows the amount of FPGA resources required for the various chaotic
systems. The Lorenz systems are denoted LRZ-RK and LRZ-FFRK, while the
Elwaki’s systems are denoted EWK1-RK, EWK2-RK, EWK1-FFRK, and EWK2-FFRK; the
RK and FFRK suffixes respectively indicate whether the system uses either the
traditional RK implementation or our proposed FFRK implementation. From Fig.
11, it can be seen that compared to traditional RK methods [28, 7], the
proposed FFRK requires 4 times less DSPs, under half the number of look-up
tables (LUTs), and approximately 15% less FFs.

Figure 11: Amount of FPGA resources in terms of LUTs, flip-flops, and DSPs
required by various continuous chaotic systems, using the traditional RK
implementation and the proposed FFRK.

Fig. 11 also includes results for the proposed 5D chaotic system, both with
the traditional RK method (denoted 5D-RK) and the proposed FFRK method
(denoted 5D-FFRK). Using the proposed FFRK leads to a significant reduction in
the amount of LUTs and DSPs required, approximately the same as that seen for
the other systems. There is a small increase, ($\sim 25\%$), in the amount of
FFs required.

TABLE I: FPGA resource comparison for various implementations of continuous chaos-based PRNGs. pseudo-random number generator (PRNG) | [9] | [11] | [29] | [30] | [15] | [31] | Proposed
---|---|---|---|---|---|---|---
Chaotic system | 3D | 3D (Liu) | 3D | 3D | 3D(Lorenz+Lu) | 3D (opt. t/p) | 4D (opt. t/p) | 3D | 5D
Implementation | fourth-order RK | GL | GL | ADM | Euler | Euler | Euler | Euler | FFRK
Slices | 12430 | 5688 | 2659 | – | 119 | – | – | – | 920
LUTs | 43732 | – | – | 1220 | 494 | 1169 | 1882 | 1355 | 2017
FFs | 42092 | 4962 | 3360 | 192 | 118 | 416 | 480 | 1318 | 3458
DSPs | – | 99 | 59 | 8 | 8 | – | – | – | 8
Max. Freq. (MHz) | 293 | 38 | 46 | 87 | 78 | 107 | 112 | 464 | 113
Throughput (Mbps) | 58.7 | 1554 | 1849 | – | 1872 | 1178 | 1869 | 464 | 6780

In addition to the above, Table I compares our proposed hyperchaos-based PRNG
implementation against the implementations of chaos-based PRNGs from the
literature. In this context, we implement our system using the FFRK method
using the default configurations for synthesis and implementation strategies.
It is noted that changing the synthesis tool optimization policy to favor
speed will lie a small improvement. However, the greatest gain will come from
duplicating the design along with some straightforward circuitry to merge
their results, effectively doubling the throughput at the cost of doubling the
amount of required resources without any effect on latency.

In terms of resource requirements, it can be observed that the 3D chaotic
system using the fourth-order RK algorithm in [9], which utilized up to 32% of
the available resources from the Virtex-6 XC6VLX240T-1-FF1156 FPGA chip,
requires more resources than the other systems. Compared to our work, the
Grünwald-Letnikov (GL)-based algorithms of [11] require a similar number of
FFs and from $3\times$ to $6.4\times$ the number of slices. While the
implementation of the 3D-fractional-order chaotic system using the Adomian
decomposition method (ADM) from [29] requires a smaller number of slices, the
random number generator is only performed in software design. The 3D-chaos
based PRNG, which is itself based on Lorenz and Lu chaotic systems from [30]
using the Euler method, uses the least number of resources among all systems,
but offers an achievable throughput of $1872$ Mbps, which is much lower than
that of our proposed system. The 4D-chaos based system from [15], which is
optimized for throughput, is the one that comes the closest to our system in
terms of the number of LUTs, where our implementation achieves $3.6\times$
greater throughput. The system in [31] uses the highest number of resources
among all the 3D-chaos based PRNG implementations using the Euler method. We
recall that the hardware implementation of the Euler method requires fewer
device resources, but achieves less accuracy as compared to the fourth-order
Runge-Kutta algorithm [7].

Our design seeks to reach a high throughput, and a high level of randomness
while requiring a modest amount of resources. The Euler method implementation
provides less accuracy, and therefore the average entropy per bit is decreased
and so is the randomness level for the output random bits. We also note that
for the targeted FPGA, the proposed system requires under 4% of the available
LUTs, FFs, or DSP blocks. The design works at a maximum frequency of 113 MHz,
which enables a maximum throughput of 6.78 Gbps for the generator. The
estimated power consumption of the random number generator is<EMAIL_ADDRESS>Therefore, the estimated energy efficiency is 10.8 pJ/b. The FFRK
implementation for the chaotic system has a delay of 65 clock cycles while
there is only one latency at the random outputs.

Figure 12: 12b-truncated chaotic signals in a) X-Y distribution points b)
Histogram diagram. Before post-processing and c) Histogram of output
bitstreams - After post-processing.

##  VI Randomness Performance Evaluation

In this section, we present multiple tests performed on the harvested binary
bitstreams before and after applying post processing. These bitstreams are
collected by hardware co-simulation of Vivado System Generator and Matlab. The
PRNG core runs on the FPGA board connected to the computer and the output data
are collected through JTAG. The experimental setup is illustrated in Fig. 13.
The random binary outputs are observed using an oscilloscope.

Figure 13: The experimental validation of the proposed chaos-based PRNG at the
50MHz clock frequency.

###  VI-A Histogram

First, $100K$ length of 12b-truncated chaotic signals are collected and used
in distribution and histogram analysis. As shown in Fig. 12-a), since there
are no obvious patterns in their distribution, 12b-truncated chaotic signals
have good randomness characteristics. Fig. 12-b) shows the histogram of the
signal $\tilde{X_{n}}$ before post-processing where several peaks that can be
observed. The histogram of the output bitstream $B_{1}$ after post-processing
is illustrated in Fig. 12-c), which indicates a random distribution, where no
apparent pattern can be discerned. Similar results are observed for $B_{2}$ to
$B_{5}$.

###  VI-B Standard statistical random testsuites

In order to be used in cryptography applications, the random binary bits
should be tested using statistical tests which require long bitstreams. The
NIST SP800-22 test-suite, which is developed and introduced by The National
Institute of Standards and Technology (NIST), is commonly used [32]. The NIST
test results are presented in Table II. Moreover, Table III provides the
results for the Diehard tests which consist of 12 sub-tests. TestU01 is
another commonly used test-suite providing a variety of statistical tests for
random bit generators. We applied two battery tests including the Rabbit and
the Alphabit tests, to the binary sequences from the proposed random bit
generator. The Rabbit test includes 39 subtests and the Alphabit test includes
17 subtests for a bitstream of $2^{25}$ bit length. The battery test results
are presented in Table IV. Most of the tests are passed, except for the
multinomial test of bitstream $B_{3}$, with a P-value = 0.00097, which is very
close to the threshold (0.001). It should not be understood that if a RNG
fails in some statistical tests, it cannot be used in practical problems; even
the best commercial RNGs experience failures in very complicated tests [33].

TABLE II: NIST SP800-22 tests for the proposed PRNG, $N=10^{3}$ sequences of
the length of 1Mb.

NIST SP-800.22 | P-value(*) | Proportion
---|---|---
Monobit test | 0.323668 | 998/1000
Frequency within block test | 0.763677 | 997/1000
Runs test | 0.164425 | 992/1000
Longest run 1’s test | 0.630872 | 1000/1000
Rank test | 0.575608 | 990/1000
DFT test | 0.655608 | 996/1000
Cumulative sum | 0.924844 | 997/1000
Overlapping template | 0.130366 | 990/1000
Non-overlapping template | 0.999999 | 997/1000
Linear complexity test | 0.663130 | 998/1000
Maurers universal test | 0.999142 | 998/1000
Approximate entropy | 0.956970 | 1000/1000
Serial | 0.740523 | 1000/1000
Random excursions | 0.304039 | 997/1000
Random excursion variant | 0.049499 | 988/1000

(*): the average value.

TABLE III: Diehard tests for the proposed PRNG, sequence of 5M bits. Dieharder | P-value | Result
---|---|---
Birthdays | 0.0528 | PASSED
OPERM5 | 0.0000 | FAILED
Rank 6x8 | 0.0042 | WEAK
Bitstream | 0.9566 | PASSED
OPSO | 0.0109 | PASSED
DNA | 0.1829 | PASSED
Count 1s string | 0.9989 | PASSED
Parking lot | 0.2452 | PASSED
2d sphere | 0.1491 | PASSED
3d sphere | 0.1581 | PASSED
Sums | 0.0222 | PASSED
Runs | 0.9288 | PASSED
Craps | 0.4094 | PASSED
TABLE IV: Battery test results for $2^{25}$ output binary bits. TestU01 | Rabbit | Alphabit
---|---|---
$B_{1}$ | 39/39 | 17/17
$B_{2}$ | 39/39 | 17/17
$B_{3}$ | 38/39 | 17/17
$B_{4}$ | 39/39 | 17/17
$B_{5}$ | 39/39 | 17/17

##  VII Conclusion

We have presented a novel 5D hyperchaotic system along with its implementation
on FPGA. In this work, we proposed a novel method, which reduces device
resource usage, to implement the fourth-order Runge-Kutta numerical technique
to solve differential equations in chaotic systems. The PRNG, which is
designed to take advantage of the high dimension and the high accuracy of the
hyperchaotic system’s hardware implementation, reaches a maximum throughput of
6.78 Gbps. The design, which is deployed in the Xilinx FPGA chip, occupies
under 4% of the available FPGA hardware resources. The random generated bits
have been verified in modern statistical tests confirming that our proposed
system is suitable for industrial applications. Future work includes
integrating the proposed PRNG in real-time communication applications such as
the secure streaming of images and videos.

## References

  * [1] B. Li, X. Liao, and Y. Jiang, “A novel image encryption scheme based on improved random number generator and its implementation,” _Nonlinear Dyn._ , vol. 95, no. 3, pp. 1781–1805, Feb. 2019.
  * [2] Y. Liu, Z. Qin, X. Liao, and J. Wu, “Cryptanalysis and enhancement of an image encryption scheme based on a 1-D coupled Sine map,” _Nonlinear Dyn._ , vol. 100, no. 3, pp. 2917–2931, May. 2020.
  * [3] L. Moysis, A. Tutueva, C. Volos, and D. Butusov, “A Chaos Based Pseudo-Random Bit Generator Using Multiple Digits Comparison,” _Chaos Theory and Applications_ , vol. 2, no. 2, p. 11, Jul. 2020.
  * [4] N. Nguyen, L. Pham-Nguyen, M. B. Nguyen, and G. Kaddoum, “A low power circuit design for chaos-key based data encryption,” _IEEE Access_ , vol. 8, pp. 104 432–104 444, 2020.
  * [5] Y. Liu and X. Tong, “Hyperchaotic system-based pseudorandom number generator,” _IET Information Security_ , vol. 10, no. 6, pp. 433–441, Nov. 2016.
  * [6] M. Irfan, A. Ali, M. A. Khan, M. Ehatisham-ul Haq, S. N. Mehmood Shah, A. Saboor, and W. Ahmad, “Pseudorandom Number Generator (PRNG) Design Using Hyper-Chaotic Modified Robust Logistic Map (HC-MRLM),” _Electronics_ , vol. 9, no. 1, p. 104, Jan. 2020\.
  * [7] M. A. Zidan, A. G. Radwan, and K. N. Salama, “The effect of numerical techniques on differential equation based chaotic generators,” in _ICM 2011 Proceeding_ , pp. 1–4, Hammamet, Tunisia, Dec. 2011.
  * [8] A. Akgul, H. Calgan, I. Koyuncu, I. Pehlivan, and A. Istanbullu, “Chaos-based engineering applications with a 3D chaotic system without equilibrium points,” _Nonlinear Dyn._ , vol. 84, no. 2, pp. 481–495, Apr. 2016.
  * [9] I. Koyuncu and A. Turan Özcerit, “The design and realization of a new high speed FPGA-based chaotic true random number generator,” _Computers & Electrical Engineering_, vol. 58, pp. 203–214, Feb. 2017.
  * [10] L. Zhang, “System generator model-based FPGA design optimization and hardware co-simulation for Lorenz chaotic generator,” in _2017 2nd Asia-Pacific Conf. on Intelligent Robot Systems (ACIRS)_ , pp. 170–174. Wuhan, China: IEEE, Jun. 2017.
  * [11] M. F. Tolba, A. M. AbdelAty, N. S. Soliman, L. A. Said, A. H. Madian, A. T. Azar, and A. G. Radwan, “FPGA implementation of two fractional order chaotic systems,” _AEU - Int J Electron Commun_ , vol. 78, pp. 162–172, Aug. 2017.
  * [12] I. Koyuncu, A. T. Ozcerit, and I. Pehlivan, “An analog circuit design and FPGA-based implementation of the Burke-Shaw chaotic system,” _Optoelectron. Adv. Mater. Rapid Commun._ , p. 5, 2013.
  * [13] E. Tlelo-Cuautle, J. Rangel-Magdaleno, A. Pano-Azucena, P. Obeso-Rodelo, and J. Nunez-Perez, “Fpga realization of multi-scroll chaotic oscillators,” _Communications in Nonlinear Science and Numerical Simulation_ , vol. 27, no. 1-3, pp. 66–80, Oct. 2015.
  * [14] E. Avaroglu, I. Koyuncu, A. B. Ozer, and M. Turk, “Hybrid pseudo-random number generator for cryptographic systems,” _Nonlinear Dyn._ , vol. 82, no. 1-2, pp. 239–248, Oct. 2015.
  * [15] T. Bonny, R. Al Debsi, S. Majzoub, and A. S. Elwakil, “Hardware Optimized FPGA Implementations of High-Speed True Random Bit Generators Based on Switching-Type Chaotic Oscillators,” _Circuits, Systems, and Signal Processing_ , vol. 38, no. 3, pp. 1342–1359, Mar. 2019.
  * [16] A. Akgul, C. Arslan, and B. Aricioglu, “Design of an Interface for Random Number Generators based on Integer and Fractional Order Chaotic Systems,” _Chaos Theory and Applications_ , vol. 1, no. 1, pp. 1–18, Nov. 2019, number: 1.
  * [17] K. Barati, S. Jafari, J. C. Sprott, and V.-T. Pham, “Simple chaotic flows with a curve of equilibria,” _Int J Bifurcat Chaos_ , vol. 26, no. 12, p. 1630034, Nov. 2016.
  * [18] T. Fozin Fonzin, K. Srinivasan, J. Kengne, and F. Pelap, “Coexisting bifurcations in a memristive hyperchaotic oscillator,” _AEU - Int J Electron Commun_ , vol. 90, pp. 110–122, Jun. 2018.
  * [19] H. Bao, N. Wang, B. Bao, M. Chen, P. Jin, and G. Wang, “Initial condition-dependent dynamics and transient period in memristor-based hypogenetic jerk system with four line equilibria,” _Communications in Nonlinear Science and Numerical Simulation_ , vol. 57, pp. 264–275, Apr. 2018\.
  * [20] Ü. Çavuşoğlu, S. Panahi, A. Akgül, S. Jafari, and S. Kacar, “A new chaotic system with hidden attractor and its engineering applications: analog circuit realization and image encryption,” _Analog Integr Circ Sig Process_ , vol. 98, no. 1, pp. 85–99, Jan. 2019.
  * [21] V.-T. Pham, S. Vaidyanathan, C. K. Volos, S. Jafari, and X. Wang, “A chaotic hyperjerk system based on memristive device,” vol. 636, p. 21.
  * [22] K. E. Chlouverakis and J. Sprott, “A comparison of correlation and Lyapunov dimensions,” _Physica D: Nonlinear Phenomena_ , vol. 200, no. 1-2, pp. 156–164, Jan. 2005.
  * [23] Z. Hua, B. Zhou, and Y. Zhou, “Sine Chaotification Model for Enhancing Chaos and Its Hardware Implementation,” vol. 66, no. 2, pp. 1273–1284, Feb. 2019.
  * [24] Y. Mao, L. Cao, and W. Liu, “Design and FPGA Implementation of a Pseudo-Random Bit Sequence Generator Using Spatiotemporal Chaos,” in _2006 Int. Conf. on Communications, Circuits and Systems_ , pp. 2114–2118. Guilin, Guangzi, China: IEEE, Jun. 2006.
  * [25] F. Pareschi, R. Rovatti, and G. Setti, “Simple and effective post-processing stage for random stream generated by a chaos-based RNG,” in _The 2006 Int. Symp. on Nonlinear Theory and its Applications (NOLTA2006)_ , p. 5, 2006\.
  * [26] V. Rozic, B. Yang, W. Dehaene, and I. Verbauwhede, “Iterating von neumann’s post-processing under hardware constraints,” in _2016 IEEE Int. Symp. on Hardware Oriented Security and Trust (HOST)_ , pp. 37–42. IEEE, 2016-05.
  * [27] A. Elwakil and M. Kennedy, “Construction of classes of circuit-independent chaotic oscillators using passive-only nonlinear devices,” vol. 48, no. 3, pp. 289–307, Mar. 2001.
  * [28] M. S. Azzaz, C. Tanougast, S. Sadoudi, and A. Dandache, “Real-time FPGA implementation of Lorenz’s chaotic generator for ciphering telecommunications,” in _2009 Joint IEEE North-East Workshop on Circuits and Systems and TAISA Conf._ , pp. 1–4. Toulouse, France: IEEE, Jun. 2009.
  * [29] K. Rajagopal, A. Akgul, S. Jafari, and B. Aricioglu, “A chaotic memcapacitor oscillator with two unstable equilibriums and its fractional form with engineering applications,” _Nonlinear Dyn_ , vol. 91, no. 2, pp. 957–974, Jan. 2018.
  * [30] A. A. Rezk, A. H. Madian, A. G. Radwan, and A. M. Soliman, “Reconfigurable chaotic pseudo random number generator based on FPGA,” _AEU - Int J Electron Commun._ , vol. 98, pp. 174–180, Jan. 2019.
  * [31] I. Koyuncu, M. Tuna, I. Pehlivan, C. B. Fidan, and M. Alcin, “Design, FPGA implementation and statistical analysis of chaos-ring based dual entropy core true random number generator,” _Analog Integr Circ Sig Process_ , vol. 102, no. 2, pp. 445–456, Feb. 2020.
  * [32] A. Rukhin _et al._ , _A statistical test suite for random and pseudorandom number generators for cryptographic applications_ , National Institute of Standards and Technology (NIST) Special Publication 800-22, Rev. 1a, Apr. 2010.
  * [33] P. L’Ecuyer and R. Simard, “Testu01: A c library for empirical testing of random number generators,” _ACM Trans. Math. Softw._ , vol. 33, no. 4, Aug. 2007.