vol. XXX:XXX 2023 XXX 10.46298/dmtcs.XXX 2023-1-27XXX

# Algorithms for ranking and unranking
the combinatorial set of RNA secondary structures

Yuriy Shablya1  Dmitry Kruchinin2  Tomsk State University of Control Systems
and Radioelectronics, Tomsk, Russia

###### Abstract

In this paper, we study the combinatorial set of RNA secondary structures of
length $n$ with $m$ base-pairs. For a compact representation, we encode an RNA
secondary structure by the corresponding Motzkin word. For this combinatorial
set, we construct an AND/OR tree structure, find a bijection between the
combinatorial set and the set of variants of the AND/OR tree, and develop
algorithms for ranking and unranking the variants of the AND/OR tree. The
developed ranking and unranking algorithms have polynomial time complexity
$O(m^{2}(n-m))$ for $m<n-2m$ and $O(m(n-m)^{2})$ for $m>n-2m$. In contrast to
the existing algorithms, the new algorithms do not require preprocessing steps
and have better time complexity.

###### keywords:

combinatorial generation, AND/OR tree, RNA secondary structure, Motzkin word,
ranking, unranking

##  1 Introduction

Combinatorial generation is a branch of science that lies at the intersection
of computer science and combinatorics and develops methods and algorithms for
ranking and unranking combinatorial sets. Fundamentals of combinatorial
generation can be found in Kreher and Stinson (1999), Ruskey (2003) and Knuth
(2011). Many information objects can be represented as a combinatorial set,
for example, in the form of permutations, partitions, trees, paths, and
others. In this case, it is possible to use combinatorial generation methods
for processing such sets, their encoding, generation, and storage. This is
especially true when processing large combinatorial sets (when studying
different chemical and biological structures).

In this article, we propose to consider the structure of RNA (ribonucleic
acid). RNA is one of the main and most important elements of the cellular
structure of any organism. An RNA molecule is a chain of nucleotides
containing four nitrogenous bases: A (adenine), C (cytosine), G (guanine), and
U (uracil). The sequence of nucleotides in the RNA chain forms the RNA primary
structure. In mathematical terms, an RNA primary structure is a sequence of
elements taken from the set $\\{A,C,G,U\\}$, i.e. $a=(a_{1},a_{2},\ldots)$
where $a_{i}\in\\{A,C,G,U\\}$.

In addition, hydrogen bonds can be formed between two bases in the RNA chain
(for example, Watson-Crick base-pairs for A–U and C–G). The formation of
hydrogen bonds in an RNA molecule leads to different loops in it, which
complicates the structure of the RNA molecule. The obtained shape of the RNA
chain with hydrogen bonds forms an RNA secondary structure. In mathematical
terms, an RNA secondary structure can be represented in the form of a graph,
where the nodes of the graph are nucleotides and the edges of the graph show
the connections between bases. For more details on the features of RNA
structures, see Saenger (1984), Reidys (2011) and Waterman (2018).

RNA molecules play an important role in cellular nucleic acid processing and
gene expression. According to Picardi (2015), a deeper understanding of these
processes can be obtained through the detailed study of RNA structures. For
example, one such research direction is the prediction of RNA secondary
structures: Zuker and Sankoff (1984); Oluoch et al. (2019); Sato et al.
(2021); Sato and Kato (2022). Moreover, the improvement of such research works
through the development of computational algorithms is one of the main
directions of modern bioinformatics.

Let us consider the combinatorial set of RNA primary structures, which will be
denoted by $A_{n}$. The set of RNA primary structures is described by one
parameter $n$ (the number of nucleotides in an RNA chain). This combinatorial
set can be represented as the set of $n$-permutations of $m$ elements with
repetitions where $m=4$. That is, $n$ elements are selected from the set of
$m=4$ elements $\\{A,C,G,U\\}$ and these elements can be repeated. Thus,
algorithms for ranking and unranking the set of such permutations can also be
used for the combinatorial set of RNA primary structures.

Next, we consider the combinatorial set of RNA secondary structures, which
will be denoted by $A_{n,m}$. This combinatorial set is described by two
parameters: $n$ is the number of nucleotides in an RNA chain, $m$ is the
number of base-pairs. There are also researches related to the development of
combinatorial generation algorithms for the set of RNA secondary structures:
Nebel et al. (2011); Seyedi-Tabari et al. (2010); Mohammadi and Nowzari-Dalini
(2018).

Nebel et al. (2011) presented algorithms for generating a set of $m$ random
RNA secondary structures of length $n$ (without specifying the number of base-
pairs). These combinatorial generation algorithms were developed using the
method of Flajolet et al. (1994). As well as in Viennot and Vauchaussade de
Chaumont (1985), an RNA secondary structure was encoded by a correct bracket
sequence with the possibility of adding spaces (or a Motzkin word). The
developed algorithms have time complexity $O(mn\log{n})$ and require a
preprocessing step with time complexity $O(n^{2})$. The main disadvantage of
this approach is the lack of control over the number of base-pairs in the
generated RNA secondary structures, as well as the absence of ranking and
unranking algorithms.

Seyedi-Tabari et al. (2010) presented algorithms for listing, ranking, and
unranking the combinatorial set of RNA secondary structures of length $n$ with
$m$ base-pairs. Then Mohammadi and Nowzari-Dalini (2018) presented a parallel
version of the listing algorithm. These combinatorial generation algorithms
were developed using an approach based on the calculation of sequence
prefixes, that is similar to the method of Ryabko (1998). In this case, an RNA
secondary structure was encoded by a tree with $n$ internal nodes and $m$
external nodes, which was also encoded by an $E$-sequence. The ranking and
unranking algorithms have polynomial time complexity $O(nm-m^{2})$ and
$O(n^{2}+m^{2})$, but both of these algorithms require a preprocessing step
with time complexity $O((n-m)|A_{n,m}|)$, where $|A_{n,m}|$ is the cardinality
of the set $A_{n,m}$. Since $|A_{n,m}|$ can be calculated using binomial
coefficients (see Equation (2)), then the preprocessing step has exponential
time complexity. Thus, the preprocessing step is the main disadvantage of this
approach, since it has exponential time complexity and requires a large memory
space to store the results of preprocessing calculations.

In this paper, the main purpose of our research is to develop new effective
algorithms (with polynomial time complexity and without preprocessing steps)
for ranking and unranking the combinatorial set of RNA secondary structures of
length $n$ with $m$ base-pairs.

##  2 Main results

###  2.1 Research method

There are different approaches for developing combinatorial generation
algorithms. In this paper, we propose to apply the method for developing
combinatorial generation algorithms that is based on using AND/OR trees and
described in Shablya et al. (2020). To apply this method for developing
ranking and unranking algorithms, it is necessary to represent a given
combinatorial set in the form of an AND/OR tree structure for which the total
number of its variants is equal to the value of the cardinality function of
the combinatorial set.

An AND/OR tree is a tree structure that contains nodes of two types: AND nodes
and OR nodes (edges to the sons of an AND node are connected by an additional
curve in Figure 2). A variant of an AND/OR tree is a tree structure obtained
by removing all edges except one for each OR node. If we know the cardinality
function $f$ of a combinatorial set $A$ and $f$ belongs to the algebra
$\\{\mathbb{N},+,\times,R\\}$ (usage of only natural numbers, addition and
product operations, and the primitive recursion operator), then we can
construct the corresponding AND/OR tree structure. Hence, we need to know the
required form of the cardinality function for the considered combinatorial set
of RNA secondary structures.

###  2.2 Combinatorial set

For a compact representation, we encode an RNA secondary structure by the
corresponding Motzkin word. Figure 1 presents an example of an RNA secondary
structure of length $n=27$ with $m=9$ base-pairs.

Figure 1: An example of representing an RNA secondary structure of length
$n=27$ with $m=9$ base-pairs.

Figure 1(a) shows the ”real” shape of the RNA secondary structure as a chain
of nucleotides (nodes connected by solid lines) with base-pairs (nodes
connected by dashed lines). In Figure 1(b), the RNA secondary structure is
presented as the straight line of the RNA chain. In Figure 1(c), this RNA
secondary structure is represented as the corresponding Motzkin word in the
form of a correct bracket sequence with the possibility of adding spaces
(denoted by ’*’). Hence, we can encode an RNA secondary structure by a
sequence $a~{}=~{}(a_{1},\ldots,a_{n})$, where $a_{i}~{}\in~{}\\{$ ’(’, ’)’,
’*’ $\\}$. In this case, ’(’ and ’)’ denote a base-pair, ’*’ denotes a
nucleotide without hydrogen bonds.

The cardinality function of the considered combinatorial set of RNA secondary
structures $A_{n,m}$ is defined by Schmitt and Waterman (1994):

$f(n,m)=|A_{n,m}|=S_{n}^{m}=S_{n-1}^{m}+\sum\limits_{i=0}^{m-1}{\sum\limits_{j=2i}^{n-2(m-i)-1}{S_{n-2-j}^{m-1-i}S_{j}^{i}}},$
(1)

where $S_{n}^{0}=1$ for $n\geq 0$ and $S_{n}^{m}=0$ for $m\geq\frac{n}{2}$.

The values of $S_{n}^{m}$ form the integer sequence $A089732$ in OEIS
Foundation Inc. (2022) and can be calculated using the following explicit
formula by Schmitt and Waterman (1994):

$S_{n}^{m}=\frac{1}{n-m}\binom{n-m}{m}\binom{n-m}{m+1}$ (2)

for $n>m\geq 0$ and $S_{0}^{0}=1$. Table 1 presents the first values of
$S_{n}^{m}$.

Table 1: Several first values of $S_{n}^{m}$. | $m$
---|---
$n$ |  0 |  1 |  2 |  3 |  4
0 |  1 |  |  |  |
1 |  1 |  |  |  |
2 |  1 |  |  |  |
3 |  1 |  1 |  |  |
4 |  1 |  3 |  |  |
5 |  1 |  6 |  1 |  |
6 |  1 |  10 |  6 |  |
7 |  1 |  15 |  20 |  1 |
8 |  1 |  21 |  50 |  10 |
9 |  1 |  28 |  105 |  50 |  1
10 |  1 |  36 |  196 |  175 |  15

Generating functions play an important role in the study of integer sequences.
The following bivariate generating function of the numbers $S_{n}^{m}$ is
known from OEIS Foundation Inc. (2022):

$S(x,y)=\sum\limits_{n\geq 0}{\sum\limits_{m\geq
0}{S_{n}^{m}x^{n}y^{m}}}=\frac{1-x+x^{2}y-\sqrt{(1-x-x^{2}y)^{2}-4x^{3}y}}{2x^{2}y}.$

There is a relationship between the numbers $S_{n}^{m}$ and the Narayana
numbers $N_{n}^{m}$. It is also worth noting that the Narayana numbers
$N_{n}^{m}$ numbers are related to the well-known Catalan numbers $C_{n}$ (for
more information see Petersen (2015)). In this case, the Narayana numbers
$N_{n}^{m}$ divide combinatorial sets defined by the Catalan numbers $C_{n}$
into subsets in accordance with some additional parameter $m$. For example,
one combinatorial interpretation of the Catalan numbers is the number of Dyck
$n$-paths, then the Narayana numbers show the number of Dyck $n$-paths with
$m$ peaks (see Kruchinin et al. (2022)).

The following bivariate generating function of the Narayana numbers
$N_{n}^{m}$ is known from OEIS Foundation Inc. (2022):

$N(x,y)=\sum\limits_{n>0}{\sum\limits_{m>0}{N_{n}^{m}x^{n}y^{m}}}=\frac{1-x-xy-\sqrt{(1-x-xy)^{2}-4x^{2}y}}{2x},$

and for $n\geq m>0$ we have

$N_{n}^{m}=\frac{1}{n}\binom{n}{m-1}\binom{n}{m}.$

