# Uniform Asymptotic Approximation Method with Pöschl-Teller Potential

Rui Pana<EMAIL_ADDRESS>John Joseph Marchettab
<EMAIL_ADDRESS>Jamal Saeeda<EMAIL_ADDRESS>Gerald
Cleaverb<EMAIL_ADDRESS>Bao-Fei Lic,d<EMAIL_ADDRESS>Anzhong Wanga 111The corresponding author<EMAIL_ADDRESS>the
corresponding author Tao Zhuc,d<EMAIL_ADDRESS>a GCAP-CASPER, Physics
Department, Baylor University, Waco, TX 76798-7316, USA
b EUCOS-CASPER, Physics Department, Baylor University, Waco, TX 76798-7316,
USA
c Institute for Advanced Physics $\&$ Mathematics, Zhejiang University of
Technology, Hangzhou, 310032, China
d United Center for Gravitational Wave Physics (UCGWP), Zhejiang University of
Technology, Hangzhou, 310032, China

###### Abstract

In this paper, we study analytical approximate solutions of the second-order
homogeneous differential equations with the existence of only two turning
points (but without poles), by using the uniform asymptotic approximation
(UAA) method. To be more concrete, we consider the Pöschl-Teller (PT)
potential, for which analytical solutions are known. Depending on the values
of the parameters involved in the PT potential, we find that the upper bounds
of the errors of the approximate solutions in general are $\lesssim 0.15\%\sim
10\%$, to the first-order approximation of the UAA method. The approximations
can be easily extended to high-order, with which the errors are expected to be
much smaller. Such obtained analytical solutions can be used to study
cosmological perturbations in the framework of quantum cosmology, as well as
quasi-normal modes of black holes.

###### pacs:

98.80.Cq, 98.80.Qc, 04.50.Kd, 04.60.Bc

##  I Introduction

A century after the first claim by Einstein that general relativity (GR) needs
to be quantized, the unification of Quantum Mechanics and GR still remains an
open question, despite enormous efforts QGs . Such a theory is necessary not
only for conceptual reasons but also for the understanding of fundamental
issues, such as the big bang and black hole singularities. Various theories
have been proposed and among them, string/M-Theory and Loop Quantum Gravity
(LQG) have been extensively investigated string ; LQG . Differences between
the two approaches are described in TTLQGS ; RHLQGS .

LQG was initially based on a canonical approach to quantum gravity (QG)
introduced earlier by Dirac, Bergmann, Wheeler, and DeWitt WDW . However,
instead of using metrics as the quantized objects WDW , LQG is formulated in
terms of densitized triads and connections, and is a non-perturbative and
background-independent quantization of GR Ashtekar86 . The gravitational
sector is described by the SU(2)-valued Ashtekar connection and its associated
conjugate momentum, the densitized triad, from which one defines the holonomy
of Ashtekar’s connection and the flux of the densitized triad. Then, one can
construct the full kinematical Hilbert space in a rigorous and well-defined
way LQG . An open question of LQG is its semiclassical limit, that is, are
there solutions of LQG that closely approximate those of GR in the
semiclassical limit?

Although the above question still remains open, concrete examples can be found
in the context of loop quantum cosmology (LQC) (For recent reviews of LQC, see
LQC_rew1 ; LQC_rew2 ; LQC_rew3 ; LQC_rew4 ; LQC_rew5 ; LQC_rew6 ; LQC_rew7 ;
LQC_rew8 ; LQC_rew9 ; LQC_rew10  and references therein). Physical
implications of LQC have also been studied using the effective descriptions of
the quantum spacetimes derived from coherent states Taveras08 , whose validity
has been verified numerically for various spacetimes AG15 ; numlsu-1 ,
especially for states sharply peaked on classical trajectories at late times
KKL20 . The effective dynamics provide a definitive answer on the resolution
of the big bang singularity generic ; SVV06 ; ZL07 ; AS10 ; AS11 ; CZ15 ,
replaced by a quantum bounce when the energy density of matter reaches a
maximum value determined purely by the underlying quantum geometry.

To connect LQC with observations, cosmological perturbations in LQC have been
also investigated intensively in the past decade, and a variety of different
approaches to extend LQC to include cosmological perturbations have been
developed. These include the dressed metric Agullo:2012sh ; Agullo:2012fc ;
Agullo:2013ai , hybrid Fernandez-Mendez:2012poe ; Fernandez-Mendez:2014raa ;
Martinez:2016hmn ; ElizagaNavascues:2020uyf , deformed algebra Bojowald:2007hv
; Bojowald:2007cd ; Bojowald:2008gz ; Bojowald:2008jv  and separate universe
Wilson-Ewing:2012dcf ; Wilson-Ewing:2015sfx  approaches. For a brief review on
each of these approaches, we refer readers to LQC_rew9 .

One of the major challenges in the studies of cosmological perturbations in
LQC is how to solve for the mode functions $\mu_{k}$ from the modified
Mukhanov-Sasaki equation. So far, it has mainly been done numerically LQC_rew1
; LQC_rew2 ; LQC_rew3 ; LQC_rew4 ; LQC_rew5 ; LQC_rew6 ; LQC_rew7 ; LQC_rew8 ;
LQC_rew9 ; LQC_rew10 . However, this is often required to be conducted with
high-performance computational resources AM15 , which are not accessible to
general audience.

In the past decade, we have systematically developed the uniform asymptotic
approximation (UAA) method initially proposed by Olver Olver97 ; Olver56 ;
Olver75 , and applied it successfully to various circumstances P0 ; P1 ; P2 ;
P3a ; P3 ; P4 ; P5 ; P6 ; P7 ; P8 ; P8a ; P10 ; P11 ; P12 ; P13 ; P14 ; P15 ;
Zhu:2022dfq ; Li:2022xww  222It should be noted that the first application of
the UAA method to cosmological perturbations in GR was carried by Habib et al,
Habib:2002yi ; Habib:2004kc .. In this paper, we shall continuously work on it
by considering the case in which the effective potential has only zero points
but without singularities. To be more concrete, we shall consider the Pöschl-
Teller (PT) potential, for which analytical solutions are known dong_wave_2011
. The consideration of this potential is also motivated from the studies of
cosmological perturbations in dressed metric and hybrid approaches ZWCKS17 ;
WZW18 , in which it was shown explicitly that the potentials for the mode
functions can be well approximated by the PT potential with different choices
of the PT parameters.

The rest of the paper is organized as follows: In Sec. II we provide a brief
review of the UAA method with two turning points, and show that the first-
order approximate solution will be described by the parabolic cylindrical
functions. In Sec. III we construct the explicit approximate analytical
solutions with the PT potential, and find that the parameter space can be
divided into three different cases: A) $k^{2}\gg\beta^{2}$, B)
$k^{2}\simeq\beta^{2}$, and C) $k^{2}\ll\beta^{2}$, where $k$ and $\beta$ are
real constants. After working out the error control function $\mathscr{T}$
[cf. Appendix C] in each case, we are able to determine the parameter $q_{0}$,
introduced in the process of the UAA method in order to minimize the errors.
Then, we show the upper bounds of errors of our approximate solutions with
respect to the exact one, given in Appendix B. In particular, in Case A), the
upper bounds are $\lesssim 0.15\%$, while in Case B) they are no larger than
$10\%$. In Case C), the errors are also very small, except the minimal points
[cf. Fig. 10], at which the approximate solutions deviate significantly from
the analytical one. The causes of such large errors are not known, and still
under our investigations. In each of these three cases, we also develop our
numerical codes, and find that the numerical solutions trace the exact one
very well, and the upper bounds of errors are always less than $10^{-4}\%$ in
each of the three cases. The paper is ended in Sec. IV, in which our main
conclusions are summarized. There are also three appendices, A, B, and C, in
which some mathematical formulas are presented.

##  II The Uniform Asymptotic Approximation Method

Let us start with the following second-order differential equation

$\displaystyle\frac{d^{2}\mu_{k}(y)}{dy^{2}}=f(y)\mu_{k}(y).$ (2.1)

It should be noted that all second-order linear homogeneous ordinary
differential equations (ODEs) can be written in the above form by properly
choosing the variable $y$ and $\mu_{k}(y)$. Instead of working with the above
form, we introduce two functions $g(y)$ and $q(y)$, so that the function
$f(y)$ takes the form Olver97

$\displaystyle f(y)=\lambda^{2}g(y)+q(y),$ (2.2)

where $\lambda$ is a large positive dimensionless constant and serves as a
bookmark, so we can expand $\mu_{k}(y)$ as

$\mu_{k}(y)=\sum_{n=0}^{\infty}{\frac{\mu^{(n)}_{k}(y)}{\lambda^{n}}}.$ (2.3)

After all the calculations are done, one can always set $\lambda=1$ by simply
absorbing the factor $\lambda^{-n}$ into $\mu^{(n)}_{k}(y)$. It should be
noted that there exist cases in which the above expansion does not converge,
and in these cases we shall expand $\mu_{k}(y)$ only to finite terms, say,
${\cal{N}}$, so that $\mu_{k}(y)$ is well approximated by the sum of these
${\cal{N}}$ terms. On the other hand, the main reason to introduce two
functions $g(y)$ and $q(y)$, instead of only $f(y)$, is to minimize errors by
properly choosing $g(y)$ and $q(y)$.

In general, the function $g(y)$ has singularities and/or zeros in the interval
of our interest. We call the zeros and singularities of $g(y)$ as turning
points and poles, respectively. The uniform asymptotic approximate (UAA)
solutions of $\mu_{k}(y)$ depend on the properties of $g(y)$ around their
poles and turning points Olver56 ; Olver97 ; Olver75 . The cases in which
$g(y)$ has both poles and turning points were studied in detail in P2 ; P3 ;
P8 , so in this paper we shall focus ourselves on the cases where
singularities are absent and only turning points exist. As to be shown below,
the treatments of these cases will be different from the ones considered in P2
; P3 ; P8 . In particular, in our previous studies the function $q(y)$ was
uniquely determined by requiring that the error control function be finite and
minimized at the poles, while in the current cases no such poles exist. So, to
fix $q(y)$, other analyses of the error control function must be carried out.

Figure 1: The function $\zeta(y)$ vs $y$, which is assumed to be always an
increasing function of $y$.

###  II.1 The UAA Method

The UAA method includes three major steps: (i) the Liouville transformations;
(ii) the minimization of the error control function; and (iii) the choice of
the function $y(\zeta)$, where $\zeta$ is a new variable. In the following, we
shall consider each of them separately.

####  II.1.1 The Liouville Transformations

The Liouville transformations consist of introducing a new variable
$\zeta(y)$, which is assumed that the inverse $y=y(\zeta)$ always exists and
is thrice-differentiable. Without loss of the generality, we also assume that
$y(\zeta)$ is a monotonically increasing function [cf. Fig. 1]. Then, in terms
of $U(\zeta)$, which is defined by

$U(\zeta)\equiv\dot{y}^{-1/2}\mu_{k},$ (2.4)

Eq.(2.1) takes the form,

$\frac{d^{2}U(\zeta)}{d\zeta^{2}}=\left[\lambda^{2}\dot{y}^{2}g+\psi(\zeta)\right]U(\zeta),$
(2.5)

where

$\displaystyle\dot{y}\equiv\frac{dy(\zeta)}{d\zeta}>0,\quad\zeta^{\prime}(y)\equiv\frac{d\zeta(y)}{dy}=\frac{1}{\dot{y}},$
(2.6)

and

$\displaystyle\psi(\zeta)$ $\displaystyle\equiv$
$\displaystyle\dot{y}^{2}q+\dot{y}^{1/2}\frac{d^{2}}{d\zeta^{2}}\left(\dot{y}^{-1/2}\right)$
(2.7) $\displaystyle=$
$\displaystyle\dot{y}^{2}q-\dot{y}^{3/2}\frac{d^{2}}{dy^{2}}\left(\dot{y}^{1/2}\right)\equiv\psi(y).$

It should be noted that Eqs.(2.1) and (2.5) are completely equivalent, and so
far no approximations are taken. However, the advantage of the form of
Eq.(2.5) is that, by properly choosing $q(y)$, the term
$\left|\psi(\zeta)\right|$ can be much smaller than
$\left|\lambda^{2}\dot{y}^{2}g\right|$, that is,

$\displaystyle\left|\frac{\psi}{\lambda^{2}\dot{y}^{2}g}\right|\ll 1,$ (2.8)

so that the exact solution of Eq.(2.1) can be well approximated by the first-
order solution of Eq.(2.5) with $\psi(\zeta)=0$. This immediately raises the
question: how to choose $q(y)$ so that the condition (2.8) holds. To explain
this in detail, let us move onto the next subsection.

####  II.1.2 Minimization of Errors

To minimize the errors, let us first introduce the error control function
Olver56 ; Olver97 ; Olver75 ; P2 ; P3 ; P8

$\mathscr{T}(\zeta)\equiv-\int\frac{\psi(\zeta)}{{\left|\dot{y}^{2}g\right|^{1/2}}}d\zeta.$
(2.9)

Then, introducing the free parameters $a_{n}$ and $b_{n}$ into the functions
$g(y)$ and $q(y)$, so we have

$g(y)=g\left(y,a_{n}\right),\quad q(y)=q\left(y,b_{n}\right),$ (2.10)

where $n=1,2,...,N$, with $N$ being an integer. It is clear that for such
chosen $g(y)$ and $q(y)$, the error control function $\mathscr{T}(\zeta)$ will
also depend on $a_{n}$ and $b_{n}$. To minimize the errors, one way is to
minimize the error control function by properly choosing $a_{n}$ and $b_{n}$,
so that

$\displaystyle\frac{\partial\mathscr{T}\left(\zeta,a_{n},b_{n}\right)}{\partial
a_{n}}=0,\;\;\;$
$\displaystyle\frac{\partial\mathscr{T}\left(\zeta,a_{n},b_{n}\right)}{\partial
b_{n}}=0,$ (2.11) $\displaystyle~{}~{}~{}~{}~{}(n=1,2,...,N).$

####  II.1.3 Choice of $y(\zeta)$

On the other hand, the errors also depend on the choice of $y(\zeta)$, which
in turn sensitively depends on the properties of the functions $g(y)$ and
$q(y)$ near their poles and turning points. In addition, it must be chosen so
that the resulting equation of the first-order approximation (obtained by
seting $\psi(\zeta)=0$) can be solved explicitly (in terms of known
functions). Considering all the above, it has been found that $y(\zeta)$ can
be chosen as Olver56 ; Olver97 ; Olver75 ; P2 ; P3 ; P8

$\displaystyle\dot{y}^{2}g=\begin{cases}\text{sgn}(g),&\text{zero turning
point},\cr\zeta,&\text{one turning
point},\cr\zeta_{0}^{2}-\zeta^{2},&\text{two turning
point},\cr\end{cases}~{}~{}~{}~{}$ (2.12)

in the cases with zero, one and two turning points, respectively. Here
$\text{sgn}(g)=1$ for $g>0$ and $\text{sgn}(g)=-1$ for $g<0$.

In the rest of this paper, we shall consider only the cases with two turning
points.

###  II.2 UAA Method for Two Turning Points

For the cases with two turning points, we can always write $g(y)$ as

$\displaystyle g(y)=p(y)(y-y_{1})(y-y_{2}),$ (2.13)

where $y_{1}$ and $y_{2}$ are the two turning points, and $p(y)$ is a function
of $y$ with $p(y_{i})\not=0,\;(i=1,2)$. In general, according to the
properties of $y_{1}$ and $y_{2}$, we can divide all the cases into three
different subclasses:

  1. 1.

$y_{1}$ and $y_{2}$ are two distinct real roots of $g(y)=0$;

  2. 2.

$y_{1}=y_{2}$, a double real root of $g(y)=0$; and

  3. 3.

$y_{1}$ and $y_{2}$ are two complex roots of $g(y)=0$. Since $g(y)$ is real,
in this case these two roots must be complex conjugate, $y_{1}=y_{2}^{*}$.

To apply the UAA method to Eq.(2.6), we assume that the following conditions
are satisfied P2 ; P3 ; P8 :

  * •

When far away from any of the two turning points, we require

$\left|\frac{q(y)}{g(y)}\right|\ll 1.$ (2.14)

  * •

When near any of these two points, we require

$\left|\frac{q(y)(y-y_{i})}{g(y)}\right|\ll 1,\;(i=1,2),$ (2.15)

provided that the two turning points are far away from each other, that is,
when $\left|y_{1}-y_{2}\right|\gg 1$.

  * •

If the two turning points are close to each other,
$\left|y_{1}-y_{2}\right|\simeq 0$, then near these points we require

$\left|\frac{q(y)(y-y_{1})(y-y_{2})}{g(y)}\right|\ll 1.$ (2.16)

Figure 2: The function ${g}(y)$ defined by Eq. (3.3) for different choices of
$k$ and $\beta$. In particular, the dotted black line denotes the case
$k^{2}<\beta^{2}$, and the solid blue line denotes the case $k^{2}=\beta^{2}$,
while the dash-dotted red line denotes the case $k^{2}>\beta^{2}$.

It should be noted that, when $\left|y_{2}-y_{1}\right|\gg 1$, the two turning
points are far away, and each of them can be treated as an isolated single
turning point Olver56 ; Olver97 . In addition, without loss of generality, we
assume that $g(y)<0$ for $y>y_{2}$ or $y<y_{1}$, when $y_{1}$ and $y_{2}$ are
real. When $y_{2}$ and $y_{1}$ are complex conjugate, we assume that $g(y)<0$
[cf. Fig. 2]. Then, in this case we adopt a method to treat all these three
classes listed above together Olver75 ; P2 ; P3 ; P8 . In particular, we
choose $\dot{y}^{2}g$ as

$\displaystyle\dot{y}^{2}g=\zeta_{0}^{2}-\zeta^{2}\begin{cases}>0,&g>0,\cr=0,&g=0,\cr<0,&g<0,\cr\end{cases}$
(2.17)

so that $\zeta$ is an increasing function of $y$ [cf. Fig. 1] and

$\sqrt{|g(y)|}\;dy=\sqrt{\left|\zeta_{0}^{2}-\zeta^{2}\right|}\;d\zeta.$
(2.18)

When we integrate the above equation, without loss of the generality, we shall
choose the integration constants so that

$\zeta(y_{1})=-\zeta_{0},\quad\zeta(y_{2})=\zeta_{0}.$ (2.19)

Then, we find that

$\displaystyle\zeta_{0}^{2}=\begin{cases}>0,&\text{$y_{1,2}$ real, and
$y_{1}\not=y_{2}$},\cr=0,&\text{$y_{1,2}$ real, and
$y_{1}=y_{2}$},\cr<0,&\text{$y_{1,2}$ complex},\cr\end{cases}$ (2.20)

with

$\displaystyle\zeta_{0}^{2}$ $\displaystyle=$
$\displaystyle\pm\frac{2}{\pi}\int_{y_{1}}^{y_{2}}\sqrt{|g(y)|}dy$ (2.21)
$\displaystyle=$
$\displaystyle\pm\frac{2}{\pi}\int_{-\zeta_{0}}^{\zeta_{0}}\sqrt{\left|\zeta_{0}^{2}-\zeta^{2}\right|}d\zeta,$

where $``+"$ corresponds to the cases that the two turning points $y_{1}$ and
$y_{2}$ are both real, and $``-"$ to the cases that the two turning points
$y_{1}$ and $y_{2}$ are complex conjugate. When $y_{1}$ and $y_{2}$ are
complex conjugate, the integration of Eq.(2.21) is along the imaginary axis
Olver75 . When the two real roots are equal, we have $\zeta_{0}=0$.

To proceed further, let us first derive the relation between $\zeta(y)$ and
$y$ by first integrating the right-hand side of Eq.(2.18). To this goal, it is
found easier to distinguish the case in which $y_{1}$ and $y_{2}$ are real
from the one in which they are complex conjugate.

####  II.2.1 When $y_{1,2}$ Are Real

Let us first consider the case when $y_{1}$ and $y_{2}$ are real. Then, when
$y>y_{2}$, we have $\zeta(y)>\zeta_{0}$ [cf. Fig.1]. Hence, from Eq. (2.18) we
find

$\displaystyle\int_{y_{2}}^{y}\sqrt{-g(y^{\prime})}dy^{\prime}=\int_{\zeta_{0}}^{\zeta}\sqrt{v^{2}-\zeta_{0}^{2}}dv$
$\displaystyle\;\;\;=\frac{1}{2}\zeta\sqrt{\zeta^{2}-\zeta_{0}^{2}}-\frac{\zeta_{0}^{2}}{2}\ln\left(\frac{\zeta+\sqrt{\zeta^{2}-\zeta_{0}^{2}}}{\zeta_{0}}\right)$
$\displaystyle\;\;\;=\frac{1}{2}\zeta\sqrt{\zeta^{2}-\zeta_{0}^{2}}-\frac{\zeta_{0}^{2}}{2}\operatorname{arcosh}{\left(\frac{\zeta}{\zeta_{0}}\right)},\;(y\geq
y_{2}).~{}~{}~{}~{}~{}~{}$ (2.22)

When $y\leq y_{1}$, we have $\zeta(y)\leq-\zeta_{0}$. Then, from Eq. (2.18) we
find