Hence, we can show the relationship of the numbers $S_{n}^{m}$ with the
Narayana numbers $N_{n}^{m}$ through the relationship of their generating
functions by

$S(x,y)=1+\frac{N(x,xy)}{xy}$

or their explicit formulas by

$S_{n}^{m}=N_{n-m}^{m+1},$

where $n>m\geq 0$ and $S_{0}^{0}=1$.

###  2.3 AND/OR tree

Since Equation (1) belongs to the algebra $\\{\mathbb{N},+,\times,R\\}$, we
can construct the AND/OR tree structure for $S_{n}^{m}$, which is presented in
Figure 2. In this AND/OR tree:

  * •

each addition in Equation (1) is represented as an OR node where all terms are
represented as its sons;

  * •

each multiplication in Equation (1) is represented as an AND node where all
factors are represented as its sons;

  * •

each recursive operation in Equation (1) is represented as recursion in the
AND/OR tree (denoted by a node with a triangle).

Figure 2: An AND/OR tree for $S_{n}^{m}$.

For this AND/OR tree structure, there are the following initial conditions:

  * •

if we get a node labeled $S_{n}^{0}$, then it is a leaf node in the AND/OR
tree where $n\geq 0$

  * •

if we get a node labeled $S_{n}^{m}$ where $m\geq\frac{n}{2}$, then we need to
remove this node and also remove the subtree containing this node until the
closest son of an OR node.

Figure 3 presents an example of an AND/OR tree structure for $S_{n}^{m}$ where
$n=6$ and $m=2$. The total number of its variants is equal to $S_{6}^{2}=6$.

Figure 3: The AND/OR tree for $S_{6}^{2}$.

###### Theorem 1

There is a bijection between the set of RNA secondary structures of length $n$
with $m$ base-pairs and the set of variants of the AND/OR tree for
$S_{n}^{m}$.

######  Proof 2.2.

Let represent RNA secondary structures of length $n$ with $m$ base-pairs as
correct bracket sequences of length $n$ with $m$ pairs of brackets with the
possibility of adding spaces. The total number of such RNA secondary
structures is equal to $S_{n}^{m}$. The total number of variants of the AND/OR
tree for $S_{n}^{m}$ presented in Figure 2 is also equal to $S_{n}^{m}$.
Hence, a bijection between the set of RNA secondary structures of length $n$
with $m$ base-pairs and the set of variants of the AND/OR tree for $S_{n}^{m}$
is defined by the following rules:

  * •

each selected left son of the OR node labeled $S_{n}^{m}$ determines the
addition of a space to the sequence $a=(a_{1},\ldots,a_{n})$: the sequence $a$
is represented as $a=\texttt{concat}$((’*’),$b$), i.e. the first element of
the sequence is the added space ($a_{1}=$ ’*’) and the remaining part of the
sequence denoted by $b=(a_{2},\ldots,a_{n})$ is a correct bracket sequence of
length $n-1$ with $m$ pairs of brackets (it corresponds to the subtree of the
node labeled $S_{n-1}^{m}$);

  * •

each selected right son of the OR node labeled $S_{n}^{m}$ determines the
addition of a pair of brackets to the sequence $a=(a_{1},\ldots,a_{n})$: the
sequence $a$ is represented as $a=\texttt{concat}$((’(’),$b$,(’)’),$c$), where
$b=(a_{2},\ldots,a_{n-j-1})$ is a correct bracket sequence of length $n-2-j$
with $m-1-i$ pairs of brackets (it corresponds to the subtree of the node
labeled $S_{n-2-j}^{m-1-i}$) and $c=(a_{n-j+1},\ldots,a_{n})$ is a correct
bracket sequence of length $j$ with $i$ pairs of brackets (it corresponds to
the subtree of the node labeled $S_{j}^{i}$), and the values of $i$ and $j$
are determined by the selected sons in the AND/OR tree.

  * •

each leaf node labeled $S_{n}^{0}$ determines the addition of $n$ spaces to
the sequence $a=(a_{1},\ldots,a_{n})$.

In the obtained rules for the bijection, the function concat denotes merging
sequences, i.e. for two given sequences $a=(a_{1},\ldots,a_{n})$ and
$b=(b_{1},\ldots,b_{m})$ we get the following sequence:

$c=\texttt{concat}(a,b)=(a_{1},\ldots,a_{n},b_{1},\ldots,b_{m})=(c_{1},\ldots,c_{n+m}).$

For a compact representation, a variant of the AND/OR tree for $S_{n}^{m}$
will be encoded by a sequence $v=(v_{1},v_{2})$, where:

  * •

an empty sequence $v=(\,)$ corresponds to the selection of a leaf node labeled
$S_{n}^{0}$;

  * •

$v_{1}$ corresponds to the selection of the left son ($v_{1}=0$) or the right
son ($v_{1}=1$) of the OR node labeled $S_{n}^{m}$;

  * •

if $v_{1}=0$, then $v_{2}$ corresponds to the variant of the subtree of the
node labeled $S_{n-1}^{m}$;

  * •

if $v_{1}=1$, then $v_{2}$ is encoded by a sequence $v_{2}=(I,J,vl,vr)$,
where:

    * –

$I$ corresponds to the selected value of $i$ in the AND/OR tree;

    * –

$J$ corresponds to the selected value of $j$ in the AND/OR tree;

    * –

$vl$ corresponds to the variant of the subtree of the node labeled
$S_{n-2-j}^{m-1-i}$;

    * –

$vr$ corresponds to the variant of the subtree of the node labeled
$S_{j}^{i}$.

Algorithm 1 and Algorithm 2 present the developed rules for the bijection
between the set of RNA secondary structures and the set of variants of the
corresponding AND/OR tree. In these algorithms, $()$ denotes an empty
sequence.

RNASSToVariant ($a=(a_{1},\ldots,a_{n})$, $n$, $m$)

begin

if  _$m=0$_ then return _$()$_

if  _$a_{1}=$ ’*’_ then $v:=(0,$ RNASSToVariant ($(a_{2},\ldots,a_{n})$,
$n-1$, $m$)$)$

else

$s_{1}:=0$

$s_{2}:=0$

for  _$i:=1$ to $n$_ do

if  _$a_{i}=$ ’(’_ then $s_{1}:=s_{1}+1$

if  _$a_{i}=$ ’)’_ then $s_{2}:=s_{2}+1$

if  _$s_{1}=s_{2}$_ then

$I:=m-s_{1}$

$J:=n-i$

$b:=(a_{2},\ldots,a_{n-J-1})$

$c:=(a_{n-J+1},\ldots,a_{n})$

break

end if

end for

$vl:=$ RNASSToVariant ($b$, $n-2-J$, $m-1-I$)

$vr:=$ RNASSToVariant ($c$, $J$, $I$)

$v:=(1,(I,J,vl,vr))$

end if

return _$v$_

end

Algorithm 1 An algorithm for converting a correct bracket sequence with the
possibility of adding spaces into a variant of the AND/OR tree for
$S_{n}^{m}$.

VariantToRNASS ($v=(v_{1},v_{2})$, $n$, $m$)

begin

if  _$m=0$_ then return _$(a_{1},\ldots,a_{n}):=($ ’*’$,\ldots,$’*’$)$_

if  _$v_{1}:=0$_ then $a:=$ concat ($($’*’$)$, VariantToRNASS ($v_{2}$, $n-1$,
$m$))

else

$(I,J,vl,vr):=v_{2}$

$b:=$ VariantToRNASS ($vl$, $n-2-J$, $m-1-I$)

$c:=$ VariantToRNASS ($vr$, $J$, $I$)

$a:=$ concat ($($’(’$)$, $b$, $($’)’$)$, $c$)

end if

return _$a$_

end

Algorithm 2 An algorithm for converting a variant of the AND/OR tree for
$S_{n}^{m}$ into a correct bracket sequence with the possibility of adding
spaces.

Let us consider the computational complexity of these algorithms in more
detail:

  * •

In Algorithm 1, there are $n-2m$ recursive calls where $a_{1}=$ ’*’ (each such
recursive call requires making a new sequence of length $n-1$) and there are
$m$ recursive calls where $a_{1}\neq$ ’*’ (each such recursive call requires
making calculations maximum $n$ times in the for-loop and making new sequences
of maximum length $n-2$). Hence, assuming algebraic operations with numbers in
$O(1)$, Algorithm 1 has polynomial time complexity $O(n(n-m)+nm)$;

  * •

In Algorithm 2, there are maximum $n-2m$ recursive calls where $v_{1}=0$ (each
such recursive call requires making a new sequence of length $n$) and there
are $m$ recursive calls where $v_{1}=1$ (each such recursive call requires
making a new sequence of length $n$). Hence, assuming algebraic operations
with numbers in $O(1)$, Algorithm 2 has polynomial time complexity
$O(n(n-m)+nm)$.

Thus, both presented algorithms have polynomial time complexity
$O(n(n-m)+nm)$. Also, taking into account the constraint $m<\frac{n}{2}$, we
get $O(n^{2})$.

The obtained bijection make it possible to apply the combinatorial generation
method mentioned above to the development algorithms for ranking and unranking
the variants of the constructed AND/OR tree structure for $S_{n}^{m}$.
Combining these combinatorial generation algorithms with the obtained
bijection, we will get algorithms for ranking and unranking the considered
combinatorial set of RNA secondary structures.

###  2.4 Ranking and unranking algorithms

Applying the general approach described in Shablya et al. (2020), we develop
algorithms for ranking and unranking the variants of the constructed AND/OR
tree structure. Algorithm 3 presents the required actions for ranking the
variants of the AND/OR tree for $S_{n}^{m}$, Algorithm 4 presents the required
actions for unranking its variants.

Let us consider the computational complexity of these algorithms in more
detail:

  * •

In Algorithm 3, there are maximum $n-2m$ recursive calls where $v_{1}=0$ (each
such recursive call requires making only one assignment) and there are $m$
recursive calls where $v_{1}=1$ (each such recursive call requires making
calculations of $S_{j}^{i}$ maximum $m$ times for $i$ and $n-2m$ times for
$j$). Applying Equation (2), the calculation of $S_{n}^{m}$ has polynomial
time complexity $O(m)$ for $m<n-2m$ and $O(n-m)$ for $m>n-2m$. Hence, assuming
algebraic operations with numbers in $O(1)$, Algorithm 3 has polynomial time
complexity $O(m^{2}(n-m))$ for $m<n-2m$ and $O(m(n-m)^{2})$ for $m>n-2m$.

  * •

In Algorithm 4, there are maximum $n-2m$ recursive calls where $r<S_{n-1}^{m}$
(each such recursive call requires making only one assignment) and there are
$m$ recursive calls where $r\geq S_{n-1}^{m}$ (each such recursive call
requires making calculations of $S_{J}^{I}$ maximum $m$ times for $I$ and
$n-2m$ times for $J$). Hence, assuming algebraic operations with numbers in
$O(1)$, Algorithm 4 has polynomial time complexity $O(m^{2}(n-m))$ for
$m<n-2m$ and $O(m(n-m)^{2})$ for $m>n-2m$.

Taking into account the constraint $m<\frac{n}{2}$, we get polynomial time
complexity $O(n^{3})$ for both presented algorithms. For comparison, the
ranking and unranking algorithms of Seyedi-Tabari et al. (2010) have better
polynomial time complexity $O(n^{2})$, however they require a preprocessing
step that has exponential time and space complexity. Since the developed
algorithms do not require preprocessing steps, they have better overall time
and space complexity. In addition, the developed algorithms can be
supplemented with a preprocessing step with the calculation of all required
values of $S_{n}^{m}$. As a result, we get the preprocessing step with
polynomial time complexity $O(n^{2})$ and the new ranking and unranking
algorithms with the same polynomial time complexity $O(n^{2})$.

RankVariant_RNASS ($v=(v_{1},v_{2})$, $n$, $m$)