$\displaystyle\int^{y_{1}}_{y}\sqrt{-g(y^{\prime})}dy^{\prime}=\int^{-\zeta_{0}}_{\zeta}\sqrt{v^{2}-\zeta_{0}^{2}}dv$
$\displaystyle\;\;\;=-\frac{1}{2}\zeta\sqrt{\zeta^{2}-\zeta_{0}^{2}}+\frac{\zeta_{0}^{2}}{2}\ln\left(\frac{-\zeta-\sqrt{\zeta^{2}-\zeta_{0}^{2}}}{\zeta_{0}}\right)$
$\displaystyle\;\;\;=-\frac{1}{2}\zeta\sqrt{\zeta^{2}-\zeta_{0}^{2}}-\frac{\zeta_{0}^{2}}{2}\ln\left(\frac{-\zeta+\sqrt{\zeta^{2}-\zeta_{0}^{2}}}{\zeta_{0}}\right)$
$\displaystyle\;\;\;=-\frac{1}{2}\zeta\sqrt{\zeta^{2}-\zeta_{0}^{2}}-\frac{\zeta_{0}^{2}}{2}\operatorname{arcosh}{\left(-\frac{\zeta}{\zeta_{0}}\right)},(y\leq
y_{1}).~{}~{}~{}~{}~{}~{}~{}$ (2.23)

When $y_{1}\leq y\leq y_{2}$, we have $-\zeta_{0}<\zeta(y)<\zeta_{0}$, and

$\displaystyle\int_{y_{1}}^{y}\sqrt{g(y^{\prime})}dy^{\prime}=\int_{-\zeta_{0}}^{\zeta}\sqrt{\zeta_{0}^{2}-v^{2}}dv=\frac{1}{2}\zeta\sqrt{\zeta_{0}^{2}-\zeta^{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\frac{\zeta_{0}^{2}}{2}\arccos\left(-\frac{\zeta}{\zeta_{0}}\right),\;(y_{1}\leq
y\leq y_{2}).$ (2.24)

####  II.2.2 When $y_{1,2}$ Are Complex Conjugate

Now let us turn to consider the case when $y_{1}$ and $y_{2}$ are complex. For
this case $\zeta_{0}^{2}$ is always negative, $\zeta_{0}^{2}<0$, thus from Eq.
(2.4) we find Olver75

$\displaystyle\int_{0}^{y}\sqrt{-g(y^{\prime})}dy^{\prime}=\int_{0}^{\zeta}\sqrt{\zeta^{2}-\zeta_{0}^{2}}d\zeta$
$\displaystyle\;\;\;=\frac{1}{2}\zeta\sqrt{\zeta^{2}-\zeta_{0}^{2}}-\frac{\zeta_{0}^{2}}{2}\ln\left(\frac{\zeta+\sqrt{\zeta^{2}-\zeta_{0}^{2}}}{|\zeta_{0}|}\right).$
(2.25)

####  II.2.3 The First-order Approximate Solutions

With the choice of Eq.(2.17), we find that Eq. (2.6) reduces to

$\displaystyle\frac{d^{2}U}{d\zeta^{2}}=\Big{[}\lambda^{2}\left(\zeta_{0}^{2}-\zeta^{2}\right)+\psi(\zeta)\Big{]}U,$
(2.26)

where we assume that $\zeta\in(-\zeta_{2},\zeta_{2})$, with $\zeta_{2}$ being
a real and positive constant, which can be arbitrarily large
$\zeta_{2}\rightarrow\infty$.

Neglecting the $\psi(\zeta)$ term, we find that the approximate solutions can
be expressed in terms of the parabolic cylinder functions
$W(\frac{1}{2}\lambda\zeta_{0}^{2},\pm\sqrt{2\lambda}\zeta)$ Olver75 , and are
given by

$\displaystyle U(\zeta)$ $\displaystyle=$
$\displaystyle\alpha_{k}\Bigg{\\{}W\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)+\epsilon_{1}\Bigg{\\}}$
(2.27)
$\displaystyle+\beta_{k}\Bigg{\\{}W\left(\frac{1}{2}\lambda\zeta_{0}^{2},-\sqrt{2\lambda}\zeta\right)+\epsilon_{2}\Bigg{\\}},$

from which we have

$\displaystyle\mu_{k}(y)$ $\displaystyle=$
$\displaystyle\alpha_{k}\left(\frac{\zeta^{2}-\zeta_{0}^{2}}{-g(y)}\right)^{\frac{1}{4}}\left[W\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)+\epsilon_{1}\right]$
$\displaystyle+\beta_{k}\left(\frac{\zeta^{2}-\zeta_{0}^{2}}{-g(y)}\right)^{\frac{1}{4}}\left[W\left(\frac{1}{2}\lambda\zeta_{0}^{2},-\sqrt{2\lambda}\zeta\right)+\epsilon_{2}\right],$

where $\alpha_{k}$ and $\beta_{k}$ are two integration constants,
$\epsilon_{1}$ and $\epsilon_{2}$ are the errors of the corresponding
approximate solutions, whose upper bounds are given by Eqs.(Appendix A. Upper
Bounds of Errors) and (A.2) in Appendix A.

For the choice of Eq.(2.17), we find that the associated error control
function defined by Eq.(2.9) now takes the form

$\displaystyle\mathscr{T}(\zeta)$ $\displaystyle=$
$\displaystyle-\int^{\zeta}\left\\{\frac{q}{g}-\frac{5}{16}\frac{g^{\prime
2}}{g^{3}}+\frac{1}{4}\frac{g^{\prime\prime}}{g^{2}}\right\\}\sqrt{v^{2}-\zeta_{0}^{2}}dv$
$\displaystyle+\int^{\zeta}\left\\{\frac{5\zeta_{0}^{2}}{4(v^{2}-\zeta_{0}^{2})^{3}}+\frac{3}{4(v^{2}-\zeta_{0}^{2})^{2}}\right\\}\sqrt{v^{2}-\zeta_{0}^{2}}dv$
$\displaystyle=$
$\displaystyle-\int^{y}\left\\{\frac{q}{g}-\frac{5}{16}\frac{g^{\prime
2}}{g^{3}}+\frac{1}{4}\frac{g^{\prime\prime}}{g^{2}}\right\\}\sqrt{-g}dy^{\prime}$
$\displaystyle+\int^{\zeta}\left\\{\frac{5\zeta_{0}^{2}}{4(v^{2}-\zeta_{0}^{2})^{5/2}}+\frac{3}{4(v^{2}-\zeta_{0}^{2})^{3/2}}\right\\}dv,$

for $g<0$, and

$\displaystyle\mathscr{T}(\zeta)$ $\displaystyle=$
$\displaystyle\int^{\zeta}\left\\{\frac{q}{g}-\frac{5}{16}\frac{g^{\prime
2}}{g^{3}}+\frac{1}{4}\frac{g^{\prime\prime}}{g^{2}}\right\\}\sqrt{\zeta_{0}^{2}-v^{2}}dv$
$\displaystyle-\int^{\zeta}\left\\{\frac{5\zeta_{0}^{2}}{4(v^{2}-\zeta_{0}^{2})^{3}}+\frac{3}{4(v^{2}-\zeta_{0}^{2})^{2}}\right\\}\sqrt{\zeta_{0}^{2}-v^{2}}dv$
$\displaystyle=$
$\displaystyle\int^{y}\left\\{\frac{q}{g}-\frac{5}{16}\frac{g^{\prime
2}}{g^{3}}+\frac{1}{4}\frac{g^{\prime\prime}}{g^{2}}\right\\}\sqrt{g}dy^{\prime}$
$\displaystyle+\int^{\zeta}\left\\{\frac{5\zeta_{0}^{2}}{4(\zeta_{0}^{2}-v^{2})^{5/2}}-\frac{3}{4(\zeta_{0}^{2}-v^{2})^{3/2}}\right\\}dv,$

for $g>0$.

##  III UAA Solutions with the Pöschl-Teller Potential

To study the case in which only turning points exist, in this paper we
consider the second-order differential equation (2.1) with a Pöschl-Teller
(PT) potential ZWCKS17 ; WZW18

$\left(\lambda^{2}g+q\right)=-\left(k^{2}-\frac{\beta_{0}^{2}}{\cosh^{2}(\alpha
y)}\right),$ (3.1)

as in this case exact solutions exist, where $k$ is the comoving wavenumber,
and $\beta_{0}$ is a real and positive constant. The two parameters
$\beta_{0}$ and $\alpha$ determine the height and the spread of the PT
potential, respectively. Under the rescaling $\alpha y\rightarrow y$, the
$\alpha$ parameter can be absorbed into the wavenumber $k$ and $\beta_{0}$ by
redefining $\left(k/\alpha\rightarrow
k,\beta_{0}/\alpha\rightarrow\beta_{0}\right)$. As a result, there is no loss
of generality to set $\alpha=1$ from now on. Then, the exact solutions in this
case exist, and are presented in Appendix B.

On the other hand, to apply the UAA method to this case, and to minimize the
errors of the analytic approximate solutions, we tentatively choose $q$ as

$q=\frac{q^{2}_{0}}{\cosh^{2}(y)},$ (3.2)

where $q_{0}$ is a free parameter, to be determined below by minimizing the
error control function (2.9) with the choice of $\dot{y}^{2}g$ given by
Eq.(2.17). Then, we have

$g(y)=\frac{\beta^{2}}{\cosh^{2}(y)}-k^{2},$ (3.3)

where $\beta\equiv\sqrt{\beta^{2}_{0}-q^{2}_{0}}$. In this paper, without loss
of generality, we shall choose $q_{0}$ so that $\beta$ is always real, that is

$\beta^{2}\equiv\beta_{0}^{2}-q_{0}^{2}>0.$ (3.4)

Thus, from $g(y)=0$ we find that the two roots are given by

$y_{i}=\pm\cosh^{-1}\frac{\beta}{k}=\pm\cosh^{-1}\frac{\sqrt{\beta^{2}_{0}-q^{2}_{0}}}{k}.$
(3.5)

It is clear that, depending on the relative magnitudes of $\beta_{0}$ and $k$,
as well as the choices of $q_{0}$, two turning points can be either complex or
real. In Fig. 2, we plot out the three different cases, $k^{2}<\beta^{2}$,
$k^{2}=\beta^{2}$, and $k^{2}>\beta^{2}$, from which it can be seen clearly
that the two turning points are real and different for $k^{2}<\beta^{2}$, real
and equal for $k^{2}=\beta^{2}$, and complex conjugate for $k^{2}>\beta^{2}$,
respectively. Then, from Eqs.(3.2) and (3.3), we find that

$\left|\frac{q(y)}{g(y)}\right|=\left|\frac{q_{0}^{2}}{\beta^{2}-k^{2}\cosh^{2}(y)}\right|\simeq
q_{0}^{2}e^{-2|y|},$ (3.6)