begin

if  _$m=0$_ then return _$0$_

if  _$v_{1}=0$_ then $r:=$ RankVariant_RNASS ($v_{2}$, $n-1$, $m$)

else

$(I,J,vl,vr):=v_{2}$

$l_{1}:=$ RankVariant_RNASS ($vl$, $n-2-J$, $m-1-I$)

$l_{2}:=$ RankVariant_RNASS ($vr$, $J$, $I$)

$r:=l_{1}+S_{n-2-J}^{m-1-I}l_{2}+S_{n-1}^{m}+\sum\limits_{i=0}^{I-1}{\sum\limits_{j=2i}^{n-2(m-i)-1}{S_{n-2-j}^{m-1-i}S_{j}^{i}}}+\sum\limits_{j=2I}^{J-1}{S_{n-2-j}^{m-1-I}S_{j}^{I}}$

end if

return _$r$_

end

Algorithm 3 An algorithm for ranking a variant of the AND/OR tree for
$S_{n}^{m}$.

UnrankVariant_RNASS ($r$, $n$, $m$)

begin

if  _$m=0$_ then return _$()$_

if  _$r <S_{n-1}^{m}$_ then $v:=(0,$ UnrankVariant_RNASS ($r$, $n-1$, $m$)$)$

else

$r:=r-S_{n-1}^{m}$

$sum:=0$

$I:=0$

$J:=0$

while  _$sum+S_{n-2-J}^{m-1-I}S_{J}^{I}\leq r$_ do

$sum:=sum+S_{n-2-J}^{m-1-I}S_{J}^{I}$

if  _$J\leq n-2(m-I)-1$_ then $J:=J+1$

else

$I:=I+1$

$J:=2I$

end if

end while

$r:=r-sum$

$l_{1}:=r\mod{S_{n-2-J}^{m-1-I}}$

$l_{2}:=\left\lfloor{r}\,/\,{S_{n-2-J}^{m-1-I}}\right\rfloor$

$vl:=$ UnrankVariant_RNASS ($l_{1}$, $n-2-J$, $m-1-I$)

$vr:=$ UnrankVariant_RNASS ($l_{2}$, $J$, $I$)

$v:=(1,(I,J,vl,vr))$

end if

return _$v$_

end

Algorithm 4 An algorithm for unranking a variant of the AND/OR tree for
$S_{n}^{m}$.

Combining Algorithm 1 with Algorithm 3, we get the algorithm for ranking the
combinatorial set of RNA secondary structures of length $n$ with $m$ base-
pairs. The combination of Algorithm 4 with Algorithm 2 performs the inverse
operation (unranking). Table 2 and Table 3 present examples of ranking the
combinatorial set of RNA secondary structures of length $n$ with $m$ base-
pairs (that are represented as correct bracket sequences of length $n$ with
$m$ pairs of brackets with the possibility of adding spaces).

Table 2: An example of ranking the set of RNA secondary structures of length $n=6$ with $m=2$ base-pairs. Correct bracket sequence |  Variant of AND/OR tree | Rank
---|---|---
$*\,(\,(\,*\,)\,)$ |  $(0,(1,(0,0,(1,(0,0,(),())),())))$ | $0$
$(\,*\,(\,*\,)\,)$ |  $(1,(0,0,(0,(1,(0,0,(),()))),()))$ | $1$
$(\,(\,*\,*\,)\,)$ |  $(1,(0,0,(1,(0,0,(),())),()))$ | $2$
$(\,(\,*\,)\,*\,)$ |  $(1,(0,0,(1,(0,1,(),())),()))$ | $3$
$(\,(\,*\,)\,)\,*$ |  $(1,(0,1,(1,(0,0,(),())),()))$ | $4$
$(\,*\,)\,(\,*\,)$ |  $(1,(1,3,(),(1,(0,0,(),()))))$ | $5$

Table 3: An example of ranking the set of RNA secondary structures of length $n=8$ with $m=3$ base-pairs. Correct bracket sequence |  Variant of AND/OR tree | Rank
---|---|---
$*\,(\,(\,(\,*\,)\,)\,)$ |  $(0,(1,(0,0,(1,(0,0,(1,(0,0,(),())),())),())))$ | $0$
$(\,*\,(\,(\,*\,)\,)\,)$ |  $(1,(0,0,(0,(1,(0,0,(1,(0,0,(),())),()))),()))$ | $1$
$(\,(\,*\,(\,*\,)\,)\,)$ |  $(1,(0,0,(1,(0,0,(0,(1,(0,0,(),()))),())),()))$ | $2$
$(\,(\,(\,*\,*\,)\,)\,)$ |  $(1,(0,0,(1,(0,0,(1,(0,0,(),())),())),()))$ | $3$
$(\,(\,(\,*\,)\,*\,)\,)$ |  $(1,(0,0,(1,(0,0,(1,(0,1,(),())),())),()))$ | $4$
$(\,(\,(\,*\,)\,)\,*\,)$ |  $(1,(0,0,(1,(0,1,(1,(0,0,(),())),())),()))$ | $5$
$(\,(\,*\,)\,(\,*\,)\,)$ |  $(1,(0,0,(1,(1,3,(),(1,(0,0,(),())))),()))$ | $6$
$(\,(\,(\,*\,)\,)\,)\,*$ |  $(1,(0,1,(1,(0,0,(1,(0,0,(),())),())),()))$ | $7$
$(\,(\,*\,)\,)\,(\,*\,)$ |  $(1,(1,3,(1,(0,0,(),())),(1,(0,0,(),()))))$ | $8$
$(\,*\,)\,(\,(\,*\,)\,)$ |  $(1,(2,5,(),(1,(0,0,(1,(0,0,(),())),()))))$ | $9$

##  3 Conclusions

In this paper, we study the combinatorial set of RNA secondary structures of
length $n$ with $m$ base-pairs. For a compact representation, an RNA secondary
structure is encoded by a correct bracket sequence with the possibility of
adding spaces. For this combinatorial set, we construct an AND/OR tree
structure, find a bijection between the combinatorial set and the set of
variants of the AND/OR tree, and develop algorithms for ranking and unranking
the variants of the AND/OR tree.

The developed ranking and unranking algorithms have polynomial time complexity
$O(m^{2}(n-m))$ for $m<n-2m$ and $O(m(n-m)^{2})$ for $m>n-2m$. The conducted
computational experiments also confirmed the derived time complexity of the
algorithms.

The existing algorithms of Seyedi-Tabari et al. for ranking and unranking the
considered combinatorial set of RNA secondary structures of length $n$ with
$m$ base-pairs also have polynomial time complexity. However, the operation of
these ranking and unranking algorithms requires a preprocessing step that has
exponential time and space complexity. In contrast to the algorithms of
Seyedi-Tabari et al., the new algorithms do not require preprocessing steps.
Thus, the developed ranking and unranking algorithms have better time
complexity (the resulting complexity is polynomial).

## References

  * Flajolet et al. (1994) P. Flajolet, P. Zimmerman, and B. Cutsem.  A calculus for the random generation of combinatorial structures.  _Theoretical Computer Science_ , 132(1–2):1–35, 1994.  10.1016/0304-3975(94)90226-7.
  * Knuth (2011) D. E. Knuth.  _The Art of Computer Programming, Volume 4A: Combinatorial Algorithms, Part 1_.  Addison-Wesley Professional, USA, 2011.
  * Kreher and Stinson (1999) D. L. Kreher and D. R. Stinson.  _Combinatorial Algorithms: Generation, Enumeration, and Search_.  CRC Press, USA, 1999.
  * Kruchinin et al. (2022) D. Kruchinin, V. Kruchinin, and Y. Shablya.  On some properties of generalized Narayana numbers.  _Quaestiones Mathematicae_ , 45(12):1949–1963, 2022.  10.2989/16073606.2021.1980448.
  * Mohammadi and Nowzari-Dalini (2018) S. Mohammadi and A. Nowzari-Dalini.  A parallel algorithm for generation of RNA secondary structures with length $n$ and $k$ base-pairs.  _Iran Journal of Computer Science_ , 1:11–17, 2018.  10.1007/s42044-017-0001-2.
  * Nebel et al. (2011) M. E. Nebel, A. Scheid, and F. Weinberg.  Random generation of RNA secondary structures according to native distributions.  _Algorithms for Molecular Biology_ , 6(1):24, 2011\.  10.1186/1748-7188-6-24.
  * OEIS Foundation Inc. (2022) OEIS Foundation Inc.  The on-line encyclopedia of integer sequences, 2022.  URL https://page.math.tu-berlin.de/~felsner/SemWS17-18/Ruskey-Comb-Gen.pdf.
  * Oluoch et al. (2019) I. K. Oluoch, A. Akalin, Y. Vural, and Y. Canbay.  A review on RNA secondary structure prediction algorithms.  In _IEEE International Congress on Big Data, Deep Learning and Fighting Cyber Terrorism (IBIGDELFT 2018)_ , pages 18–23, 2019.  10.1109/ibigdelft.2018.8625347.
  * Petersen (2015) T. K. Petersen.  _Eulerian Numbers_.  Birkhauser Advanced Texts Basler Lehrbucher, USA, 2015.  10.1007/978-1-4939-3091-3.
  * Picardi (2015) E. Picardi.  _RNA Bioinformatics_.  Humana Press, USA, 2015.  10.1007/978-1-4939-2291-8.
  * Reidys (2011) C. Reidys.  _Combinatorial Computational Biology of RNA: Pseudoknots and Neutral Networks_.  Springer New York, USA, 2011.  10.1007/978-0-387-76731-4.
  * Ruskey (2003) F. Ruskey.  Combinatorial generation. working version (1j-csc 425/520), 2003.  URL https://page.math.tu-berlin.de/~felsner/SemWS17-18/Ruskey-Comb-Gen.pdf.
  * Ryabko (1998) B. Y. Ryabko.  Fast enumeration of combinatorial objects.  _Discrete Mathematics and Applications_ , 8(2):163–182, 1998.  10.4213/dm426.
  * Saenger (1984) W. Saenger.  _Principles of Nucleic Acid Structure_.  Springer New York, USA, 1984.  10.1007/978-1-4612-5190-3.
  * Sato and Kato (2022) K. Sato and Y. Kato.  Prediction of RNA secondary structure including pseudoknots for long sequences.  _Briefings in Bioinformatics_ , 23(1):bbab395, 2022.  10.1093/bib/bbab395.
  * Sato et al. (2021) K. Sato, M. Akiyama, and Y. Sakakibara.  RNA secondary structure prediction using deep learning with thermodynamic integration.  _Nature Communications_ , 12:941, 2021.  10.1038/s41467-021-21194-4.
  * Schmitt and Waterman (1994) W. R. Schmitt and M. S. Waterman.  Linear trees and RNA secondary structure.  _Discrete Applied Mathematics_ , 51(3):317–323, 1994.  10.1016/0166-218x(92)00038-n.
  * Seyedi-Tabari et al. (2010) E. Seyedi-Tabari, H. Ahrabian, and A. Nowzari-Dalini.  A new algorithm for generation of different types of RNA.  _International Journal of Computer Mathematics_ , 87(6):1197–1207, 2010.  10.1080/00207160802140049.
  * Shablya et al. (2020) Y. Shablya, D. Kruchinin, and V. Kruchinin.  Method for developing combinatorial generation algorithms based on AND/OR trees and its application.  _Mathematics_ , 8(6):962, 2020.  10.3390/math8060962.
  * Viennot and Vauchaussade de Chaumont (1985) G. Viennot and M. Vauchaussade de Chaumont.  Enumeration of RNA secondary structures by complexity.  _Mathematics in Biology and Medicine_ , pages 360–365, 1985.  10.1007/978-3-642-93287-8_50.
  * Waterman (2018) M. S. Waterman.  _Introduction to Computational Biology_.  CRC Press, USA, 2018.
  * Zuker and Sankoff (1984) M. Zuker and D. Sankoff.  RNA secondary structures and their prediction.  _Bulletin of Mathematical Biology_ , 46(4):591–621, 1984.  10.1007/bf02459506.