for $|y|\gg 0$, and

$\left|\frac{q(y)(y-y_{i})}{g(y)}\right|\simeq\frac{q_{0}^{2}}{y+y_{j}},(i\not=j),$
(3.7)

for $|y|\simeq|y_{i}|$ and $|y_{1}-y_{2}|\gg 1$, and

$\left|\frac{q(y)(y-y_{1})(y-y_{2})}{g(y)}\right|\simeq q_{0}^{2},$ (3.8)

for $|y|\simeq|y_{1}|$ and $|y_{1}-y_{2}|\simeq 0$. In the following, let us
consider the three cases: (a) $k^{2}\gg\beta^{2}$; (b) $k^{2}\simeq\beta^{2}$;
and (c) $\beta^{2}\gg k^{2}$, separately.

###  III.1 $k^{2}\gg\beta^{2}$

In this case, we have $g(y)$ is always negative, $g(y)<0$, so that the two
turning points of $g(y)=0$ are complex conjugate and are given by

$y_{1}=y_{2}^{*}=-i\cos^{-1}\left(\frac{\beta}{k}\right)\simeq-\frac{i\pi}{2}.$
(3.9)

As discussed in the last section, now $\zeta_{0}^{2}<0$, for which Eq.(2.26)
can be cast in the form

$\frac{d^{2}W(\zeta)}{d^{2}\zeta}=\Big{\\{}-\lambda^{2}\left(\zeta^{2}+\hat{\zeta}^{2}_{0}\right)+\psi\Big{\\}}W(\zeta),$
(3.10)

where $\hat{\zeta}^{2}_{0}\equiv-\zeta^{2}_{0}>0$. Note that in writing down
the above equation, we had replaced $U$ by $W$. In addition, the new variable
$\zeta$ is related to $y$ via

$\displaystyle\int^{y}_{0}\sqrt{-g(y)}dy=\int^{\zeta}_{0}\sqrt{v^{2}+\hat{\zeta}^{2}_{0}}dv=\frac{1}{2}\hat{\zeta}^{2}_{0}\ln\left(\zeta+\sqrt{\zeta^{2}+\hat{\zeta}^{2}_{0}}\right)+\frac{1}{2}\zeta\sqrt{\zeta^{2}+\hat{\zeta}^{2}_{0}}-\frac{1}{2}\hat{\zeta}^{2}_{0}\ln\hat{\zeta}_{0},$
(3.11)

from which we find that $\hat{\zeta}_{0}$ is given explicitly by

$\hat{\zeta}^{2}_{0}=2\left(k-\beta\right)>0.$ (3.12)

Moreover, in the case of the PT potential, the integration of Eq. (3.11) can
be carried out explicitly, giving

$\displaystyle\int^{y}_{0}dy\sqrt{-g}=\epsilon_{y}\sqrt{1-x^{2}}\sqrt{k^{2}-\beta^{2}}\times\mathrm{AppellF_{1}}\left(\frac{1}{2},-\frac{1}{2},1,\frac{3}{2};\frac{1-x^{2}}{1-k^{2}/\beta^{2}},1-x^{2}\right),$
(3.13)

where $\epsilon_{y}$ denotes the sign of $y$ with $x\equiv 1/\cosh(y)$, and
$\mathrm{AppellF_{1}}$ is the Appell hypergeometric function. Ignoring the
$\psi$ term in Eq. (3.10), we find the general solution

$\displaystyle\mu_{k}(y)=\left(\frac{\zeta^{2}+\hat{\zeta}^{2}_{0}}{-g(y)}\right)^{1/4}\Bigg{\\{}\alpha_{k}W\left(-\frac{\hat{\zeta}^{2}_{0}}{2},\sqrt{2}\zeta\right)+\beta_{k}W\left(-\frac{\hat{\zeta}^{2}_{0}}{2},-\sqrt{2}\zeta\right)\Bigg{\\}},$
(3.14)

where $W$ denotes the Weber parabolic cylinder function WeberF , and
$\alpha_{k}$ and $\beta_{k}$ are two integration parameters which generally
depend on the comoving wavenumber $k$.

The validity of the analytic solution (3.14) depends on the criteria given by
Eqs.(2.14) - (2.16), while its accuracy can be predicted by the error control
function $\mathscr{T}$. In the current case, we find that $\mathscr{T}$ of
Eq.(II.2.3) can be written as a combination of three terms as that given by
Eqs.(C.3) 333In this case, the associated error control function is
$\mathscr{V}_{\zeta_{1},\zeta}(\mathscr{T})$ for any given $\zeta_{1}$, where
$\zeta_{1}\in\left(-\infty,\infty\right)$ Olver75 . In this paper, we choose
$\zeta_{1}=0$, so the integrations will be carried out in the interval
$\zeta\in[0,\infty)$, corresponding to $y\in[0,\infty)$. Due to the symmetry
of the equation, one can easily obtain the solutions for the region
$y\in(-\infty,0]$ by simply replacing $y$ by $-y$ (or $\zeta$ by $-\zeta$).,
where

$\displaystyle\mathscr{T}_{1}$ $\displaystyle=$
$\displaystyle\int^{y}_{0}\frac{q}{\sqrt{-g}}dy=\frac{q^{2}_{0}\epsilon_{y}}{\beta}\ln\left(\frac{\sqrt{1-x^{2}}\beta+\sqrt{k^{2}-\beta^{2}x^{2}}}{\sqrt{k^{2}-\beta^{2}}}\right),$
$\displaystyle\mathscr{T}_{2}$ $\displaystyle=$
$\displaystyle\int^{y}_{0}\left(\frac{5g^{\prime
2}}{16g^{3}}-\frac{g^{{}^{\prime\prime}}}{4g^{2}}\right)\sqrt{-g}dy=-\epsilon_{y}\Bigg{\\{}\frac{1}{4\beta}\ln\left(\frac{\sqrt{1-x^{2}}\beta+\sqrt{k^{2}-x^{2}\beta^{2}}}{\sqrt{k^{2}-\beta^{2}}}\right)-\frac{\sqrt{1-x^{2}}A}{12(k^{2}-\beta^{2})(k^{2}-\beta^{2}x^{2})^{3/2}}\Bigg{\\}},$
$\displaystyle\mathscr{T}_{3}$ $\displaystyle=$
$\displaystyle\int^{\zeta}_{0}\left(\frac{-5\hat{\zeta}^{2}_{0}}{4\left(v^{2}+\hat{\zeta}^{2}_{0}\right)^{5/2}}+\frac{3}{4\left(v^{2}+\hat{\zeta}^{2}_{0}\right)^{3/2}}\right)dv=-\frac{\zeta\left(\zeta^{2}+6\hat{\zeta}^{2}_{0}\right)}{12\hat{\zeta}^{2}_{0}\left(\zeta^{2}+\hat{\zeta}^{2}_{0}\right)^{3/2}},$
(3.15)

where $A$ is given by Eq.(C.5). It should be noted that $\mathscr{T}_{1}$,
$\mathscr{T}_{2}$ and $\mathscr{T}_{3}$ given in Eq. (III.1) all vanish when
$y=0$ (for which we have $x=1$ and $\zeta=0$), that is,

$\mathscr{T}(\zeta=0)=0.$ (3.16)

Besides, as the PT potential is an even function, the error control function
is antisymmetric about the origin, namely, $\mathscr{T}(-y)=-\mathscr{T}(y)$.
As a result, we will study its behavior only on the positive $y$ axis, $y\geq
0$. With the help of Eq. (3.11), the numeric value of the error control
function at any point $y>0$ can be found from Eq.(III.1). In particular, for
$\beta/k\ll 1$, we find that

$\displaystyle\mathscr{T}=\frac{q^{2}_{0}}{k}\sqrt{1-x^{2}}-\frac{\zeta\left(\zeta^{2}+6\hat{\zeta}^{2}_{0}\right)}{12\hat{\zeta}^{2}_{0}\left(\zeta^{2}+\hat{\zeta}^{2}_{0}\right)^{3/2}}+{\cal{O}}\left(x^{2},\frac{\beta^{2}}{k^{3}}\right)\rightarrow\frac{1}{24k}\left[\left(24q_{0}^{2}-1\right)-\left(\frac{\beta}{k}\right)+{\cal{O}}\left(\frac{\beta^{2}}{k^{2}}\right)\right],$
(3.17)

as $x\rightarrow 0$ (or $y\rightarrow\infty$). Note that
$\zeta\rightarrow\infty$ as $y\rightarrow\infty$, which can be seen clearly
from Eq.(3.11). Thus, to minimize the error control function for very large
values of $y$, we must choose

$q_{0}^{2}=\frac{1}{24}\simeq 4.167\times 10^{-2}.$ (3.18)

Figure 3: Plots of the quantities $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$,
and the error control function ${\cal{T}}$ for $k=5.0$, $\beta=1.0$ and
$q_{0}=1/\sqrt{24}$, for which we have $y_{2}=y_{1}^{*}=1.36944i$.

In Fig. 3 we plot the functions, $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, and
$\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$, together with the error control
function defined by Eqs.(C.3) - (C.5) for $(k,\beta)=(5.0,1.0)$, with $q_{0}$
being given by Eq.(3.18). (Recall
$\beta_{0}\equiv\sqrt{\beta^{2}+q_{0}^{2}}$). From these figures it is clear
that the conditions (2.14) - (2.16) are well satisfied, and the error control
function remains small all the time. In particular, it decreases as $\beta/k$
decreases.

Figure 4: Plots of the mode functions
$\mu_{k}(y),\;\mu^{N}_{k}(y),\mu^{E}_{k}(y)$ and their relative differences
$\delta^{N}(y),\;\delta^{E}(y)$ and $\epsilon(y)$ for $k=5.0$, $\beta=1.0$ and
$q_{0}=1/\sqrt{24}$, for which we have $y_{2}=y_{1}^{*}=1.36944i$.

In Fig. 4, we plot the mode functions $\mu_{k}(y)$, $\mu^{\text{N}}_{k}(y)$,
$\mu^{E}_{k}(y)$, and the relative difference $\delta^{\text{A}}(y)$ defined
by

$\displaystyle\delta^{\text{A}}(y)$ $\displaystyle\equiv$
$\displaystyle\left|\frac{\left|\mu_{k}(y)\right|-\left|\mu^{\text{A}}_{k}(y)\right|}{\mu^{\text{A}}_{k}(y)}\right|,$
(3.19)

where $A=(N,E)$, $\mu_{k}(y)$ denotes the mode function obtained by the UAA
method given by Eq.(3.14), $\mu^{\text{N}}_{k}(y)$ is the numerical solution
obtained by integrating Eq.(2.1) directly with the same initial conditions,
while $\mu^{\text{E}}_{k}(y)$ is the exact solution given by Eq.(Appendix B.
Exact Solutions with the Pöschl-Teller Potential). From these figures we can
see that the maximal errors occur in the region near $y=0$, but the upper
bound is no larger than $0.15\%$ at any given $y$, including the region near
$y\simeq 0$.

It is interesting to note that this analytical approximate solution is only up
to the first-order approximation of the UAA method. With higher order
approximations, the relative errors are even smaller.

To check our numerical solutions, in Fig. 4, we also plot the relative
differences $\epsilon(y)$ between $\mu^{\text{N}}_{k}(y)$ and
$\mu^{\text{E}}_{k}(y)$, defined by

$\displaystyle\epsilon(y)$ $\displaystyle\equiv$
$\displaystyle\left|\frac{\left|\mu^{N}_{k}(y)\right|-\left|\mu^{\text{E}}_{k}(y)\right|}{\mu^{\text{E}}_{k}(y)}\right|.$
(3.20)

From these figures it can be seen that $\epsilon(y)$ is no larger than
$10^{-7}$, and our numerical code is well tested and justified.

It is also interesting to note that the mode functions are oscillating for
$y\lesssim-10$, and these fine features are captured in all three mode
functions, although there are some differences in the details. Again, as shown
by their relative variations, these differences are very small. In addition,
we also consider other choices of $\beta$ and $k$, and find that they all have
similar properties, as long as the condition $k^{2}\gg\beta^{2}$.

Figure 5: Plots of the quantities $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$,
and the error control function ${\cal{T}}$ for $k=5.0$, $\beta=4.9$ and
$q_{0}=1/2$, for which we have $y_{1}=y_{2}^{*}=0.200335i$.

Figure 6: Plots of the mode functions
$\mu_{k}(y),\;\mu^{N}_{k}(y),\mu^{E}_{k}(y)$ and their relative differences
$\delta^{N}(y),\;\delta^{E}(y)$ and $\epsilon(y)$ for $k=5$, $\beta=4.9$ and
$q_{0}=1/2$.

###  III.2 $\beta^{2}\simeq k^{2}$

In this case, depending on $k\gtrsim\beta$ or $k\lesssim\beta$, the function
$g(y)$ has different properties, as shown in Fig. 2. Therefore, in the
following subsections let us consider them separately.

####  III.2.1 $k\gtrsim\beta$

When $k\gtrsim\beta$ the function $g(y)$ is always non-positive for
$y\in(-\infty,\infty)$. Then, from Eqs.(C.3) and (III.1) we find that

$\displaystyle\mathscr{T}(y)\simeq\frac{q_{0}^{2}-1/4}{2\beta}\ln\left(\frac{2}{\epsilon}\right)+\frac{9}{48k}+{\cal{O}}\left(\epsilon\right),$
(3.21)

as $y\rightarrow\infty$, but now with $\epsilon\equiv(k-\beta)/k$. Thus, to
have the error control function be finite at $y=\infty$, now we must set

$q_{0}^{2}=\frac{1}{4},$ (3.22)

instead of the value given by Eq.(3.18) for the case $k\gg\beta^{2}$. In Fig.
5, we plot the quantities $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$,
and the error control function ${\cal{T}}$ for $k=5.0$, $\beta=4.9$ and
$q_{0}=1/2$, for which we have $y_{1}=y_{2}^{*}=0.200335i$. From these figures
we can see clearly that the conditions (2.14) - (2.16) are well satisfied, and
the error control function remains small all the time. Then, the corresponding
quantities $\mu_{k}(y)$, $\mu^{\text{N}}_{k}(y)$, $\mu^{E}_{k}(y)$,
$\delta^{\text{A}}(y)$ and $\epsilon(y)$ are plotted in Fig. 6. From the
curves of $\delta^{\text{N}}(y)$ and $\delta^{\text{E}}(y)$ we can see that
now the errors of the first-order UAA solution are $\leq 4\%$, which are
larger than those of the last subcase. This is mainly because of the fast
oscillations of the solution in the region $y<0$. Therefore, in order to
obtain solutions with high precision, high-order approximations for this case
are needed. However, we do like to note that our numerical solution still
matches to the exact one very well, as shown by the curve of $\epsilon(y)$,
which is no larger than $6.0\times 10^{-6}$.

Figure 7: Plots of the quantities $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$,
and the error control function ${\cal{T}}$ for $k=5.0$, $\beta=5.1$ and
$q_{0}=1/2$, for which we have $y_{1}=-y_{2}=-0.199668$.

####  III.2.2 $k\lesssim\beta$

In this case, we find that

$\zeta_{0}^{2}=\frac{2}{\pi}\left|\int_{y_{1}}^{y_{2}}{\sqrt{g(y)}\;dy}\right|=2\left|k-\beta\right|.$
(3.23)

On the other hand, from Eqs.(II.2.3), (C.3) and (.2) we find that

$\displaystyle\mathscr{T}(y)\simeq\begin{cases}\frac{\zeta(0)\left(6\zeta^{2}_{0}-\zeta^{2}(0)\right)}{12\zeta^{2}_{0}\left(\zeta^{2}_{0}-\zeta^{2}(0)\right)^{3/2}},&y\rightarrow
0,\cr\frac{\pi\left(q_{0}^{2}-1/4\right)}{2\beta},&y\rightarrow
y_{2},\cr\end{cases}$ (3.24)

where $\zeta(0)\equiv\left.\zeta(y)\right|_{y=0}<\zeta_{0}$. Note that in
calculating the error control function near the turning point $y\simeq y_{2}$,
we have used the relation

$\frac{\beta}{k^{2}\sqrt{\beta^{2}-k^{2}}}\left(\beta^{2}x^{2}-k^{2}\right)^{3/2}\simeq\frac{1}{\zeta_{0}}(\zeta_{0}^{2}-\zeta^{2})^{3/2},$
(3.25)

so that the divergence of the second term of $\mathscr{T}_{2}$ cancels exactly
with that of $\mathscr{T}_{3}$. Eq.(3.25) can be obtained directly from the
relation $\sqrt{g}dy=\sqrt{\zeta_{0}^{2}-\zeta^{2}}d\zeta$ for the case $g\geq
0$. Similarly, it can be shown that

$\displaystyle\mathscr{T}(y)\simeq\frac{q_{0}^{2}-1/4}{2\beta}\ln\left(\frac{2}{\epsilon}\right),\;\;\;\;y\rightarrow\infty.$
(3.26)

It is clear that to minimize the errors, in the present case $q_{0}^{2}$ must
be also chosen to be

$q_{0}^{2}=\frac{1}{4},$ (3.27)

as that given by Eq.(3.22). In Fig. 7, we plot the quantities
$\left|{q}/{g}\right|$, $\left|{q(y-y_{1})}/{g}\right|$,
$\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$, and the error control function
${\cal{T}}$ for $k=5.0$, $\beta=5.1$ and $q_{0}=1/2$, for which we have
$y_{1}=-y_{2}=-0.199668$. It is clear that in this case the two turning points
are very close, and the conditions $\left|{q}/{g}\right|\ll 1$ and
$\left|{q(y-y_{1})}/{g}\right|\ll 1$ are violated near these points. But, the
condition $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|\ll 1$ holds near them.
So, the conditions (2.14) - (2.16) are also satisfied, and the error control
function remains small all the time.

Then, the corresponding quantities $\mu_{k}(y)$, $\mu^{\text{N}}_{k}(y)$,
$\mu^{E}_{k}(y)$, $\delta^{\text{A}}(y)$ and $\epsilon(y)$ are plotted in Fig.
8. From the curves of $\delta^{\text{N}}(y)$ and $\delta^{\text{E}}(y)$ we can
see that now the errors of the first-order UAA solution are $\lesssim 10\%$.
Similar to the last subcase, this is mainly because of the fast oscillations
of the solution in the region $y<0$. Therefore, in order to obtain high
precision, high-order approximations for this case are needed, too. In
addition, our numerical solution still matches well to the exact one, as shown
by the curve of $\epsilon(y)$, which is no larger than $2.0\times 10^{-6}$.

Figure 8: Plots of the mode functions
$\mu_{k}(y),\;\mu^{N}_{k}(y),\mu^{E}_{k}(y)$ and their relative differences
$\delta^{N}(y),\;\delta^{E}(y)$ and $\epsilon(y)$ for $k=5$, $\beta=5.1$ and
$q_{0}=1/2$.

Figure 9: Plots of the quantities $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$,
and the error control function ${\cal{T}}$ for $k=0.6$, $\beta=4.0$ and
$q_{0}=1/2$, for which we have $y_{1}=-y_{2}\simeq-2.58459$.

(a) (b)

(c) (d)

(e) (f)

Figure 10: Plots of the mode functions
$\mu_{k}(y),\;\mu^{N}_{k}(y),\mu^{E}_{k}(y)$ and their relative differences
$\delta^{N}(y),\;\delta^{E}(y)$ and $\epsilon(y)$ for $k=0.6$, $\beta=4.0$ and
$q_{0}=1/2$, for which we have $y_{1}=-y_{2}\simeq-2.58459$.

###  III.3 $\beta^{2}\gg k^{2}$

In this case, two real turning points appear, given, respectively by

$\displaystyle y_{1}=-y_{2}=-\cosh^{-1}\left(\frac{\beta^{2}}{k^{2}}\right).$
(3.28)

Then, we find that Eqs.(3.23) and (3.24) still hold in the current case, while
Eq.(3.25) is replaced by

$\displaystyle\mathscr{T}(y)\rightarrow\frac{q_{0}^{2}-1/4}{2\beta}\ln\left(\frac{1+\epsilon}{1-\epsilon}\right)+\frac{4+\epsilon-5\epsilon^{2}}{24k(1-\epsilon^{2})},~{}~{}~{}$
(3.29)

as $y\rightarrow\infty$, but now with $\epsilon\equiv k/\beta$. Combining
Eqs.(3.23), (3.24) and (3.29), we find that currently the proper choice of
$q_{0}$ is still that given by $q_{0}=1/2$, as those given in the last two
subcases.

In Fig. 9, we plot the quantities $\left|{q}/{g}\right|$,
$\left|{q(y-y_{1})}/{g}\right|$, $\left|{q(y-y_{1})(y-y_{2})}/{{{g}}}\right|$,
and the error control function ${\cal{T}}$ for $k=0.6$, $\beta=4.0$ and
$q_{0}=1/2$, for which we have $y_{1}=-y_{2}\simeq-2.58459$. From this figure
we can see that the preconditions (2.14)-(2.16) are well satisfied. Then, to
the first-order approximation of the UAA method, the solution can be
approximated by Eq.(II.2.3), where $\zeta_{0}^{2}$ is given by Eq.(3.23),
$\alpha_{k}$ and $\beta_{k}$ are two integration constants, and $\epsilon_{1}$
and $\epsilon_{2}$ the errors of the corresponding approximate solutions,
whose upper bounds are given by Eqs.(Appendix A. Upper Bounds of Errors) and
(A.2) in Appendix A.