# Semiclassical propagation of coherent states and wave packets: hidden
saddles

Huichao Wang and Steven Tomsovic  Department of Physics and Astronomy,
Washington State University, Pullman, WA. USA 99164-2814

###### Abstract

Semiclassical methods are extremely important in the subjects of wave packet
and coherent state dynamics. Unfortunately, these essentially saddle point
approximations are considered nearly impossible to carry out in detail for
systems with multiple degrees of freedom due to the difficulties of solving
the resulting two-point boundary value problems. However, recent developments
have extended the applicability to a broader range of systems and
circumstances. The most important advances are first to generate a set of real
reference trajectories using appropriately reduced dimensional spaces of
initial conditions, and second to feed that set into a Newton-Raphson search
scheme to locate the exposed complex saddle trajectories. The arguments for
this approach were based mostly on intuition and numerical verification. In
this paper, the methods are put on a firmer theoretical foundation and then
extended to incorporate saddles hidden from Newton-Raphson searches initiated
with real trajectories. This hidden class of saddles is relevant to tunneling-
type processes, but a hidden saddle can sometimes contribute just as much as
or more than an exposed one. The distinctions between hidden and exposed
saddles clarifies the interpretation of what constitutes tunneling for wave
packets and coherent states in the time domain.

##  I Introduction

The evolutions of Glauber coherent states for bosonic many-body systems and
the mathematically nearly-identically related Gaussian wave packets arise in a
multitude of physical contexts. Examples abound in quantum optics [1, 2], far
out-of-equilibrium dynamics in bosonic many-body systems [3, 4], molecular
spectroscopy [5, 6], femto-chemistry [7], and attosecond physics [8].
Theoretical work encompasses, for example, coherent state representations of
path integrals [9] and in the context of many bosons or the short wavelength
limit, the semiclassical approximation [9, 10, 11]. This approximation is
fundamental to studies of quantum-classical correspondence [12], and pre- and
post-Ehrenfest-time-scale dynamics [13, 14, 15] as well.

A complete semiclassical approximation for coherent state dynamics may be
obtained by the saddle point approximation applied to coherent state path
integrals [9, 11]. In the context of wave packets the essentially identical
approximation is known as generalized Gaussian wave packet dynamics (GGWPD)
and it has been proven to be equivalent to a complexified version of time-
dependent WKB theory [10]. In either case, implementation in practice leads to
considerable technical difficulties for any system possessing nonlinear
dynamics.

To begin with, the classical trajectories that define the saddles’ properties
are solutions of a two-point boundary value problem, which is highly
nontrivial for nonlinear dynamical systems. If there are many degrees of
freedom, this problem may be effectively impossible to solve. It appears that
for any propagation time $t\neq 0$, there exists an infinity of saddles; see
Figure 1 ahead. At the shortest evolution time scale, for minimum uncertainty
wave packets (coherent states) only one saddle is physically relevant. For a
chaotic dynamical system, this time scale is extremely short and has a
logarithmic dependence on $\hbar$. The time scale is much longer for an
integrable system, scaling algebraically in $\hbar$. Nevertheless, in either
situation this defines an Ehrenfest time regime [16]. Parenthetically,
physically relevant means: i) the saddle is on the correct side of the Stokes
lines; and ii) its contribution to the wave function at the relevant position
is significant enough to be larger than any of the errors due to the saddle
point approximations from other saddles contributing to the same position. As
time increases beyond the Ehrenfest time, the number of physically relevant
saddles increases, necessarily creating quantum interferences, but the number
is always finite at some fixed propagation time. One might say that a set of
measure zero saddles from the infinite set must be selected as the physically
relevant saddles. Furthermore, these trajectories necessarily involve
complexified position and momenta. This analytic continuation of real
classical dynamics has many nontrivial features, such as branch cuts
associated with singular runaway trajectories [10], and prefactor
singularities that pose new problems for determining the phase (see the
companion paper to this one ahead [17]).

Considerable progress has been made in developing a practical method of
identifying physically relevant saddles directly without encountering any of
the irrelevant saddles [18, 19, 20, 21]. The basic idea builds on earlier
techniques of identifying a single real reference trajectory for each
classical transport pathway that exists in the nonlinear dynamics of the
system [22, 23, 24, 25]. Then using a Newton-Raphson algorithm, a unique
saddle point is identified for each classically allowed transport pathway and
it accounts for that pathway’s contributions to the dynamics. The techniques
have been developed using Gaussian wave packets, but it applies equally well
to bosonic coherent states expressed in quadratures; see [15]. In many cases
the method can be extended to systems with many degrees of freedom by
identifying and neglecting directions in phase space that do not lead to
diverging initial conditions [21].

Using these techniques, a complex saddle is identified with each classically
allowed transport pathway. Let us dub these “exposed saddles”. Amongst the
infinity of exposed saddles, almost all contribute too little to be concerned
with, but it is straightforward to restrict the search for real transport
pathways that have sufficient amplitude to be relevant. Nevertheless, there is
still a great deal more to be done. These works were justified intuitively and
left open the question of how to identify physically relevant saddles
associated with classically non-allowed transport pathways. In contrast to
exposed saddles, consider these “hidden saddles”, in part because they are not
directly discoverable using real trajectory input into a Newton-Raphson search
scheme. The focus of this work is to examine these techniques in greater
detail, add additional justifications where possible, and develop a method to
locate the physically relevant hidden saddles. This extends the previous
practical techniques [20, 15, 21] to incorporate tunneling-like phenomena and
provides an interpretation of exactly what constitutes tunneling in the time
domain for wave packet and coherent state propagation.

The structure of the paper is as follows, Sect. II defines the critical
quantities, sets notations, introduces the purely quartic oscillator for
simple illustrations, and discusses the background of various semiclassical
approximation methods and their interrelationships. This is followed in Sect.
III by a discussion of the justification of the existing methods, and develops
a technique for identifying hidden saddles. The paper concludes with a summary
of the work and possible future lines of related research.

##  II Background

In order to elucidate the method for implementing a search for hidden saddles
and justify certain procedures, it is helpful to start with some background on
the complete semiclassical approximation along with certain partial or
imperfect versions of the semiclassical approximation. As the mathematical
results and manipulations can be made to appear essentially identical through
the application of quadratures, it is unnecessary to treat evolution of
Glauber coherent states and wave packets separately. The discussion is
presented in the language of wave packets, but it is understood that all the
results carry over to coherent states. As a final note, for simplicity the
discussion and equations given here are reduced to their single degree of
freedom forms. Although, the interest is in multi-degree-of-freedom systems,
it is much easier to illuminate the basic ideas with a simple example and the
equations reduced to their one degree of freedom forms. All the necessary
multi-degree-of-freedom equations can be found elsewhere, for example in [20,
21], and the additional ideas presented here extend to the many degrees of
freedom case.

###  II.1 Coherent states

A Glauber coherent state describing a bosonic many-body system takes the
normalized form

$|z\rangle=\exp\left(-\frac{\left|z\right|^{2}}{2}+z\hat{a}^{\dagger}\right)|0\rangle$
(1)

with the properties that

$\hat{a}|z\rangle=z|z\rangle\
,\quad\left[\hat{a},\hat{a}^{\dagger}\right]=\mathbb{1}\ ,\quad\langle
z|\hat{a}^{\dagger}\hat{a}|z\rangle=|z|^{2}\ .$ (2)

Hence, a coherent state’s mean number of bosons $n=\langle
z|\hat{n}|z\rangle=|z|^{2}$. Its evolved form $|z(t)\rangle=U(t)|z\rangle$
does not remain a coherent state, but the overlap with the bra vector version
of another coherent state $\langle z^{\prime}|$ can be viewed as a matrix
element of a coherent state path integral, $\langle z^{\prime}|U(t)|z\rangle$.
In the language of wave packets, this is often termed a correlation function
or transport coefficient.

Introducing the quadrature (Hermitian) operators normalized as follows

$\displaystyle\hat{p}$ $\displaystyle=$ $\displaystyle
i\sqrt{\frac{\tilde{\hbar}}{2}}\left(\hat{a}^{\dagger}-\hat{a}\right)$
$\displaystyle\hat{x}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\tilde{\hbar}}{2}}\left(\hat{a}^{\dagger}+\hat{a}\right)$

the operator pair behaves mathematically exactly like the quantized version of
a canonical pair of momentum and position variables, even though they cannot
be interpreted as some boson’s momentum and position. Setting $\langle
z|\hat{x}|z\rangle=q_{\alpha}$ and $\langle z|\hat{p}|z\rangle=p_{\alpha}$,
then $z=(q_{\alpha}+ip_{\alpha})/\sqrt{2\tilde{\hbar}}$. The $\tilde{\hbar}$
constant is written to look like the usual $\hbar$ to emphasize the analogy,
but is rather usefully set to the inverse number of bosons making the
$n\rightarrow\infty$ limit equivalent to the $\tilde{\hbar}\rightarrow 0$
limit. With the quadrature operator $\hat{x}$’s eigenbasis representation
$\hat{x}|x\rangle=x|x\rangle$, the coherent state takes the form

$\displaystyle\langle
x|z\rangle=\left(\frac{1}{\pi\tilde{\hbar}}\right)^{1/4}$
$\displaystyle\exp\left[-\frac{\left(\
x-q_{\alpha}\right)^{2}}{2\tilde{\hbar}}+\frac{ip_{\alpha}}{\tilde{\hbar}}\left(x-q_{\alpha}\right)\right.$
$\displaystyle\qquad\qquad\qquad\left.+\frac{i}{2\tilde{\hbar}}p_{\alpha}q_{\alpha}\right]\
.$ (4)

Thus, an appropriate projection of the coherent state results in a Gaussian
wave packet form [1] with parameters which can be straightforwardly mapped
onto those of a wave packet. Compared with Eq. (II.2) just ahead, the forms
are equivalent with the replacements $\tilde{\hbar}\rightarrow\hbar$ and
$1\rightarrow b_{\alpha}$. One exposition of the parameter mapping is given in
Appendix A of [21]. There are sufficient wave packet parameters in Eq. (II.2)
to squeeze the coherent states and rotate the quadrature operators. The focus
now turns to Gaussian wave packets.

###  II.2 Gaussian wave packets

A normalized Gaussian wave packet may be parameterized as follows:

$\displaystyle\phi_{\alpha}(x)=$
$\displaystyle\exp\left[-\frac{b_{\alpha}}{2\hbar}\left(\
x-q_{\alpha}\right)^{2}+\frac{ip_{\alpha}}{\hbar}\left(x-q_{\alpha}\right)+\frac{i}{2\hbar}p_{\alpha}q_{\alpha}\right]$
$\displaystyle\qquad\times\left[\frac{b_{\alpha}+b^{*}_{\alpha}}{2\pi\hbar}\right]^{1/4},$
(5)

where the subscript $\alpha$ is a label for the parameters that define the
particular wave packet, $x$ is the position variable for the quantum system,
and $(q,p)$ are the canonically conjugate position-momentum phase space
variables for the analogous classical system. The real centroid is given by
$(q_{\alpha},p_{\alpha})$ and the width by $b_{\alpha}$ (if $b_{\alpha}$ is
complex, the wave packet has a chirped phase dependence). All the calculations
in this paper, use

$(q_{\alpha},p_{\alpha})=(0,20),\>b_{\alpha}=32$ (6)