In Fig. 10 (a), we plot out our first-order approximate solution, while Fig.
10 (b) to compare the approximate solution with the exact one, we plot both of
them. In particular, the solid line represents the exact solution, while the
red dotted line the approximate solution. From this figure it can be seen that
except the minimal points, the two solutions match well. However, at these
extreme minimal points, they deviate significantly from each. The causes of
such errors are not clear, and we hope to come back to this issue in another
occasion.

Finally, similar to all other cases, our numerical solution still matches well
to the exact one, as shown by the curve of $\epsilon(y)$, which is no larger
than $8.0\times 10^{-6}$.

##  IV Conclusions

In this paper, we have applied the UAA method to the mode function $\mu_{k}$
with a PT potential, for which it satisfies the second-order differential
equation

$\displaystyle\frac{d^{2}\mu_{k}(y)}{dy^{2}}+\left(k^{2}-\frac{\beta_{0}^{2}}{\cosh
y}\right)\mu_{k}(y)=0,$ (4.1)

where $k$ and $\beta_{0}$ are real constants. In this case, the exact solution
is known and given by Eq.(Appendix B. Exact Solutions with the Pöschl-Teller
Potential). The implementation of the UAA method includes the introduction of
an auxiliary function $q(y)$, which is taken as

$\displaystyle q(y)=\frac{q_{0}^{2}}{\cosh y},$ (4.2)

where $q_{0}$ is a free parameter. Then, we carry out the integration of the
error control function, defined by

$\displaystyle\mathscr{T}(\zeta)$ $\displaystyle\equiv$
$\displaystyle-\int\frac{\psi(\zeta)}{{\left|\dot{y}^{2}g\right|^{1/2}}}d\zeta,$
(4.3)

where

$\displaystyle\psi(\zeta)$ $\displaystyle\equiv$
$\displaystyle\dot{y}^{2}q+\dot{y}^{1/2}\frac{d^{2}}{d\zeta^{2}}\left(\dot{y}^{-1/2}\right),$
$\displaystyle\dot{y}^{2}g$ $\displaystyle=$
$\displaystyle\zeta_{0}^{2}-\zeta^{2}.$ (4.4)

Clearly, the error control function $\mathscr{T}(\zeta)$ will depend on
$q_{0}$. After working out the details, we find that it is convenient to
distinguish the three cases: A) $k^{2}\gg\beta^{2}$, B)
$k^{2}\simeq\beta^{2}$, and C) $k^{2}\ll\beta^{2}$, where
$\beta^{2}\equiv\beta_{0}^{2}-q_{0}^{2}>0$. In particular, in Case A), a
proper choice of $q_{0}$ is $q_{0}=1/\sqrt{24}$, while in Cases B) and C), it
is $q_{0}=1/2$.

Once $q_{0}$ is fixed, the analytical approximate solutions are uniquely
determined by the linear combination of the two parabolic cylindrical
functions $W(\zeta_{0}^{2}/2,\pm\sqrt{2}\zeta)$, as shown by Eq.(II.2.3). In
particular, in Case A) the upper bounds of errors are $\lesssim 0.15\%$, as
shown in Fig. 4. In Case B), two subcases are considered, one with
$k\gtrsim\beta$, and the other with $k\lesssim\beta$. In the first case, the
upper bounds of errors are $\lesssim 4\%$, while in the second case they are
$\lesssim 10\%$, as shown, respectively, by Figs. 6 and 8. In Case C), the
approximate solutions trace also very well to the exact one, except the
minimal points, as shown in Fig. 10. This might be caused by the fact that at
these points the mode function $\mu_{k}$ is almost zero, and very small non-
zero values will cause significantly deviations. We are still working on this
case, and hope to come back to this point in another occasion.

As mentioned in the Introduction, the potentials of the mode functions in both
dressed metric and hybrid approaches can be well modeled by PT potentials.
Therefore, the current analysis on the choice of the function $q(y)$ and the
minimization of the error control function shall shed great light on how to
carry out similar analyses in order to obtain more accurate approximate
solutions in these models. We have been working on it recently, and wish to
report our results soon in another occasion.

In addition, the differential equations for the quasi-normal modes of black
holes usually also take the form of Eq.(2.1) with potentials that have no
singularities 444Recall the inner boundaries of black hole perturbations are
the horizons, at which the potentials are usually finite and non-singular.,
but normally do have turning points BCS09 ; KZ11 . Thus, one can equally apply
the analysis presented here to the studies of quasi-normal modes of black
holes.

## Acknowledgements

RP and AW were partially supported by the US Natural Science Foundation (NSF)
with the Grant No. PHY2308845, and JJM and JS were supported through the
Baylor Physics graduate program. BFL and TZ were supported in part by the
National Key Research and Development Program of China under Grant No.
2020YFC2201503, the National Natural Science Foundation of China under Grant
Nos. 11975203, 11675143, 12205254, 12275238, and 12005186, the Zhejiang
Provincial Natural Science Foundation of China under Grant Nos. LR21A050001
and LY20A050002, and the Fundamental Research Funds for the Provincial
Universities of Zhejiang in China under Grant No. RF-A2019015.

## Appendix A. Upper Bounds of Errors

The upper bounds of the errors $\epsilon_{1}$ and $\epsilon_{2}$ appearing in
Eq.(2.27) are given by

$\displaystyle\frac{|\epsilon_{1}|}{M\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)},\;\frac{|\partial\epsilon_{1}/\partial\zeta|}{\sqrt{2}N\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)}\leq\frac{\kappa}{\lambda_{0}E\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)}\Bigg{\\{}\exp{\Big{(}\lambda\mathscr{V}_{\zeta,\zeta_{2}}(\mathscr{T})\Big{)}}-1\Bigg{\\}}.~{}~{}~{}~{}~{}~{}$
$\displaystyle\frac{|\epsilon_{2}|}{M\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)},\;\frac{|\partial\epsilon_{2}/\partial\zeta|}{\sqrt{2}N\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)}\leq\frac{\kappa
E\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)}{\lambda}\Bigg{\\{}\exp{\Big{(}\lambda_{0}\mathscr{V}_{0,\zeta}(\mathscr{T})\Big{)}}-1\Bigg{\\}},~{}~{}~{}~{}$
(A.1)

where $M\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)$,
$N\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)$, and
$E\left(\frac{1}{2}\lambda\zeta_{0}^{2},\sqrt{2\lambda}\zeta\right)$ are
auxiliary functions of the parabolic cylinder functions defined explicitly in
Olver75 , and 555This corresponds to choosing the function $\Omega(x)$
introduced by Olver in Olver75  as
$\Omega(x)=\sqrt{\left|x^{2}-\zeta_{0}^{2}\right|}$, which satisfies the
requirement $\Omega(x)={\cal{O}}(x)$, as $x\rightarrow\pm\infty$. For more
details, see Olver75 .

$\displaystyle\mathscr{V}_{\zeta_{1},\zeta_{2}}\equiv\int_{\zeta_{1}}^{\zeta_{2}}{\frac{\left|\psi(\zeta)\right|}{\sqrt{|\zeta^{2}-\zeta_{0}^{2}|}}d\zeta},$
(A.2)

is the associated error control function.

## Appendix B. Exact Solutions with the Pöschl-Teller Potential

Let us consider the case with the Pöschl-Teller Potential given by

$\left(\lambda^{2}g+q\right)=-\left(k^{2}-\frac{\beta_{0}^{2}}{\cosh^{2}y}\right).$
(B.1)

Then, introducing the two new variables $x$ and ${\cal{Y}}$ via the relations

$\displaystyle
x=\frac{1}{1+e^{-2y}},\quad{\cal{Y}}(x)=[x(1-x)]^{ik/2}\mu_{k},$ (B.2)

we find that Eq.(2.1) with the above PT potential reads

$\displaystyle
x(1-x)\frac{d^{2}{\cal{Y}}}{dx^{2}}+\left[a_{3}-\left(a_{1}+a_{2}+1\right)x\right]\frac{d{\cal{Y}}}{dx}-a_{1}a_{2}{\cal{Y}}=0,$

where

$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\sqrt{1-4\beta_{0}^{2}})-ik,$ $\displaystyle
a_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1-\sqrt{1-4\beta_{0}^{2}})-ik,$ $\displaystyle
a_{3}$ $\displaystyle=$ $\displaystyle 1-ik.$ (B.4)

Eq.(Appendix B. Exact Solutions with the Pöschl-Teller Potential) is the
standard hypergeometric equation, and has the general solution ZWCKS17

$\displaystyle\mu^{\text{E}}_{k}(\eta)$ $\displaystyle=$ $\displaystyle
a_{k}\left(\frac{x}{1-x}\right)^{ik/2}$
$\displaystyle\times\;_{2}F_{1}(a_{1}-a_{3}+1,a_{2}-a_{3}+1,2-a_{3},x)$
$\displaystyle+\frac{b_{k}}{[x(1-x)]^{ik/2}}\;_{2}F_{1}(a_{1},a_{2},a_{3},x).$

Here ${}_{2}F_{1}(a_{1},a_{2},a_{3},x)$ denotes the hypergeometric function,
and $a_{k}$ and $b_{k}$ are two independent integration constants, and are
uniquely determined by the initial conditions.

## Appendix C. Computing the error control function

In this appendix, we collect some useful formulae for working out the error
control function explicitly. Considering the particular form of the PT
potential, it is easier to compute the error control function by using the new
variable $x=\mathrm{sech}(y)$, thus

$dy=-\frac{\epsilon_{y}~{}dx}{x\sqrt{1-x^{2}}},$ (C.1)

where $\epsilon_{y}$ denotes the sign of $y$. In terms of the new variable,

$\displaystyle q=q^{2}_{0}x^{2},\quad\quad g=\beta^{2}x^{2}-k^{2}.$ (C.2)

To calculate the error control function explicitly, let us consider the cases
$g<0$ and $g>0$ separately.

###  .1 $g<0$

In this case, the error control function is defined by Eq.(II.2.3), which can
be written as

$\mathscr{T}(\zeta)=\mathscr{T}_{1}(\zeta)+\mathscr{T}_{2}(\zeta)+\mathscr{T}_{3}(\zeta),$
(C.3)

where

$\displaystyle\mathscr{T}_{1}$ $\displaystyle\equiv$
$\displaystyle\int\frac{q}{\sqrt{-g}}dy=-q^{2}_{0}\epsilon_{y}\int\frac{xdx}{\sqrt{1-x^{2}}\sqrt{k^{2}-\beta^{2}x^{2}}}=\frac{q^{2}_{0}\epsilon_{y}}{\beta}\ln\left(\frac{\sqrt{1-x^{2}}\beta+\sqrt{k^{2}-\beta^{2}x^{2}}}{\sqrt{|k^{2}-\beta^{2}|}}\right),$
$\displaystyle\mathscr{T}_{2}$ $\displaystyle\equiv$
$\displaystyle\int\left(\frac{5g^{\prime
2}}{16g^{3}}-\frac{g^{{}^{\prime\prime}}}{4g^{2}}\right)\sqrt{-g}dy=\epsilon_{y}\int
dx\left(\frac{5\beta^{4}(x^{3}-x^{5})}{4\sqrt{1-x^{2}}(k^{2}-\beta^{2}x^{2})^{5/2}}+\frac{\beta^{2}(2x-3x^{3})}{2\sqrt{1-x^{2}}(k^{2}-\beta^{2}x^{2})^{3/2}}\right)$
$\displaystyle=$
$\displaystyle\epsilon_{y}\Bigg{\\{}-\frac{1}{4\beta}\ln\left(\frac{\sqrt{1-x^{2}}\beta+\sqrt{k^{2}-x^{2}\beta^{2}}}{\sqrt{|k^{2}-\beta^{2}|}}\right)+\frac{\sqrt{1-x^{2}}A}{12(k^{2}-\beta^{2})(k^{2}-\beta^{2}x^{2})^{3/2}}\Bigg{\\}},$
$\displaystyle\mathscr{T}_{3}$ $\displaystyle\equiv$
$\displaystyle\int^{\zeta}\left\\{\frac{5\zeta_{0}^{2}}{4(v^{2}-\zeta_{0}^{2})^{5/2}}+\frac{3}{4(v^{2}-\zeta_{0}^{2})^{3/2}}\right\\}dv=\frac{\zeta^{3}-6\zeta\zeta^{2}_{0}}{12\zeta^{2}_{0}(\zeta^{2}-\zeta^{2}_{0})^{3/2}},$
(C.4)

where

$A(x)\equiv 3k^{4}+2k^{2}\beta^{2}\left(x^{2}-1\right)-3x^{2}\beta^{4}.$ (C.5)

###  .2 $g>0$

In this case, the error control function is defined by Eq.(II.2.3), which can
be also written as Eq.(C.3), but now with

$\displaystyle\mathscr{T}_{1}$ $\displaystyle\equiv$
$\displaystyle\int\frac{q}{\sqrt{g}}dy=\epsilon_{y}\frac{q^{2}_{0}}{\beta}\arcsin\left(\frac{\beta\sqrt{1-x^{2}}}{\sqrt{\beta^{2}-k^{2}}}\right),$
$\displaystyle\mathscr{T}_{2}$ $\displaystyle\equiv$
$\displaystyle\int\left(-\frac{5g^{\prime
2}}{16g^{3}}+\frac{g^{{}^{\prime\prime}}}{4g^{2}}\right)\sqrt{g}dy=\epsilon_{y}\int
dx\left(\frac{5\beta^{4}(x^{3}-x^{5})}{4\sqrt{1-x^{2}}(\beta^{2}x^{2}-k^{2})^{5/2}}-\frac{\beta^{2}(2x-3x^{3})}{2\sqrt{1-x^{2}}(\beta^{2}x^{2}-k^{2})^{3/2}}\right)$
$\displaystyle=$
$\displaystyle\epsilon_{y}\Bigg{\\{}-\frac{1}{4\beta}\arcsin\left(\frac{\sqrt{1-x^{2}}\beta}{\sqrt{\beta^{2}-k^{2}}}\right)+\frac{\sqrt{1-x^{2}}A}{12(\beta^{2}-k^{2})(\beta^{2}x^{2}-k^{2})^{3/2}}\Bigg{\\}},$
$\displaystyle\mathscr{T}_{3}$ $\displaystyle\equiv$
$\displaystyle\int^{\zeta}\left\\{\frac{5\zeta_{0}^{2}}{4(\zeta_{0}^{2}-v^{2})^{5/2}}-\frac{3}{4(\zeta_{0}^{2}-v^{2})^{3/2}}\right\\}dv=\frac{6\zeta\zeta^{2}_{0}-\zeta^{3}}{12\zeta^{2}_{0}(\zeta^{2}_{0}-\zeta^{2})^{3/2}}.$
(C.6)