as a convenient set for illustration ahead. This form has the advantage that
$\hbar$ does not explicitly appear in the equations for the two-point boundary
value problem given ahead; i.e. the equations become purely classical. Note
that this wave packet form is given here so as to have the exact same global
phase convention as the coherent state of Eq. (1), but the global phase term
could just be dropped. It also leaves the overall shape of its Wigner
transform independent of $\hbar$, other than the volume (overall scale).
Similarly to a coherent state, for nonlinear dynamical systems the evolved
wave packet $\phi_{\alpha}(x;t)$ is generally not Gaussian. Note that a matrix
element of the coherent state path integral could be expressed using
quadratures analogously to

${\cal A}_{\beta\alpha}(t)=\int_{-\infty}^{\infty}{\rm d}x\
\phi^{*}_{\beta}(x)\phi_{\alpha}(x;t)$ (7)

for wave packets, and which could be thought of as a transport coefficient. If
$\beta=\alpha$, it would be a diagonal element or a return amplitude.

####  II.2.1 Lagrangian manifolds

Lagrangian manifolds play a central role in semiclassical approximations, i.e.
WKB theory [26]. They provide a very geometric picture of the application of
the saddle point approximation. Using the parameterization of Eq. (II.2), the
appropriate manifold for a Gaussian wave packet was identified in [10] as

$b_{\alpha}\left({\cal q}-q_{\alpha}\right)+i\left({\cal
p}-p_{\alpha}\right)=0,$ (8)

where $({\cal p},{\cal q})$ are chosen from the sets of all complex positions
and conjugate momenta satisfying this equation. This generates a complex line
in a two dimensional complex phase space (complex position and momentum),
somewhat akin to a plane embedded in a four dimensional space. Likewise for
the complex conjugate wave packet (dual version), the equation is

$b^{*}_{\alpha}\left({\cal q}-q_{\alpha}\right)-i\left({\cal
p}-p_{\alpha}\right)=0.$ (9)

For a position ket or bra vector, the Lagrangian manifold would be the value
of the position and the set of all complex momenta.

One consequence of the Lagrangian manifold being necessarily complex is that
it blurs the distinction between classically allowed and non-allowed
processes. In ordinary WKB, the tori used for constructing wave functions are
real Lagrangian manifolds, the intersections of two manifolds generate
stationary phase points (not saddles), and the exponentially decaying tails of
wave functions require propagation with imaginary time or some analytic
continuation, for example, the introduction of imaginary momenta [27]. Here it
is the hidden saddles which correspond to classically non-allowed processes.
Since both hidden and exposed saddles are linked to complex trajectories,
“real” versus “complex” is not the distinguishing factor. There is however the
intuitive notion of complex phase points close or near to being real; see [18]
and references therein. Since distance is not defined in phase spaces, a
priori, the concept of “near” is not well defined. Even worse, a saddle can
have an arbitrarily large imaginary component in its initial condition and
still be exposed. Part of our work here is to add some precision to this
concept, which ends up helping to categorize the saddles into exposed or
hidden groups. Just as the distinction of real and complex trajectories is
physically significant in the case of ordinary WKB theory, so is the
exposed/hidden saddle distinction important physically. It also alters
considerably how one must search for the saddles, which is treated ahead.

####  II.2.2 Two-point boundary value problem

The two-point boundary value problem that arises in the saddle point
approximation can be succinctly expressed using these manifolds. In essence,
the manifold corresponding to the ket vector is propagated to time $t$
according to the complexified classical dynamical equations of motion, and its
intersections with the manifold associated with bra vector give the
trajectories that define the saddle points’ properties. Letting $({\cal
p}_{0},{\cal q}_{0})$ represent the complex initial condition for a trajectory
and $({\cal p}_{t},{\cal q}_{t})$ the complex phase point that is generated by
evolving this initial condition to time $t$, the complete two-point boundary
value problem can be written as

$\displaystyle b_{\alpha}\left({\cal q}_{0}-q_{\alpha}\right)+i\left({\cal
p}_{0}-p_{\alpha}\right)$ $\displaystyle=$ $\displaystyle 0,$ $\displaystyle
b^{*}_{\beta}\left({\cal q}_{t}-q_{\beta}\right)-i\left({\cal
p}_{t}-p_{\beta}\right)$ $\displaystyle=$ $\displaystyle 0.$ (10)

Each solution of these equations is a potential saddle point in the theory for
transport coefficients. If instead, one is interested in the evolution of the
wave packet itself, the equations would be given by

$\displaystyle b_{\alpha}\left({\cal q}_{0}-q_{\alpha}\right)+i\left({\cal
p}_{0}-p_{\alpha}\right)$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle{\cal q}_{t}-x$ $\displaystyle=$ $\displaystyle 0,$ (11)

where $x$ is the position argument of the evolved wave packet (a momentum
representation is also possible, but not used here). For either boundary value
problem, nonlinear dynamical systems lead to an infinite number of solutions
to these equations. Almost all of them must be thrown away due to boundary
conditions or due to their contributions being vastly smaller than the errors
inherent in the semiclassical approximation. Excluding them from the search
algorithm by design is highly desirable and the basis of the works [20, 21].

###  II.3 Exposed and hidden saddles

The solutions of Eqs. (II.2.2) and (II.2.2) contain as a subset all the
necessary (i.e. physically relevant) saddle points of the corresponding
integrals. On the other hand, taken alone as an abstract set of equations
there is no indication at face value which of the infinity of solutions must
be thrown away, which must be kept, which can be related to real classical
transport processes, which cannot, etc… A priori, there is no interpretation
of tunneling processes in the coherent state/wave packet propagation as all
solutions are complex. From a mathematical perspective, it is neither trivial
to locate the solutions nor unpack how to separate out the relevant saddles
and determine what their physical interpretation is. For example, excluding
simple standard illustrations of Stokes phenomena, such as the Airy function
and the construction of its Stokes and anti-Stokes lines, it is highly
desirable to avoid those lines’ explicit construction and use, especially as a
function of many complex variables. The interest is in finding a direct,
practical approach to these challenges.

One procedure recently developed [20, 21] is to identify exposed saddles
directly using a Newton-Raphson scheme. Imagine to begin with, running all
possible real trajectory initial conditions, $(q_{0},p_{0})$, and letting them
serve as the initial input to the Newton-Raphson search algorithm. Excluding
cases in which convergence is impacted by singular trajectories [10] and
complex time paths are needed (see the “bald spot” problem [28]), certain
groupings of initial conditions lead to the exact same saddle point because
they are in the same basin of attraction. Thus, the real phase points can be
grouped into subregions (with negligible exceptions). All trajectories within
a single subregion lead to the same saddle point. A subregion’s trajectories
all shadow each other and behave similarly in a dynamical sense. Each
subregion thus represents a single classical transport pathway, and it is
necessary only to find a single reference trajectory to represent that pathway
as an input for the actual Newton-Raphson search. This implies that a complex
exposed saddle trajectory shadows the real trajectories associated with that
pathway. In particular, the stability matrix relation

$\left(\begin{array}[]{c}\delta\cal p_{t}\\\ \delta\cal
q_{t}\end{array}\right)=\left(\begin{array}[]{c}\bf{M_{11}}\\\
\bf{M_{21}}\end{array}\begin{array}[]{c}\bf{M_{12}}\\\
\bf{M_{22}}\end{array}\right)\left(\begin{array}[]{c}\delta\cal p_{0}\\\
\delta\cal q_{0}\end{array}\right)$ (12)

allowing for complex initial condition deviations from the real reference
trajectory would adequately describe the true exposed saddle trajectory. The
exposed saddle acquires the interpretation that it accounts for the
contributions of the classically allowed transport pathway to which the
reference trajectory belongs, and there is a one-to-one correspondence within
the collection of exposed saddles and transport pathways.

An immediate consequence of this shadowing relationship is that the exposed
saddle’s classical action would almost exclusively vary in its real part as a
function of time or other parameters leading to predominantly oscillatory
contributions, say for example, to the propagated wave function. Its imaginary
part would be positive, constrained to vary slowly, and account mainly for the
Gaussian weighting of the initial condition and/or final conditions in the
wave packets’ Wigner transform.

The above procedure necessarily leaves out the subset of physically relevant
saddles associated with tunneling and/or classically non-allowed processes,
which we classify as hidden saddles since they are hidden from the Newton-
Raphson search just described. One main purpose of this paper is to extend
that procedure so as to obtain these saddles as well; see Sect. III.2.
Evidently, these saddle trajectories behave rather differently, and cannot
shadow any of the real trajectories. One may anticipate that their classical
actions vary as a function of parameters more predominantly in their positive
imaginary parts rather than their real parts. In this way, they contribute
largely to exponentially decaying behaviors.

The complexification of the classical dynamics does have some less familiar
fundamental impacts. An example is given in [28] of saddles obscured by bald
spots. The basic idea is that singular complex trajectories create dynamical
branch cuts and to arrive at some saddle solutions, it is necessary to deform
a real time propagation path into some complex time domain path. Thus, one can
also anticipate the complication that there exists another class of saddles,
those behind bald spots, which require complex time integration paths to
follow. This class does not exist for the simple example described next nor
was it significant enough in any of the Bose-Hubbard model calculations of
[15, 21] to impact the accuracy. They will not be considered further here.

###  II.4 The purely quartic oscillator

The one degree of freedom purely quartic oscillator has a number of simple
features that makes it ideal for illustrating the ideas discussed in this
paper. Its analytically continued form [complex (${\cal p},{\cal q}$)] is
given by

$H({\cal p},{\cal q})=\frac{{\cal p}^{2}}{2}+\lambda{\cal q}^{4},$ (13)

where $\hbar=m=1$ and $\lambda=0.05$ are the values taken for all the
illustrations shown in this paper. The corresponding Schrödinger equation is
given by

$i\frac{\partial}{\partial t}\phi(x;t)=\left(-\frac{\partial^{2}}{2\partial
x^{2}}+\lambda x^{4}\right)\phi(x;t).$ (14)

Being a homogeneous Hamiltonian, the real classical dynamics of this system
leads to simple scaling relations amongst the quantities $(q,p;t;E)$. For
example, propagating an initial condition $(q_{E_{0}},p_{E_{0}})$ for a time
$t$ that belongs on the $E_{0}$ energy surface leads to a scaled trajectory
replica on any desired second energy surface $E$ according to:

$\displaystyle q_{E}(t/\gamma)$ $\displaystyle=\gamma q_{E_{0}}(t),$
$\displaystyle p_{E}(t/\gamma)$ $\displaystyle=\gamma^{2}p_{E_{0}}(t),$ (15)

where $\gamma=(E/E_{0})^{1/4}$. Classical actions, defined as $S=\int{\cal
L}{\rm d}t$, scale as $\gamma^{3}$, and the period of a closed orbit,
$\tau_{E_{0}}$, scales as $\tau_{E}=\tau_{E_{0}}/\gamma$. Therefore, if
$\gamma=n^{\pm 1}$, where $n$ is any positive integer, the two periodic orbits
on different energy surfaces would have periods that are an integer multiple
of each other. In the time that the longer period orbit closes once, the
shorter period orbit retraces itself $n$ times.

The potential generates nonlinear dynamical equations that lead to
sufficiently complicated behaviors for our purposes. As a first example, see
Fig. 1, which illustrates the placement of just the exposed saddle initial
conditions in a plane for the evolution of an initial wave packet evaluated at
$x=0$ and some fixed time $t$.

Figure 1: Initial conditions of exposed saddles satisfying Eq. (II.2.2) for a
wave packet defined by Eqs. (II.2) and (6), and propagated with the
Hamiltonian, Eq. (13). The complex momentum of a saddle’s initial conditions
follows from the complex position using Eq. (8). In the upper panel, the
propagation time is the period of the central trajectory with initial
condition $(q_{0},p_{0})=(0,20)$ and in the lower panel twice this period. The
density of exposed saddles increases with increasing time, but even in the
limit of small times, the total number of just the exposed saddles is
infinite. The only exception is for $t=0$.

There is no limiting domain within the Lagrangian manifold for the solutions,
and an infinity of solutions is implied over the manifold’s infinite domain.
In particular, note that these exposed saddles can have imaginary parts of
their initial conditions that extend out to $\pm$ infinity. What alters with
propagation time is the density of solutions. As time increases, double the
propagation time is shown in the lower panel, the density increases. At $t=0$,
there is effectively only one saddle left with all the other saddles pushed
out to infinity. Keep in mind that these exposed saddles are only a small
subset of all the saddles. There are solutions everywhere throughout the
complex ${\cal q}$ plane, not just approximately along a line, which is where
the exposed saddles find themselves.

###  II.5 Wigner transform

The Wigner transform of a wave packet (ket-like) generates a multivariate
Gaussian function, which can be thought of as a density of real classical
initial conditions that underlie a quantum wave packet and account for the
uncertainty principle. It can be used to help understand the partial
semiclassical approximation known as linearized wave packet dynamics [13], and
is essential for discussing the basis of an enhanced approximation known as an
off-center, real trajectory method [29, 23, 24, 30, 25].

The Wigner transform of a wave packet parametrized as in Eq. (II.2) is given
by

${\cal
W}(p,q)=\frac{1}{\pi\hbar}\exp\left[-\left(p-p_{\alpha},q-q_{\alpha}\right)\cdot\frac{{\bf
A}_{\alpha}}{\hbar}\cdot\left(\begin{array}[]{c}p-p_{\alpha}\\\
q-q_{\alpha}\end{array}\right)\right],$ (16)

where ${\bf A}_{\alpha}$ is

${\bf A}_{\alpha}=\left(\begin{array}[]{cc}1/c&d/c\\\
d/c&c+d^{2}/c\end{array}\right),\qquad{\rm Det}\left[{\bf
A}_{\alpha}\right]=1$ (17)

with the association

$b_{\alpha}=c+id.$ (18)

The $2\times 2$ dimensional matrix ${\bf A}_{\alpha}$ is real and symmetric.
If $b_{\alpha}$ is real, there are no correlations between $p$ and $q$ ($d$
vanishes); i.e. the wave packet is not chirped. The off-diagonal elements
(blocks in more degrees of freedom) of ${\bf A}_{\alpha}$ disappear.

###  II.6 Partial semiclassical theories

There are two important partial semiclassical approximations known as
linearized wave packet dynamics and off-center, real trajectory methods,
respectively, in order of increasing sophistication. The former linearizes the
dynamics completely, and the latter contains the nonlinear dynamical
information, which is much, much closer to GGWPD.

####  II.6.1 Linearized wave packet dynamics (LWPD)

The main idea involves solving the equations of motion for the parameters that
define an evolved Gaussian, including the center, width, and
phase/normalization factor in a global fashion. Its validity is constrained to
cases where the quadratically expanded potential around the central trajectory
is a good approximation to the dynamics within the phase space volume defined
by the Wigner transform of the evolving wave packet. As the wave packet
spreads out with increasing time, it soon extends well beyond the domain of
validity for the quadratic expansion, which leads to the breakdown of the
approximation.

The validity is limited to pre-Ehrenfest time scales. One method for
determining this time scale in practice is to consider two evolved densities.
The first is the initial wave packet’s Wigner transform evolved classically
(numerically), and the second is the Wigner transform of the LWPD evolved wave
packet. Let the overlap integral of these two densities at $t=0$ be normalized
to unity. Then LWPD remains an accurate approximation until a time such that
the overlap drops precipitously from one (this is the Ehrenfest time scale).
This is because LWPD has only linearized dynamical information about the
central trajectory whereas the density of initial conditions incorporates
nonlinear dynamical effects relative to the central trajectory. See Fig. 2 for
an illustration. The one and three standard deviation ($\sigma$) ellipses of
LWPD appear straight and are superposed on top of the $1$ and $3\sigma$ curves
of the initial state’s propagated Wigner transform. After one period of the
central trajectory’s motion some information about the wave packet central
region may still be given correctly by LWPD, but beyond the $1\sigma$ curve,
the dynamical information is no longer accurate. This illustrates the
breakdown of LWPD at short time scales for nonlinear dynamical systems.

Figure 2: Illustration of the limitations of LWPD. The vertically oriented
ellipses are contours of the initial wave packet’s Wigner transform using the
parameters of Eq. (6). This density is classically propagated for $1$ period
of the central trajectory’s motion ($\tau$), forming curvy foliations. The
LWPD produced wave packet generates stretched and rotated ellipses. Only phase
points within $1\sigma$ are still approximately linearly related to the
central trajectory by this propagation time.

####  II.6.2 Off-center, real trajectory method

The important idea underlying this method is that for any nonlinear dynamical
system over time a localized density of initial conditions disperses and
especially in a bounded system, the dispersed evolved density must fold over
or one might say create foliations. In some local region of phase space, these
foliations slice through (with increasing density as time increases), and each
one contains an infinite collection of nearly identically behaving
trajectories; i.e. a single classical transport pathway. See the study of the
production of whorls and tendrils in [12]. In other words, the stability
matrix of one member of the local set predicts rather well the behavior of its
neighbors within the foliation. On the other hand, each foliation has a
completely different dynamical character and represents a quite distinct
transport pathway. By identifying a member from each foliation, which
necessarily involves off-center initial conditions [those initial conditions
within the Wigner density, but not equal to $(q_{\alpha},p_{\alpha})$], and
following the essential prescription of LWPD, a wave is constructed whose
information is locally correct. Summing over the contributions of all the
foliations incorporates the nonlinear dynamical features. This technique
locally reconstructs the features of the evolving wave packet including
quantum interference between transport pathways, and its validity goes well
beyond the Ehrenfest time scale [29, 22, 31].

Figure 3: Illustration of distinct classical transport pathways. The
vertically oriented ellipses are contours of the initial wave packet’s Wigner
transform using the parameters of Eq. (6). The left panel is propagated
classically for $1$ period of the motion ($\tau$) and the right panel $3$
periods. They generate the distorted spirals. The insets show expanded views
of localized regions. In the $t=3\tau$ figure, the $5\sigma$ contour can be
divided into $9$ foliations (enumerated in later figures), a trajectory member
of which can be used to generate initial conditions for a Newton-Raphson
scheme to locate unique saddles. The blackened dots show the locations of the
real parts of the $9$ associated saddles’ complex final phase space points
contributing to $\phi(x=0,3\tau)$ (the initial conditions follow by back
propagating a time $3\tau$). An analogous statement applies to the $t=1\tau$
figure on the left.

The Wigner transformed initial density plays a more critical role for this
method since it is essentially being used to identify the foliations
themselves as a function of time. This is illustrated in Fig. 3 where each
foliation can be used to identify a single real reference trajectory about
which one can linearize the dynamics locally and that can be used to generate
a contribution to the quantity of interest. Furthermore, each of these off-
center, real reference trajectories can be used in a Newton-Raphson search to
locate a unique exposed saddle, which generates an important subset of the
saddles of GGWPD [20].

It also lends itself to a prescription for identifying the phase space
directions of fastest dispersion thereby indicating the directions that are
unnecessary for developing the foliations for the quantity of interest [21].
In this way, problems with large numbers of degrees of freedom can be reduced
to vastly fewer numbers. This effective dimensional reduction makes it
possible to treat fully more complicated dynamical systems. Even if many
degrees of freedom are active in this sense, they usually separate out in time
scales. Thus, one can successively add directions as time increases, but start
with just one or two. For many problems, being able to go to intermediate time
ranges beyond the Ehrenfest scale is sufficient [15].

###  II.7 Complications due to complexifying the classical dynamics

There are a number of complications arising if classical dynamics are
analytically continued to complex phase space variables. One of the most
significant considerations is the existence of singular orbits, i.e. orbits
that acquire infinite momenta in finite propagation times [10]. If time is
also allowed to be complex so that a path can be defined to avoid the singular
real time point, there may exist multi-sheeted functions describing the
dynamics [32, 33]. Numerically for the quartic oscillator, important dynamical
quantities such as position, momentum and classical action are not multi-
valued functions of time. This is consistent with it possessing the Painlevé
property in which solutions of a differential equation have no movable branch
points. In this sense, for the pure quartic oscillator the singularities are
isolated poles in a complex propagation time space.

On the other hand, these singularities must come in continuous families.
Consider a trajectory which becomes singular exactly at a propagation time
$t_{0}$. Differentially nearby will be other trajectories that become singular
at real propagation times of $t_{0}\pm\epsilon$. By following the changing
initial conditions for a singular trajectory continuously as the propagation
time varies from $t=0$ to $t=\infty$, they can be associated to a continuously
shifting set of initial conditions. This is illustrated for the quartic
oscillator in Fig. 4.

Figure 4: Contours of the final position $x$’s real part as a function of the
trajectory initial conditions on the initial wave packet’s Lagrangian
manifold. Real and imaginary parts of $\cal p_{0}$ follow using Eq. (8) with
the parameters of Eq. (6); a similar figure (Fig.3) was given in [10]. The
black dot near \raisebox{-0.9pt}{7}⃝ indicates the real initial condition of
the Wigner density centroid. The initial conditions of exposed saddles
associated with foliations \raisebox{-0.9pt}{6}⃝, \raisebox{-0.9pt}{7}⃝,
\raisebox{-0.9pt}{8}⃝ as a function of the propagated wave function’s position
$x$ appear as nearly vertical solid lines perpendicular to the contour lines
of $Re\,x$; they are located on $Im\,x=0$ contours which are not shown. The
transition from exposed to hidden occurs at the avoided crossing of two
saddles from different foliations, which occurs near the classical turning
point. The solid line turns to dashed or dotted where it transitions to
hidden. The dashes are for saddles that crossed a Stokes line and must be
excluded.

For the real and imaginary parts of positions on the initial manifold, the
real part of the final position is contoured and shown as a density plot (the
complex momentum initial condition for each point in the figure can be deduced
from Eq. (8)). The trajectory calculations fail in the immediate neighborhood
of singular trajectories. This shows up as extremely narrow cones devoid of
data, which correspond to sets of initial conditions for trajectories that
come too close to singular trajectories for any time up to the final
propagation time of $3\tau$.

Mathematically, the narrow cones would be lines. All such lines appear to
emanate from one of two limiting points and they extend closer to those
particular points as longer propagation times are considered. However, the
precise structure of the dynamics for $\mathcal{q}$ in the neighborhood of
$(0,0.625)$ is quite complicated and difficult to discern precisely. Modulo
the difficulties in the immediate neighborhood of these points for long times,
this Lagrangian manifold for the initial coherent state defined by Eq. (II.2)
and (6) can be partitioned into an infinity of subregions by delineating the
boundary of each region with the lines of singular trajectory initial
conditions. Within any given subregion, the initial conditions can be varied
continuously without encountering a singular trajectory. In the next section,
one particular subregion plays an important role. Finally, note that for
systems which do possess branch cuts (i.e. multi-sheeted functions describing
the dynamics), it would be much more difficult to delineate subregions than
for the quartic oscillator. We leave this for future study.

###  II.8 Features and accuracy of the semiclassical approximation

Even though semiclassical methods seem to be applied most often to make
qualitative physical arguments, it is important to point out that their full
implementation in detail leads to very quantitatively accurate approximations
for systems in the appropriate limits, e.g. short wavelength or large particle
number limits. For example, coherent state propagations in Bose-Hubbard
lattices with many sites are accurately reproduced well past the Ehrenfest
time scale [15].

The accuracy of the semiclassical approximation is illustrated in Fig. 5 for
the much simpler quartic oscillator where it is possible to compare the full
propagated wave