## References

  * (1) C. Kiefer, Quantum Gravity, third edition (Oxford Science Publications, Oxford University Press, 2012); H. W. Hamber, Quantum Gravity, the Feynman Path Integral Approach (Springer-Verlag Berlin Heidelberg, 2009).
  * (2) M.B. Green, J.H. Schwarz and E. Witten, Superstring Theory: Vol.1 $\&$ 2, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 1999); J. Polchinski, String Theory, Vol. 1 $\&$ 2 (Cambridge University Press, Cambridge, 2001); C. V. Johson, D-Branes, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2003); K. Becker, M. Becker, and J.H. Schwarz, String Theory and M-Theory (Cambridge University Press, Cambridge, 2007).
  * (3) A. Ashtekar and J. Lewandowski, Background independent quantum gravity: A status report, Class. Quant. Grav. 21, R53 (2004) [arXiv:gr-qc/0404018]; T. Thiemann, Modern Canonical Quantum General Relativity, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2007); C. Rovelli, Quantum gravity (Cambridge University Press, Cambridge, 2008); M. Bojowald, Canonical Gravity and Applications: Cosmology, Black Holes, and Quantum Gravity (Cambridge University Press, Cambridge, 2011); R. Gambini and J. Pullin, A First Course in Loop Quantum Gravity (Oxford University Press, Oxford, 2011); C. Rovelli, and F. Vidotto, Covariant Loop Quantum Gravity (Cambridge University Press, Cambridge, 2015); A. Ashtekar and J. Pullin, Loop Quantum Gravity, the first 30 years (World Scientific, 2017).
  * (4) T. Thiemann, The LQG string—loop quantum gravity quantization of string theory: I. Flat target space, Class. Quamtum Grav. 23 (2006) 1923.
  * (5) R.C. Helling, and P. Giuseppe, String quantization: Fock vs. LQG representations, arXiv:hep-th/0409182 (2004).
  * (6) B. S. De Witt, Quantum Theory of Gravity. I. The Canonical Theory, Phys. Rev. 160 (1967) 1113; J. A Wheeler, Relativity, Groups and Topology, eds C. M. De Witt and J. A. Wheeler, Benjamin, New York, 1968.
  * (7) A. Ashtekar, New variables for classical and quantum gravity, Phys. Rev. Lett. 57, 2244 (1986); New Hamiltonian formulation of general relativity, Phys. Rev. D36, 1587 (1987).
  * (8) M. Bojowald, Loop Quantum Cosmology, Living Rev. Relativ. 11, 4 (2008); Quantum cosmology: a review, Rep. Prog. Phys. 78 (2015) 023901 [arXiv:1501.04899].
  * (9) A. Ashtekar and P. Singh, Loop Quantum Cosmology: A Status Report, Class. Quant. Grav. 28, 213001 (2011) [arXiv:1108.0893].
  * (10) A. Ashtekar and A. Barrau, Loop quantum cosmology: From pre-inflationary dynamics to observations, Class. Quant. Grav. 32, 234001 (2015) [arXiv:1504.07559].
  * (11) P. Singh, special issue editor, Loop Quantum Cosmology, Int. J. Mod. Phys. D25, No. 8 (2016).
  * (12) I. Agullo and P. Singh Loop Quantum Cosmology, in Loop Quantum Gravity: The First 30 Years, Eds: A. Ashtekar, J. Pullin, World Scientific (2017) [arXiv:1612.01236].
  * (13) E. Wilson-Ewing, Testing loop quantum cosmology, Comptes Rendus Physique 18, 207 (2017) [arXiv:1612.04551].
  * (14) B. Elizaga Navascués, G. A. Mena Marugán, Hybrid Loop Quantum Cosmology: An Overview, Front. Astron. Space Sci. 8, 624824 (2021) [arXiv:2011.04559].
  * (15) A. Ashtekar and E. Bianchi, A short review of loop quantum gravity, Rep. Prog. Phys. 84 (2021) 042001 [arXiv:2104.04394].
  * (16) I. Agullo, A. Wang, E. Wilson-Ewing, Loop quantum cosmology: relation between theory and observations, to appear in the “Handbook of Quantum Gravity”, C. Bambi, L. Modesto, I. Shapiro (editors), Springer (2023) [arXiv:2301.10215].
  * (17) B.-F. Li, P. Singh, Loop Quantum Cosmology: Physics of Singularity Resolution and its Implications, to appear in the “Handbook of Quantum Gravity”, C. Bambi, L. Modesto, I. Shapiro (editors), Springer (2023) [arXiv:2304.05426].
  * (18) V. Taveras, Corrections to the Friedmann Equations from LQG for a Universe with a Free Scalar Field, Phys. Rev. D78, 064072 (2008) [arXiv:0807.3325].
  * (19) A. Ashtekar, B. Gupt, Generalized effective description of loop quantum cosmology, Phys. Rev. D92, 084060 (2015) [arXiv:1509.08899]; I. Agullo, A. Ashtekar, B. Gupt, Phenomenology with fluctuating quantum geometries in loop quantum cosmology, Class. Quantum Grav. 34 (2017) 074003 [arXiv:1611.09810].
  * (20) P. Singh, Glimpses of Space-Time Beyond the Singularities Using Supercomputers, Comput. Sci. Eng. 20, 26 (2018) [arXiv:1809.01747]; P. Diener, A. Joe, M. Megevand and P. Singh, Numerical simulations of loop quantum Bianchi-I spacetimes, Class. Quant. Grav. 34, 094004 (2017) [arXiv:1701.05824]; P. Diener, B. Gupt, M. Megevand and P. Singh, Numerical evolution of squeezed and non-Gaussian states in loop quantum cosmology, Class. Quant. Grav. 31, 165006 (2014) [arXiv:1406.1486]; P. Diener, B. Gupt and P. Singh, Numerical simulations of a loop quantum cosmos: robustness of the quantum bounce and the validity of effective dynamics, Class. Quant. Grav. 31, 105015 (2014) [arXiv:1402.6613]; and references therein.
  * (21) W. Kaminski, M. Kolanowski, and J. Lewandowski, Dressed metric predictions revisited, Class. Quantum Grav. 37 (2020) 095001 [arXiv:1912.02556].
  * (22) P. Singh, Are loop quantum cosmos never singular?, Class. Quant. Grav. 26, 125005 (2009) [arXiv:0901.2750]; P. Singh, Curvature invariants, geodesics and the strength of singularities in Bianchi-I loop quantum cosmology, Phys. Rev. D85, 104011 (2012) [arXiv:1112.6391]; S. Saini and P. Singh, Geodesic completeness and the lack of strong singularities in effective loop quantum Kantowski-Sachs spacetime, Class. Quant. Grav. 33, 245019 (2016) [arXiv:1606.04932]; Resolution of strong singularities and geodesic completeness in loop quantum Bianchi-II spacetimes, Class. Quant. Grav. 34, 235006 (2017) [arXiv:1707.08556]; Generic absence of strong singularities in loop quantum Bianchi-IX spacetimes, Class. Quant. Grav. 35 (2018) 065014 [arXiv:1712.09474].
  * (23) P. Singh, K. Vandersloot, and G. V. Vereshchagin, Nonsingular bouncing universes in loop quantum cosmology, Phys. Rev. D74, 043510 (2006) [arXiv:gr-qc/0606032].
  * (24) X. Zhang and Y. Ling, Inflationary universe in loop quantum cosmology, J. Cosmol. Astropart. Phys. 08 (2007) 012 [arXiv:0705.2656].
  * (25) A. Ashtekar and D. Sloan, Loop quantum cosmology and slow roll inflation, Phys. Lett. B694, 108 (2010) [arXiv:0912.4093].
  * (26) A. Ashtekar and D. Sloan, Probability of inflation in loop quantum cosmology, Gen. Relativ. Gravit. 43, 3619 (2011) [arXiv:1103.2475].
  * (27) L. Chen and J.-Y. Zhu, Loop quantum cosmology: The horizon problem and the probability of inflation, Phys. Rev. D92, 084063 (2015) [arXiv:1510.03135].
  * (28) I. Agullo, A. Ashtekar, W. Nelson, Quantum Gravity Extension of the Inflationary Scenario, Phys. Rev. Lett. 109, 251301 (2012).  10.1103/PhysRevLett.109.251301
  * (29) I. Agullo, A. Ashtekar, W. Nelson, Extension of the quantum theory of cosmological perturbations to the Planck era, Phys. Rev. D 87(4), 043507 (2013).  10.1103/PhysRevD.87.043507
  * (30) I. Agullo, A. Ashtekar, W. Nelson, The pre-inflationary dynamics of loop quantum cosmology: confronting quantum gravity with observations, Class. Quant. Grav. 30, 085014 (2013).  10.1088/0264-9381/30/8/085014
  * (31) M. Fernandez-Mendez, G.A. Mena Marugán, J. Olmedo, Hybrid quantization of an inflationary universe, Phys. Rev. D 86, 024003 (2012).  10.1103/PhysRevD.86.024003
  * (32) M. Fernández-Méndez, G.A. Mena Marugán, J. Olmedo, Effective dynamics of scalar perturbations in a flat Friedmann-Robertson-Walker spacetime in loop quantum cosmology, Phys. Rev. D 89(4), 044041 (2014).  10.1103/PhysRevD.89.044041
  * (33) F.B. Martínez, J. Olmedo, Primordial tensor modes of the early Universe, Phys. Rev. D 93(12), 124008 (2016).  10.1103/PhysRevD.93.124008
  * (34) B. Elizaga Navascués, G.A. Mena Marugán, Hybrid Loop Quantum Cosmology: An Overview , Front. Astron. Space Sci. 8, 81 (2021).  10.3389/fspas.2021.624824
  * (35) M. Bojowald, G.M. Hossain, Cosmological vector modes and quantum gravity effects,Class. Quant. Grav. 24, 4801 (2007).  10.1088/0264-9381/24/18/015
  * (36) M. Bojowald, G.M. Hossain, Loop quantum gravity corrections to gravitational wave dispersion , Phys. Rev. D 77, 023508 (2008).  10.1103/PhysRevD.77.023508
  * (37) M. Bojowald, G.M. Hossain, M. Kagan, S. Shankaranarayanan, Anomaly freedom in perturbative loop quantum gravity, Phys. Rev. D 78, 063547 (2008).  10.1103/PhysRevD.78.063547
  * (38) M. Bojowald, G.M. Hossain, M. Kagan, S. Shankaranarayanan, Gauge invariant cosmological perturbation equations with corrections from loop quantum gravity ,Phys. Rev. D 79, 043505 (2009).  10.1103/PhysRevD.79.043505.  [Erratum: Phys.Rev.D 82, 109903 (2010)]
  * (39) E. Wilson-Ewing, Lattice loop quantum cosmology: scalar perturbations ,Class. Quant. Grav. 29, 215013 (2012).  10.1088/0264-9381/29/21/215013
  * (40) E. Wilson-Ewing, Separate universes in loop quantum cosmology: Framework and applications, Int. J. Mod. Phys. D 25(08), 1642002 (2016).  10.1142/S0218271816420025
  * (41) I. Agullo and N.A. Morris, Detailed analysis of the predictions of loop quantum cosmology for the primordial power spectra, Phys. Rev. D92, 124040 (2015).
  * (42) F.W.J. Olver, Asymptotics and Special functions, (AKP Classics, Wellesley, MA 1997).
  * (43) F.W.J. Olver, The asymptotic solution of linear differential equations of the second order in a domain containing one transition point, Philos. Trans. Roy. Soc. London, A249 (1956) 65.
  * (44) F.W.J. Olver, Second-order linear differential equations with two turning points, Philos. Trans. Roy. Soc. London, A278 (1975) 137.
  * (45) A. Wang, Vector and tensor perturbations in Horava-Lifshitz cosmology, Phys. Rev. D82, 124063 (2010) [arXiv:1008.3637].
  * (46) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, Q. Sheng, Constructing analytical solutions of linear perturbations of inflation with modified dispersion relations, Int. J. Mod. Phys. A29 (2014) 1450142 [arXiv:1308.1104].
  * (47) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, Q. Sheng, Inflationary cosmology with nonlinear dispersion relations, Phys. Rev. D89, 043507 (2014) [arXiv:1308.5708].
  * (48) T. Zhu, A. Wang, Gravitational quantum effects in the light of BICEP2 results, Phys. Rev. D90, 027304 (2014) [arXiv:1403.7696].
  * (49) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, Q. Sheng, Gravitational quantum effects on power spectra and spectral indices with higher-order corrections, Phys. Rev. D90, 063503 (2014) [arXiv:1405.5301].
  * (50) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, Q. Sheng, Power spectra and spectral indices of k-inflation: high-order corrections, Phys. Rev. D90, 103517 (2014) [arXiv:1407.8011]
  * (51) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, Q. Sheng, Q. Wu, Detecting quantum gravitational effects of loop quantum cosmology in the early universe?, Astrophy. J. Lett. 807 (2015) L17 [arXiv:1503.06761].
  * (52) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, Q. Sheng, Q. Wu, Scalar and tensor perturbations in loop quantum cosmology: High-order corrections, JCAP 10 (2015) 052 [arXiv:1508.03239].
  * (53) T. Zhu, A. Wang, K. Kirsten, G. Cleaver, Q. Sheng, Q. Wu, Inflationary spectra with inverse-volume corrections in loop quantum cosmology and their observational constraints from Planck 2015 data, JCAP 03 (2016) 046 [arXiv:1510.03855].
  * (54) T. Zhu, A. Wang, K. Kirsten, G. Cleaver, Q. Sheng, High-order Primordial Perturbations with Quantum Gravitational Effects, Phys. Rev. D93, 123525 (2016) [arXiv:1604.05739].
  * (55) Q. Wu, T. Zhu, A. Wang, Primordial Spectra of slow-roll inflation at second-order with the Gauss-Bonnet correction, Phys. Rev. D97, 103502 (2018) [arXiv:1707.08020].
  * (56) J. Qiao, G.-H. Ding, Q. Wu, T. Zhu, A. Wang, Inflationary perturbation spectrum in extended effective field theory of inflation, JCAP 09 (2019) 064 [arXiv:1811.03216].
  * (57) T. Zhu, Q. Wu, A. Wang, An analytical approach to the field amplification and particle production by parametric resonance during inflation and reheating, Phys. Dark Universe 26, 100373 (2019) [arXiv:1811.12612].
  * (58) B.-F. Li, T. Zhu, A. Wang, K. Kirsten, G. Cleaver, Q. Sheng, Preinflationary perturbations from the closed algebra approach in loop quantum cosmology, Phys. Rev. D99, 103536 (2019) [arXiv:1812.11191].
  * (59) G.-H. Ding, J. Qiao, Q. Wu, T. Zhu, A. Wang, Inflationary perturbation spectra at next-to-leading slow-roll order in effective field theory of inflation, Eur. Phys. J. C79, 976 (2019) [arXiv:1907.13108].
  * (60) B.-F. Li, T. Zhu, A. Wang, Langer Modification, Quantization condition and Barrier Penetration in Quantum Mechanics, Universe 6, 90 (2020) [arXiv:1902.09675].
  * (61) J. Qiao, T. Zhu, W. Zhao, A. Wang, Polarized primordial gravitational waves in the ghost-free parity-violating gravity, Phys. Rev. D101, 043528 (2020) [arXiv:1911.01580].
  * (62) T. Zhu, W. Zhao and A. Wang, Polarized primordial gravitational waves in spatial covariant gravities, Phys. Rev. D 107, 024031 (2023) [arXiv:2210.05259 [gr-qc]].
  * (63) T. C. Li, T. Zhu and A. Wang, Power spectra of slow-roll inflation in the consistent $D\rightarrow 4$ Einstein-Gauss-Bonnet gravity, JCAP 05, 006 (2023) [arXiv:2212.08253 [gr-qc]].
  * (64) S. Habib, K. Heitmann, G. Jungman and C. Molina-Paris, The Inflationary perturbation spectrum, Phys. Rev. Lett. 89, 281301 (2002).
  * (65) S. Habib, A. Heinen, K. Heitmann, G. Jungman and C. Molina-Paris, Characterizing inflationary perturbations: The Uniform approximation, Phys. Rev. D70, 083507 (2004).
  * (66) S.-H. Dong, Wave Equations in Higher Dimensions (Springer, New York, 2011).
  * (67) T. Zhu, A. Wang, G. Cleaver, K. Kirsten, and Q. Sheng, Pre-inflationary universe in loop quantum cosmology, Phys. Rev. D96, 083520 (2017).
  * (68) Q. Wu, T. Zhu, and A. Wang, Nonadiabatic evolution of primordial perturbations and non-Gaussinity in hybrid approach of loop quantum cosmology, Phys. Rev. D98, 103528 (2018).
  * (69) F.W. J. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clark, NIST Handbook of Mathematical Functions (Cambridge University Press, Cambridge, 2010).
  * (70) E. Berti, V. Cardoso and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quant. Grav. 26, 163001 (2009) [arXiv:0905.2975 [gr-qc]].
  * (71) R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83, 793 (2011) [arXiv:1102.4014 [gr-qc]].