Figure 5: Comparison of the quantum propagation of a wave packet and the GGWPD
approximation. The particular wave packet illustrated is the one defined by
Eq. (II.2) and (6) propagated for $3\tau$. The lower two panels show the
imaginary and real parts of the classical actions, respectively. The dashed
and dotted lines indicate where saddles become hidden. The (avoided) crossing
marks the transition in how quickly the real and imaginary parts of the
classical action are changing as a function of $x$, respectively. The dashed
line indicates a saddle has crossed a Stokes line, which is indicated here by
a crossing of the real parts of two classical actions, and must be excluded.
The dotted lines indicate hidden saddles, which can be included, and may be
physically relevant if the imaginary parts of their actions are not too large.

function with its GGWPD approximation. For this example, it is easier to
ensure that all the physically relevant saddles have been identified. The
distorted $5\sigma$ spiral of Fig. 3 implies $9$ foliations, each one’s
boundaries approximately given by classical turning points; see Fig. $5$ ahead
in the companion paper where the foliations are labeled. It is also
straightforward to ensure that saddles on the wrong side of Stokes lines are
excluded. In Fig. 5, the real and imaginary parts of the contributing saddles’
classical actions are shown as a function of position and labeled by their
respective associated foliation below the comparison. Note that the foliation
label remains valid beyond the turning point caustic positions where exposed
saddles become hidden. The transition point from exposed to hidden also marks
a distinct change in the behaviors of the real and imaginary parts of the
classical actions.

It turns out that all of the potentially, physically relevant exposed and
hidden saddles reside within the central vertical subregion of Fig. 4. Let’s
call it the classical zone. Not a single physically relevant saddle, exposed
or hidden, comes from any other subregion; we suspect that this is true more
generally than just for the quartic oscillator. A curious feature of the
classical zone boundary for the quartic oscillator is that each point on it
can be thought of as a limiting or accumulation point of an infinite sequence
of points from an infinity of lines passing by getting ever closer to the
boundary. At least for the quartic oscillators, the saddles tend to be well
away from these complicated boundary zones.

Just for emphasis, all exposed saddles reside within the classical zone and
that means an infinity, even those contributing vastly too little to be
physically relevant. Since they are exposed though, it is straightforward to
implement a cutoff on their physical relevance based on the weighting of the
classical reference trajectories used to locate them. For example, the set of
saddles associated with just the nine foliations arising from within the
Wigner transform’s $5\sigma$ contour are all within just the part of the
classical zone limited by $-0.625\leq Im(q_{0})\leq 0.625$
($=5/\sqrt{2b_{\alpha}}=5/8$). Any saddles from outside this restricted
subregion, contribute much less than any of the dominant saddles from within.
This has the practical consequence of eliminating almost all of the Lagrangian
manifold from the search technique to be used to solve the two-point boundary
value problem, and therefore is a great and essential simplification.

In the central

Figure 6: Expanded view of the left side of Fig. 5. The propagated wave
function decreases quite a bit to the left of where foliations
\raisebox{-0.9pt}{8}⃝ and \raisebox{-0.9pt}{9}⃝ switch from exposed to hidden.
Just to the left of this switch, the hidden saddle \raisebox{-0.9pt}{8}⃝’s
contribution is just as significant as those of the exposed saddles
\raisebox{-0.9pt}{1\\!0}⃝ and \raisebox{-0.9pt}{1\\!1}⃝, whereas
\raisebox{-0.9pt}{9}⃝ has crossed a Stokes line and ceases to contribute.
Foliations \raisebox{-0.9pt}{1\\!0}⃝ and \raisebox{-0.9pt}{1\\!1}⃝ are from
larger-than-$5\sigma$ contours, not plotted in Fig. 3 and Fig. 5. The
semiclassical inaccuracies seen where \raisebox{-0.9pt}{9}⃝ and
\raisebox{-0.9pt}{1\\!1}⃝ cross Stokes lines can be healed (uniformized)
following Berry’s prescription [34].

region of the propagated wave function of Fig. 5, the exposed saddles dominate
the contributions of the hidden saddles. Even some of the exposed saddles
included are contributing with very small magnitudes. The exposed saddles are
the ones whose imaginary parts of the action are changing relatively slowly
(real parts relatively quickly) as a function of position. Also, the
contribution of the hidden saddles are imperceptible there, i.e. the saddles
on the correct side of Stokes lines and possessing the opposite behavior in
their real and imaginary parts as the exposed saddles’ actions.

Away from the center, the number of exposed saddles diminishes as various
turning points in the real classical dynamics are surpassed. Some of the
hidden saddles begin to get comparable in magnitude to the remaining exposed
saddles and even further out dominate the contributions to the wave function.
The left side of the wave function is magnified to illustrate these points in
Fig. 6. In the upper panel, the logarithm of the propagated wave function
magnitude is plotted. This makes the alternation of oscillatory and
exponential-like decay behaviors clearly visible. On the right, the exposed
saddles of foliations \raisebox{-0.9pt}{8}⃝ and \raisebox{-0.9pt}{9}⃝ are the
dominant contributions and they create an oscillatory “plateau” in which the
overall magnitude is barely changing. Just beyond $x=-10$, they pass a caustic
where they both convert from exposed to hidden. Foliation
\raisebox{-0.9pt}{9}⃝’s saddle contribution must be thrown away (Stokes
phenomenon) and only foliation \raisebox{-0.9pt}{8}⃝’s saddle contribution
kept. This creates a “shoulder” of exponential-like decay because its
imaginary action is changing rapidly, but not its real part. Its contribution
dominates up to the point where its decay renders it smaller than the
magnitude of the next plateau created by the exposed saddles
\raisebox{-0.9pt}{1\\!0}⃝ and \raisebox{-0.9pt}{1\\!1}⃝. This circumstance
repeats itself further to the left where \raisebox{-0.9pt}{1\\!0}⃝ and
\raisebox{-0.9pt}{1\\!1}⃝ pass a caustic and convert to being hidden. This
exponential-like decaying behavior cannot be described properly based only on
real classical transport processes. Thus, the hidden saddles are responsible
for properly describing tunneling effects in the evolution of wave packets.

This point where two saddles nearly coalesce and convert from being exposed to
hidden can be defined more precisely by where the real parts of their actions
cross [10]. Unlike in ordinary WKB where the coalescence is perfect at
classical turning point caustics, for coherent states and wave packets
generally the real and imaginary parts of the classical actions do not cross
at the same point and in that sense these are just near coalescences that
appear as avoided crossings. Even the small inaccuracies due to
discontinuously dropping a saddle can be improved upon following a
uniformization prescription due to Berry [34]. The upshot is that GGWPD
carried out in full is highly accurate for the entire propagated wave function
including the exponentially decaying tails, and even the most troublesome
locations near caustics can be improved through uniformization techniques.

##  III Justifications for implementation strategies and hidden saddles

The focus of this section is justifications of off-center, real trajectory
methods [23, 25], implementations of GGWPD [18, 19, 20, 15], and incorporating
hidden saddles. These techniques have used reasonable, physically motivated,
but heuristic arguments and numerical evidence to support their developments.
However, certain concepts such as the already mentioned idea of a complex
phase space point being near [18] a real one are somewhat vague or misleading.
It is possible to add some precision to such ideas.

###  III.1 Nearness of complex to real trajectories

We have thought of three somewhat related possible avenues of adding precision
to the concept of nearness between a complex and real trajectory.

  * 1.

The first route is to insist that “nearby” means a Newton-Raphson search
beginning from a real reference trajectory taken from some particular
foliation generated by the dynamics converges to a unique complex saddle
trajectory that can consequently be associated with that foliation. After all,
with some small technical modifications, this is how the exposed saddles are
identified [18, 20, 21]. The complex saddle trajectories do not have to
possess tiny imaginary parts for this to happen, but they do have to shadow
their real trajectory counterparts in the sense that the stability matrix
relation, Eq. (12), holds between the saddle and reference trajectories.

  * 2.

A second route, which we find a bit more compelling, is related to how the
saddle point approximation is carried out. It involves a locally quadratic
expansion of some very complicated action function about a point for which the
linear term vanishes. However, inside a domain in which this function is
behaving quadratically enough, an expansion about any point inside that domain
is possible to describe the behavior of that function. For the action
functions of dynamical systems these domains tend to be highly asymmetric. One
direction in the domain may be very "compressed", whereas another is quite
"elongated", which, of course, complicates the indeterminate idea of “near”.
There are exposed saddles with imaginary parts of their initial conditions
approaching infinity as well. Nevertheless, the quadratic expansion is used as
an approximation locally for the purposes of integration (usually with the
limits extended to $\pm\infty$). The only distinction between the two
approaches is that, except for the expansion made exactly at the saddle point,
there exists a linear term. Incorporating the linear term into the Gaussian
integral being performed would lead to the exact same value for a perfectly
quadratic function, and gives nearly the same result if the function is
sufficiently quadratic locally.

Therefore, a second definition of the complex saddle trajectory being near a
real trajectory is that the Gaussian saddle integration gives the same result
as the Gaussian integration using the related neighboring real reference
trajectory as an expansion point and incorporating the non-vanishing linear
term. Adopting this definition, the off-center real reference trajectory
method is partly a poor man’s GGWPD in which one does not bother with
performing the Newton-Raphson search to find the true saddle point, but each
contribution from the Gaussian integrals expanded about a real trajectory
gives essentially the same result as the Gaussian integral performed at the
real trajectory’s associated saddle point. It is important to note that this
does not imply that the complex saddle is near a real point. It can have very
large imaginary components of position and momentum. If a real classical
transport pathway gives an extremely small contribution, then the complex
saddle necessarily has a significant imaginary part of the complex classical
action attenuating its contribution equivalently. It is also true that this
necessarily excludes hidden saddles from the method. Therefore, one expects
accurate approximations using the techniques of [22, 23, 25] where the quantum
dynamics is not dominated by hidden saddle processes even though the true
saddles are not being used.

  * 3.

A final third route for defining nearness begins with the observation that
even though the saddles for wave packets generally do not coalesce exactly
where a caustic is encountered, there would typically be a near coalescence,
i.e. avoided crossing. As parameters are varied, a saddle remains exposed
until a caustic of some kind is encountered and it switches to hidden status
beyond the caustic. This would be expected to make the first route fail as
well since a Newton-Raphson search beginning from a real trajectory would no
longer have a unique saddle towards which it would converge. As one might have
expected, the behavior of the real and imaginary parts of the saddle’s action
function changes character beyond a caustic. Whereas, the real part of the
action is varying rapidly along the exposed part of a foliation and the
imaginary part relatively slowly, the opposite occurs for the hidden saddle,
the imaginary part varies relatively rapidly, but not the real part. The
saddle’s contribution switches from a rapidly varying phase to a rapidly
diminishing magnitude as a function of the parameter. In Figs. 5,6 it is clear
where the transition between exposed and hidden occurs as the position is
varied.

###  III.2 Incorporating hidden saddles

The semiclassical methods of [18, 20, 15, 21] are designed to locate the
physically relevant subset of exposed saddles. The question which remains is
how to identify physically relevant hidden saddles in cases where they must be
found. They are particularly relevant where exponentially decaying behavior is
dominant, such as in tunneling problems, but even in the simple quartic
oscillator example, there are locations where hidden saddles give
contributions as large as those coming from the most relevant exposed saddles,
as mentioned previously. In such cases, they also need to be incorporated into
the methods.

The technique proposed here relies on the general feature that the saddles do
not coalesce exactly in the neighborhood of caustics. The idea is to begin
with a complete set of exposed saddles, and then follow the complex saddles
continuously as a function of the natural parameters. Since the saddles do not
collide precisely, they can be followed through and beyond the caustic regions
without ambiguity. As the parameters are changed slightly, the previous saddle
is to be used as the input for a Newton-Raphson search instead of some real
trajectory. If as a parameter shifts, a particular saddle crosses the boundary
between being exposed and hidden, the Newton-Raphson search using real
trajectory input would fail whereas a search using the neighboring complex
saddle as input has no relevance to the concept of nearness to real dynamics.
Thus, there is no convergence problem as one steps through the parameter
variation. This is illustrated in Fig. 4 where the initial conditions for
saddles associated with the foliations over a range of final position are
shown. For each near approach of two exposed saddles, which is near a turning
point caustic in this case, the exposed saddles become hidden and one of them
crosses a Stokes line. The curvature of the avoided crossing introduces no
problems to a Newton-Raphson search using the neighboring saddle as a starting
point.

For the quartic oscillator this method identifies the most physically relevant
hidden saddles, they all remain within the classical zone, and they even
maintain a labelling with respect to the real classical foliations. The
homogeneity of the quartic oscillator Hamiltonian makes it a particularly
simple case. The problem of bald spots requiring complex time paths [28] does
not arise for the quartic oscillator either. Parenthetically, continuous
parameter variation would also work for locating hidden saddles behind bald
spots, but it adds an additional continuous variation variable and one has to
be aware of when a branch cut is being crossed. The idea of following exposed
saddles through caustic regions is more general than just the turning point
caustic example encountered here and extends to systems with more degrees of
freedom, whatever their dynamics, integrable, chaotic or mixed [20, 15, 21].
Nevertheless, even without considering the bald spot problem, the hidden
saddles are certainly a more complicated story in general as even the double
well creates a number of foreseeable complications. For foliations whose
trajectories pass near a caustic near the barrier multiple times in its
history, it will be necessary to consider the crossing through the caustic (or
barrier) at each instance, leading to a much larger multiplicity of hidden
saddles.

##  IV Conclusions

The semiclassical approximation can be extremely important in many physical
contexts. Of particular interest in this paper are problems in which part of
the analysis concerns the propagation of Gaussian wave packets or coherent
states in bosonic many-body systems. For nonlinear dynamical systems, the
saddle point approximation leads to a two-point boundary value problem with an
infinity of solutions, almost all of which are physically irrelevant. The
practical problems with implementing the theory only mount further as a
system’s number of degrees of freedom increases or the dynamics become
chaotic. Two techniques, i.e. off-center real trajectory methods [23, 25] and
using those resultant trajectories to identify physically relevant saddles
[18, 20, 15, 21], are placed on a firmer foundation by showing how they are
related to dynamical structure that partitions the Lagrangian manifolds
underlying wave packet and coherent state propagation.

For a simple one degree of freedom dynamical system, the pure quartic
oscillator, the wave packet / coherent state related Lagrangian manifolds can
be partitioned into an infinite number of subregions by the lines of initial
conditions of singular trajectories. One region in particular, here dubbed the
classical zone, contains all the exposed and hidden saddles discussed here.
The exposed saddles can be placed in a one-to-one correspondence with all the
real classical transport pathways, and are straightforwardly located with a
Newton-Raphson search technique fed by a single reference trajectory from each
transport pathway. A simple criterion based on the real reference trajectory’s
initial conditions and final phase point can be used to determine which
saddles contribute sufficiently for them to be considered physically relevant.
It only involves the differences between the initial phase space centroid and
the initial conditions, and the final phase point and final centroid.

The hidden saddles, where relevant, cannot be identified this way. However,
one can begin with the maximal set of exposed saddles (for the quartic
oscillator this means setting the final position to zero) and vary the
relevant parameters smoothly, here position. Following each exposed saddle
continuously with final position leads eventually to a regime where that
saddle cannot be found with a Newton-Raphson search using a real trajectory as
initial input. Nevertheless, one can follow each saddle as it moves from an
exposed region into a hidden one. Some of the hidden saddles found this way do
cross Stokes lines and need to be thrown away, but that can be evaluated by
calculating their complex action functions. For this simple dynamical example,
all the physically relevant hidden saddles necessary to calculate the
propagating wave packet / coherent state out in the exponentially decreasing
tail regions are from the classical zone and could be found through real
parameter variation of final position. This method would work in an identical
way if the quantity of interest were the overlap of the propagated wave
function with a final coherent state/ wave packet. The parameter(s) varied
would be the centroid of the final wave packet.

The dynamics of the pure quartic oscillator are extremely simple, not only
because it is a one-degree-of-freedom system, but also because the potential
is homogeneous. This gives rise to simple scaling relations between
trajectories on any pair of energy surfaces. More general dynamical systems
might be expected to give rise to further challenges not encountered in this
simple example. Indeed, one issue that has been identified is the bald spot
problem [10, 28]. This arises when a system possesses movable branch points as
a function of initial conditions, which in turn may block the Newton-Raphson
search scheme for saddles beyond some point as dynamical quantities are varied
(such as the position variable in the propagated wave function). In these
circumstances, complex time contours must be introduced in order to circumvent
the branch points and locate the saddles “behind” the bald spots [28]. Thus,
depending on the system, complex time may be unavoidable.

Further complications may arise with the introduction of multiple degrees of
freedom, where the dynamics may include KAM tori, chaotic regions, Arnol’d
diffusion, etc… The off-center trajectory and exposed saddle methods have been
shown to work in higher dimensional systems. Nevertheless, it is not known
whether the partitioning of the Lagrangian manifold found here has a
straightforward multi-degree of freedom generalization nor is it known whether
the appearance of chaotic dynamics fundamentally alters the picture of a
partitioned Lagrangian manifold. One natural extension of the current work
would be to investigate the possibility of such higher dimensional
partitionings.

## References

  * Glauber [1963] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev. 131, 2766 (1963).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, Cambridge, UK, 1997).
  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. W. Hänsch, and I. Bloch, Collapse and revival of the matter wave field of a bose-einstein condensate, Nature 419, 51 (2002).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Heller [1981] E. J. Heller, The semiclassical way to molecular spectroscopy, Acc. Chem. Res. 14, 368 (1981).
  * Gruebele and Zewail [1992] M. Gruebele and A. H. Zewail, Femtosecond wave packet spectroscopy: Coherences, the potential, and structural determination, J. Chem. Phys. 98, 883 (1992).
  * Zewail [2000] A. H. Zewail, Femtochemistry: Atomic-scale dynamics of the chemical bond, J. Phys. Chem. A 104, 5660 (2000).
  * Agostini and DiMauro [2004] P. Agostini and L. F. DiMauro, The physics of attosecond light pulses, Rep. Prog. Phys. 67, 813 (2004).
  * Klauder and Skagerstam [1985] J. R. Klauder and B.-S. Skagerstam,  _Coherent States: Applications in Physics and Mathematical Physics_ (World Scientific, Singapore, 1985).
  * Huber  _et al._ [1988] D. Huber, E. J. Heller, and R. G. Littlejohn, Generalized gaussian wave packet dynamics, schrödinger equation, and stationary phase approximation, J. Chem. Phys. 89, 2003 (1988).
  * Baranger  _et al._ [2001] M. Baranger, M. A. M. de Aguiar, F. Keck, H. J. Korsch, and B. Schellhaass, Semiclassical approximations in phase space with coherent states, J. Phys. A: Math. Gen. 34, 7227 (2001).
  * Berry and Balazs [1979] M. V. Berry and N. L. Balazs, Evolution of semiclassical quantum states in phase space, J. Phys. A 12, 625 (1979).
  * Heller [1975] E. J. Heller, Time-dependent approach to semiclassical dynamics, J. Chem. Phys. 62, 1544 (1975).
  * Heller and Tomsovic [1993] E. J. Heller and S. Tomsovic, Post-modern quantum mechanics, Physics Today 46, 38 (1993), reprinted in "Parity", Vol. 12 (8), (1993); Figures 1 & 6 reproduced in "Semiclassical Physics", Brack and Bhaduri, Addison-Wesley, (1997, USA).
  * Tomsovic  _et al._ [2018] S. Tomsovic, P. Schlagheck, D. Ullmo, J.-D. Urbina, and K. Richter, Post-ehrenfest many-body quantum interferences in ultracold atoms far-out-of-equilibrium, Phys. Rev. A 97, 061606(R) (2018), arXiv:1711.04693v2 [quant-ph].
  * Ehrenfest [1927] P. Ehrenfest, Bemerkung über die angenäherte gültigkeit der klassischen mechanik innerhalb der quantenmechanik, Zeit. Phys. 45, 455 (1927).
  * Wang and Tomsovic [2021] H. Wang and S. Tomsovic, A corrected maslov index for complex saddle trajectories, arXiv:2108.04301 [quant-ph] (2021).
  * Van Voorhis and Heller [2002] T. Van Voorhis and E. J. Heller, Nearly real trajectories in complex semiclassical dynamics, Phys. Rev. A 66, 050501(R) (2002).
  * Van Voorhis and Heller [2003] T. Van Voorhis and E. J. Heller, Similarity transformed semiclassical dynamics, J. Chem. Phys. 119, 12153 (2003).
  * Pal  _et al._ [2016] H. Pal, M. Vyas, and S. Tomsovic, Generalized gaussian wave packet dynamics: Integrable and chaotic systems, Phys. Rev. E 93, 012213 (2016), arXiv:1510.08051 [quant-ph].
  * Tomsovic [2018] S. Tomsovic, Complex saddle trajectories for multidimensional quantum wave packet/coherent state propagation: application to a many-body system, Phys. Rev. E 98, 023301 (2018), arXiv:1804.10511 [cond-mat.stat-mech].
  * O’Connor  _et al._ [1992] P. W. O’Connor, S. Tomsovic, and E. J. Heller, Semiclassical dynamics in the strongly chaotic regime - breaking the log-time barrier, Physica D 55, 340 (1992).
  * Tomsovic and Heller [1993a] S. Tomsovic and E. J. Heller, The long-time semiclassical dynamics of chaos: the stadium billiard, Phys. Rev. E 47, 282 (1993a).
  * Tomsovic and Heller [1993b] S. Tomsovic and E. J. Heller, The semiclassical construction of chaotic eigenstates, Phys. Rev. Lett. 70, 1405 (1993b).
  * Barnes  _et al._ [1994] I. M. S. Barnes, M. Nauenberg, M. Nockleby, and S. Tomsovic, Classical orbits and semiclassical wave packet propagation in the coulomb potential, J. Phys. A: Math. Gen. 27, 3299 (1994).
  * Maslov and Fedoriuk [1981] V. P. Maslov and M. V. Fedoriuk,  _Semiclassical approximation in quantum mechanics_ (Reidel Publishing Company, Dordrecht, 1981).
  * Creagh [1998] S. C. Creagh, Tunneling in two dimensions, in  _Tunneling in complex systems, Proceedings from the Institute for Nuclear Theory: Volume 5_ , edited by S. Tomsovic (World Scientific, Singapore, 1998) pp. 35–100.
  * Petersen and Kay [2014] J. Petersen and K. G. Kay, Complex time paths for semiclassical wave packet propagation with complex trajectories, J. Chem. Phys. 141, 054114 (2014).
  * Tomsovic and Heller [1991] S. Tomsovic and E. J. Heller, Semiclassical dynamics of chaotic motion: Unexpected long time accuracy, Phys. Rev. Lett. 67, 664 (1991).
  * Barnes  _et al._ [1993] I. M. S. Barnes, M. Nauenberg, M. Nockleby, and S. Tomsovic, Semiclassical theory of quantum propagation: The coulomb potential, Phys. Rev. Lett. 71, 1961 (1993).
  * Sepúlveda  _et al._ [1992] M.-A. Sepúlveda, S. Tomsovic, and E. J. Heller, Semiclassical propagation: how long can it last, Phys. Rev. Lett. 69, 402 (1992).
  * Ince [1956] E. L. Ince,  _Ordinary Differential Equations_ (Dover, Mineola, New York, 1956).
  * Polyanin and Zaitsev [2003] A. D. Polyanin and V. F. Zaitsev,  _Handbook of exact solutions for ordinary differential equations_ (Chapman and Hall/CRC, Boca Raton, 2003).
  * Berry [1989] M. V. Berry, Uniform asymptotic smoothing of stokes’ discontinuities, Proc. R. Soc. A 422, 7 (1989).