1 Introduction

The milestone in multi-messenger astronomy is marked by the detection of the gravitational waves from the merger of binary black holes [1]. Some typical characteristics of these black holes, such as their unexpectedly large masses and relatively small spins, are not consistent with usual astrophysical black holes, but are more possible to be of primordial origin [2,3,4]. Primordial black holes (PBHs) are theorized to form before the Big Bang nucleosynthesis, so they are a very powerful probe in the early Universe [5]. For instance, they could explain the strong absorption trough found in the 21-cm global spectrum [6, 7], and serve as the seeds of the supermassive black holes in galactic centers [8], supported by the recent observation from the James Webb Space Telescope on the massive galaxies at high redshifts [9]. Meanwhile, the first-order scalar perturbations that generate PBHs can also produce the second-order scalar-induced gravitational waves (SIGWs) [10,11,12,13,14,15,16], such as the possible nHz gravitational wave background discovered recently [17,18,19,20,21,22,23,24,25,26,27]. More important, PBHs are a natural and promising candidate of dark matter (DM) [28].

Unlike astrophysical black holes, the masses of PBHs can range from the Planck mass to supermassive scale (\(10^{-38}\)\(10^{10}~M_\odot ,\) with \(M_\odot =1.989\times 10^{30}~\text {kg}\) being the solar mass) [29]. This versatility enables PBHs to offer valuable insights into different cosmic conundrums [30]. The PBHs with mass \(M<5\times 10^{-19}~M_\odot \) have already evaporated because of the Hawking radiation, changing the background intensities of various cosmic rays [31]. However, those with mass \(M>5\times 10^{-19}~M_\odot \) can still exist today, acting as a stable and pressureless candidate of DM. The PBH abundance \(f_{\textrm{PBH}}\) is defined as its proportion in DM today. If \(f_{\textrm{PBH}} > rsim 0.1,\) the PBHs can be considered as an effective candidate of DM; if \(f_{\textrm{PBH}}\ll 10^{-3},\) its possibility as DM can be safely excluded in the relevant mass range. Currently, according to various constraints on the upper bounds of \(f_{\textrm{PBH}}\) in different mass ranges, there still remains an open mass window from \(10^{-17}~M_\odot \) (the asteroid mass range) to \(10^{-13}~M_\odot \) (the sub-lunar mass range) that is possible to compose all DM with PBHs [28, 29, 31, 32].

The PBH abundance \(f_{\textrm{PBH}}\) can be calculated from the power spectrum \({{\mathscr {P}}}_{{\mathscr {R}}}(k)\) of the primordial curvature perturbation \({{\mathscr {R}}}\) [33,34,35]. On large scales (e.g., the pivot scale \(k_* =0.05~{\textrm{Mpc}}^{-1}\) in the Planck satellite experiment), \({{\mathscr {P}}}_{{\mathscr {R}}}\) is precisely measured by the anisotropies of the cosmic microwave background (CMB), with an amplitude of \(2.10\times 10^{-9}\) [36]. However, to produce abundant PBHs, \({{\mathscr {R}}}\) must be large enough, so that \({{\mathscr {P}}}_{{\mathscr {R}}}\) is significantly enhanced up to \({{\mathscr {O}}}(10^{-2})\) on small scales. In the usual single-field slow-roll (SR) inflation models, such a huge enhancement of \({{\mathscr {P}}}_{{\mathscr {R}}}\) is almost impossible. Therefore, the SR conditions must be violated on small scales, and various ultraslow-roll (USR) inflation models are thus proposed [37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75]. In these models, there is always a plateau or a (near-)inflection point on the inflaton potential, where the inflaton rolls down extremely slowly, amplifying \({{\mathscr {P}}}_{{\mathscr {R}}}\) accordingly.

Albeit the amplitude and shape of \({{\mathscr {P}}}_{{\mathscr {R}}}\) can be straightforwardly obtained in various USR inflation models by numerical methods, the physical comprehension behind seems still indistinct. Previously, in Refs. [51, 76, 77], the spectral index of the steepest growth of \({{\mathscr {P}}}_{{\mathscr {R}}}\) was analytically investigated, and in Ref. [78], the authors raised an approximation method to calculate the Fourier mode \({\mathcal {R}}_k\) of the primordial curvature perturbation. However, to our knowledge, these analyses are still not sufficient in some aspects. For example, there is always a sharp dip in \({{\mathscr {P}}}_{{\mathscr {R}}}\) before its steep growth. In some early literature [79,80,81,82,83,84,85,86,87], the location of the dip was calculated numerically, but the underlying physical explanation is lacking yet. Although the detail of this dip does not affect the PBH abundance or SIGW spectrum, as they are relevant mainly to the peak of \({{\mathscr {P}}}_{{\mathscr {R}}},\) we still hope to answer whether the minimum of \({{\mathscr {P}}}_{{\mathscr {R}}}\) is zero analytically.

Consequently, this paper is dedicated to study the evolution of \({{\mathscr {R}}}_k\) and the resulting \({{\mathscr {P}}}_{{\mathscr {R}}}\) from a more theoretical perspective, and is a natural extension of our previous numerical results in Ref. [87]. We will provide an analytical approximation to explore the USR inflation, which will save us away from the mathematical tediousness and extract the physical essence to the most extent. Following our early works in Refs. [58, 59, 73, 87], by considering an antisymmetric perturbation on the background inflaton potential, inflation can be led into the USR stage smoothly. Our main improvements in this work are threefold. First, by simplifying the two parameters \(\varepsilon \) and \(\eta \) in the SR and USR stages separately, we are able to obtain the analytical solutions of \({{\mathscr {R}}}_k\) and its time derivative \({{\mathscr {R}}}_{k,N}.\) From their asymptotic forms, we achieve \(|{{\mathscr {R}}}_k|,\) \(\theta _{k,N},\) and \(\varphi _{k,N}\) (with \(\theta _{k}\) and \(\varphi _{k}\) being the arguments of \({{\mathscr {R}}}_k\) and \({{\mathscr {R}}}_{k,N}\)) and find that all these quantities possess concise exponential forms, perfectly in agreement with the numerical results in Ref. [87]. Second, we predict that, besides the revolving evolution of \({{\mathscr {R}}}_k\) around the origin in the complex plane when \(k\gg He^N,\) there appears an interesting linear evolution of \({{\mathscr {R}}}_k\) towards or away from the origin if \({{\mathscr {R}}}_k\) crosses the horizon around the start of the USR stage, naturally explaining the sharp dip and the peak in \({{\mathscr {P}}}_{{\mathscr {R}}}.\) Third, we analytically study the minimum of the dip in \({{\mathscr {P}}}_{{\mathscr {R}}}\) and prove that it cannot reach zero exactly, which is also seldom mentioned before. Altogether, we wish to give a whole picture and thorough understanding of the physical essence in the complicated evolution of \({{\mathscr {R}}}_k\) and the relevant \({{\mathscr {P}}}_{{\mathscr {R}}},\) which will be beneficial to the model building of the USR inflation and help our research on PBH and gravitational wave physics in future.

This paper is organized as follows. The basic equations of motion for \({{\mathscr {R}}}_k\) and \({{\mathscr {R}}}_{k,N}\) are presented in Sect. 2. Then, we show our analytical approximation in Sect. 3 and discuss the asymptotic solutions of \(|{{\mathscr {R}}}_k|,\) \(|{{\mathscr {R}}}_k|_{,N},\) \(\theta _{k,N},\) and \(\varphi _{k,N}\) in the sub- and super-horizon limits, respectively. In Sect. 4, the evolution of \({{\mathscr {R}}}_k\) and the resulting \({{\mathscr {P}}}_{{\mathscr {R}}}\) for five typical scales with different times of horizon-exit are systematically investigated. We conclude in Sect. 5. We work in the natural system of units and set \(c=\hbar =k_\textrm{B}=1.\)

2 Basic equations

In this section, we show the equation of motion for the primordial curvature perturbation \({{\mathscr {R}}}_k\) and its relation to the power spectrum \({{\mathscr {P}}}_{{\mathscr {R}}}(k).\) Two important variables \(\theta _k\) and \(\varphi _k\) are introduced, with their evolutions calculated in detail.

2.1 USR inflation

We start from the single-field inflation model, with the corresponding action being

$$\begin{aligned} S=\int {\textrm{d}}^{4} x\,\sqrt{-g}\left[\frac{m_{\textrm{P}}^{2}}{2}R-\frac{1}{2}\partial _{\mu }\phi \partial ^{\mu }\phi -V(\phi )\right], \end{aligned}$$

where \(m_\textrm{P}=1/\sqrt{8\pi G}\) is the reduced Planck mass, R is the Ricci scalar, \(\phi \) is the inflaton field, and \(V(\phi )\) is its potential. The inflaton potential V can be further decomposed into its background \(V_{\textrm{b}}(\phi )\) and a perturbation \(\delta V(\phi )\) on it. In this work, we choose \(V_{\textrm{b}}\) as the Kachru–Kallosh–Linde–Trivedi potential [88],

$$\begin{aligned} V_{\textrm{b}}(\phi )=V_0\frac{\phi ^2}{\phi ^2+(m_\textrm{P}/2)^2}. \end{aligned}$$

Meanwhile, we consider a perturbation \(\delta V\) on \(V_{\textrm{b}},\)

$$\begin{aligned} \delta V(\phi )=-A\frac{\phi -\phi _0}{1+ ( \phi -\phi _0)^2/ ( 2\sigma ^2)}, \end{aligned}$$

where three model parameters A\(\phi _0,\) and \(\sigma \) characterize the amplitude, position, and width of \(\delta V,\) respectively. (In principle, the specific form of \(\delta V\) is not unique, and other kinds of \(\delta V\) can be found in Refs. [58, 59, 73].) Hence, \(\delta V\) is antisymmetric around \(\phi _0,\) so it can be smoothly connected on \(V_{\textrm{b}}\) on both sides of \(\phi _0.\) In addition, we demand \(A=V_{{\textrm{b}},\phi }(\phi _0)\) for simplicity, in order to create a perfect plateau around \(\phi _0,\) leading inflation into the USR stage. Following Ref. [59], we choose the model parameters as \(V_0/m_\textrm{P}^4=10^{-10},\) \(\phi _0/m_\textrm{P}=1.31,\) and \(\sigma /m_\textrm{P}=0.0831881,\) such that there can be PBHs with mass \(M\approx 10^{-17}~M_\odot \) and abundance \(f_{\textrm{PBH}}\approx 0.1.\) The inflaton potential \(V(\phi )\) is shown in Fig. 1. Moreover, the initial conditions for inflation are set to be \(\phi /m_\textrm{P}=3.30\) and \(\phi _{,N}/m_\textrm{P}=-0.0137,\) so as to satisfy the Planck constraints on the power spectrum and the tensor-to-scalar ratio on the CMB pivot scale \(k_*\) [36].

Fig. 1
figure 1

The inflaton potential \(V(\phi )\) with the model parameters \(V_0/m_\textrm{P}^4=10^{-10},\) \(\phi _0/m_\textrm{P}=1.31,\) and \(\sigma /m_\textrm{P}=0.0831881.\) The antisymmetric perturbation \(\delta V\) creates a plateau on V and thus leads inflation into the USR stage, producing the PBHs with mass \(M\approx 10^{-17}~M_\odot \) and abundance \(f_{\textrm{PBH}}\approx 0.1.\) The CMB pivot scale corresponds to \(\phi _{\textrm{CMB}}/m_\textrm{P}=2.933,\) which is far away from the USR regime

In cosmic inflation, it is more convenient to utilize the number of e-folds N as the time variable, defined as \(\textrm{d}N=H\,\textrm{d}t=\textrm{d}\ln a,\) where t is the cosmic time, \(a(t)=e^N\) is the scale factor, and \(H(t)=\dot{a}/a\) is the Hubble expansion rate. Furthermore, to characterize the motion of the inflaton, two useful parameters are introduced as

$$\begin{aligned} \varepsilon&=-\frac{\dot{H}}{H^2}=\frac{\phi _{,N}^2}{2m_\textrm{P}^2}, \end{aligned}$$
(1)
$$\begin{aligned} \eta&=-\frac{\ddot{\phi }}{H{\dot{\phi }}}=\frac{\phi _{,N}^2}{2m_\textrm{P}^2}-\frac{\phi _{,NN}}{\phi _{,N}}. \end{aligned}$$
(2)

In terms of these parameters, the Friedmann equation for cosmic expansion becomes

$$\begin{aligned} H^2=\frac{V}{(3-\varepsilon )m^2_\textrm{P}}, \end{aligned}$$
(3)

and the Klein–Gordon equation for the inflaton field can be reexpressed as

$$\begin{aligned} \phi _{,NN}+(3-\varepsilon )\phi _{,N}+\frac{1}{H^2}V_{,\phi }=0. \end{aligned}$$
(4)

In the usual SR inflation, both \(\varepsilon \) and \(|\eta |\) are always much smaller than 1 and are thus named as the SR parameters. However, in the USR stage, \(\varepsilon \) drops almost exponentially and \(|\eta |\) may even approach \({{\mathscr {O}}}(1),\) so the SR conditions are broken. Thus, the starting and ending points of the USR stage can be determined by \(\eta (N_\mathrm{{s}})=\eta (N_\mathrm{{e}})=0.\) With the model parameters listed above, we have \(N_{\textrm{s}}=56.81\) and \(N_{\textrm{e}}=60.93.\) The evolutions of \(\varepsilon \) and \(|\eta |\) in both SR and USR stages are shown in Fig. 2. Moreover, the number of e-folds when the relevant scale crosses the horizon is denoted as \(N_{\textrm{out}}\) (i.e., when \(k=He^{N_{\textrm{out}}}\)), which is an important quantity for the subsequent discussions in Sect. 4.

Fig. 2
figure 2

The evolutions of the parameters \(\varepsilon \) and \(|\eta |\) as a function of the number of e-folds N. The starting and ending points of the USR stage are determined by \(\eta (N_\mathrm{{s}})=\eta (N_\mathrm{{e}})=0,\) with \(N_{\textrm{s}}=56.81\) and \(N_{\textrm{e}}=60.93.\) During the USR stage, \(\varepsilon \) decreases nearly exponentially, but \(|\eta |\) increases significantly and maintains a value around 3 for a sufficiently long period

2.2 \({{\mathscr {R}}}_k\) and \({{\mathscr {P}}}_{{\mathscr {R}}}\)

Now, we move on to the perturbations on the background spacetime. In Newtonian gauge, the perturbed metric reads

$$\begin{aligned} \textrm{d}s^2=-(1+2\Psi )\,\textrm{d}t^2+a^2(1-2\Psi )\delta _{ij}\,\textrm{d}x^i\textrm{d}x^j, \end{aligned}$$

where \(\Psi \) is the scalar perturbation. (Here, we have neglected the vector perturbation, tensor perturbation, and anisotropic stress, as they are irrelevant to PBH formation.) A more convenient and gauge-invariant scalar perturbation is the primordial curvature perturbation \({{\mathscr {R}}}\) defined as

$$\begin{aligned} {{\mathscr {R}}}=\Psi +\frac{H}{{\dot{\phi }}}\,\delta \phi =\Psi +\frac{\delta \phi }{\phi _{,N}}, \end{aligned}$$

and the equation of motion for its Fourier mode \({{\mathscr {R}}}_k\) is the Mukhanov–Sasaki (MS) equation [89, 90],

$$\begin{aligned} {{\mathscr {R}}}_{k,NN}+(3+\varepsilon -2\eta ){{\mathscr {R}}}_{k,N}+\left(\frac{k}{He^{N}}\right)^2{{\mathscr {R}}}_k=0. \end{aligned}$$
(5)

The initial conditions for the MS equation can be fixed by the Bunch–Davies vacuum [91]. In the very early Universe, both \(\epsilon \) and \(|\eta |\) are very small, and the asymptotic solution of \({{\mathscr {R}}}_k\) is

$$\begin{aligned} {{\mathscr {R}}}_k=\frac{e^{ik/(He^N)}}{z\sqrt{2k}}\left(1+\frac{iHe^N}{k}\right), \end{aligned}$$
(6)

where \(z=\phi _{,N}e^N\) (i.e., \(z^2=2m_\textrm{P}^2a^2\varepsilon \)). At the initial time, we have \(k\gg He^N,\) so from Eq. (6), we obtain

$$\begin{aligned} {{\mathscr {R}}}_k\bigg |_{N\rightarrow N_\mathrm{{ini}}}&=\frac{e^{ik/(He^N)}}{z\sqrt{2k}}, \end{aligned}$$
(7)
$$\begin{aligned} {{\mathscr {R}}}_{k,N}\bigg |_{N\rightarrow N_\mathrm{{ini}}}&=-\frac{e^{ik/(He^N)}}{z\sqrt{2k}}\left(\frac{z_{,N}}{z}+\frac{ik}{He^N}\right), \end{aligned}$$
(8)

where \(N_\mathrm{{ini}}\) is the initial number of e-folds.

Furthermore, the PBH abundance and the SIGW spectrum can be calculated from the power spectrum of \({{\mathscr {R}}},\) which corresponds to the two-point correlation function of \({{\mathscr {R}}}_k\) in the Fourier space. The statistical information and observational significance of \({{\mathscr {R}}}_k\) is encoded in the dimensionless power spectrum \(\mathcal{P}_{{\mathscr {R}}}(k)\) as

$$\begin{aligned} \left.{{\mathscr {P}}}_{{\mathscr {R}}}(k)=\frac{k^3}{2\pi ^2}|{{\mathscr {R}}}_{k}|^2\right|_{k\ll aH}. \end{aligned}$$

We should point out that, in the usual SR inflation, \({{\mathscr {R}}}_k\) becomes almost frozen when the relevant scale crosses the horizon, and \({{\mathscr {P}}}_{{\mathscr {R}}}\) can be effectively evaluated at the epoch of horizon-exit. Nonetheless, in the USR inflation, \({{\mathscr {R}}}_k\) can still evolve significantly after the horizon-exit, so \({{\mathscr {P}}}_{{\mathscr {R}}}\) must be determined at the end of inflation.

Generally speaking, both \({{\mathscr {R}}}_k\) and \({{\mathscr {P}}}_{{\mathscr {R}}}\) can be obtained by numerically solving Eqs. (1)–(5), but this is always a time-consuming task. Consequently, one of the aims of this work is to explore an appropriate and analytical approximation to calculate \({{\mathscr {R}}}_k\) and \({{\mathscr {P}}}_{{\mathscr {R}}}.\) Only in this way can we extract the physical essence in the complicated evolution of \({{\mathscr {R}}}_k\) and the characteristic shape of \({{\mathscr {P}}}_{\mathcal {R}}\) in the USR inflation, without getting lost in the tedious calculation program.

2.3 \(\theta _k\) and \(\varphi _k\)

As the MS equation is complex, it is more convenient to decompose \({{\mathscr {R}}}_k\) as

$$\begin{aligned} {{\mathscr {R}}}_k=|{{\mathscr {R}}}_k| e^{-i\theta _k}, \end{aligned}$$
(9)

where \(|{{\mathscr {R}}}_k|\) and \(\theta _k\) are the modulus and argument of \({{\mathscr {R}}}_k.\) Another advantage of this decomposition, rather than decomposing \({{\mathscr {R}}}_k\) into the real and imaginary parts, is that we will observe interesting revolving and linear evolutions of \({{\mathscr {R}}}_k\) in the complex plane, to be shown in Sect. 4.

Substituting Eq. (9) into Eq. (5), we arrive at the equations of motion for \(|{{\mathscr {R}}}_k|\) and \(\theta _k,\)

$$\begin{aligned}&|{{\mathscr {R}}}_k|_{,NN}+(3+\varepsilon -2\eta )|{{\mathscr {R}}}_k|_{,N} \nonumber \\&\quad +\left(\frac{k}{He^{N}}\right)^2\left[1-\left(\frac{He^{N}\theta _{k,N}}{k}\right)^2\right]|{{\mathscr {R}}}_k|=0, \end{aligned}$$
(10)
$$\begin{aligned}&\theta _{k,NN}+\left(3+\varepsilon -2\eta +2\frac{|{{\mathscr {R}}}_k|_{,N}}{|{{\mathscr {R}}}_k|}\right)\theta _{k,N}=0. \end{aligned}$$
(11)

Furthermore, from Eq. (7), the initial conditions for \(|{{\mathscr {R}}}_k|\) read

$$\begin{aligned} |{{\mathscr {R}}}_k|\bigg |_{N\rightarrow N_\mathrm{{ini}}}&=-\frac{1}{z\sqrt{2k}}, \end{aligned}$$
(12)
$$\begin{aligned} |{{\mathscr {R}}}_k|_{,N}\bigg |_{N\rightarrow N_\mathrm{{ini}}}&=\frac{z_{,N}}{z^2\sqrt{2k}}. \end{aligned}$$
(13)

When the inflaton \(\phi \) rolls down from its potential V,  we have \(z=\phi _{,N}e^N<0.\) Since \(|{{\mathscr {R}}}_k|\) is positive, we have assigned a negative sign in Eq. (12).

Below, for our analytical purpose, we provide two mathematical tricks to simplify the calculations of \(|{{\mathscr {R}}}_k|\) and \(\theta _k.\) First, from Eqs. (6) and (9), we can obtain the evolution of the argument \(\theta _k\) as [57]

$$\begin{aligned} \theta _{k,N}=\frac{1}{2z^2 He^N |{{\mathscr {R}}}_k|^2}, \end{aligned}$$
(14)

and the detailed derivation can be found in Appendix A. Hence, \(\theta _{k,N}\) is positive-definite, so \(\theta _k\) is a monotonic function.

Second, from Eq. (9), we have \({{\mathscr {R}}}_{k,N}=|{{\mathscr {R}}}_k|_{,N}e^{-i\theta _k}-i|{{\mathscr {R}}}_k|\theta _{k,N}e^{-i\theta _k},\) so the evolution of the modulus \(|{{\mathscr {R}}}_k|\) is

$$\begin{aligned} |{{\mathscr {R}}}_k|_{,N}=\sqrt{|{{\mathscr {R}}}_{k,N}|^2-(|{{\mathscr {R}}}_k|\theta _{k,N})^2}. \end{aligned}$$
(15)

By this means, we are able to achieve \(|{{\mathscr {R}}}_k|\) and \(\theta _k\) more quickly, as both Eqs. (14) and (15) are first-order differential equations, in which \({{\mathscr {R}}}_{k}\) and \({{\mathscr {R}}}_{k,N}\) can be easily obtained from the MS equation via the analytical approximation in Sect. 3. In contrast, it is impossible to solve Eq. (10) analytically, due to the complicity in the prefactor of \(|{{\mathscr {R}}}_k|.\)

Similarly, we decompose the derivative of \({{\mathscr {R}}}_k\) as

$$\begin{aligned} {{\mathscr {R}}}_{k,N}=|{{\mathscr {R}}}_{k,N}| e^{-i\varphi _k}, \end{aligned}$$

where the argument \(\varphi _k\) can be obtained from Eq. (9) as

$$\begin{aligned} \varphi _k =\theta _k+\arctan \left(\frac{|{{\mathscr {R}}}_k|}{|{{\mathscr {R}}}_k|_{,N}}\theta _{k,N}\right). \end{aligned}$$
(16)

Substituting Eq. (14) into Eq. (16), we arrive at the evolution of \(\varphi _k\) as

$$\begin{aligned} \varphi _{k,N}=\frac{1}{2z^2He^N|{{\mathscr {R}}}_k|^2}-\frac{(2z^2He^N|{{\mathscr {R}}}_k||{{\mathscr {R}}}_k|_{,N})_{,N}}{(2z^2He^N|{{\mathscr {R}}}_k||{{\mathscr {R}}}_k|_{,N})^2+1}. \end{aligned}$$
(17)

Different from \(\theta _{k,N}\) in Eq. (14), \(\varphi _{k,N}\) depends not only on \(|{{\mathscr {R}}}_k|,\) but also on \(|{{\mathscr {R}}}_k|_{,N}.\)

In summary, Eqs. (9), (14), (15), and (17) summarize the evolutions of the moduli and arguments of \({{\mathscr {R}}}_k\) and \({{\mathscr {R}}}_{k,N},\) which are the starting point of our following discussions.

3 Analytical approximation

In this section, we propose our analytical approximation to calculate \(|{{\mathscr {R}}}_k|\) and \(|{{\mathscr {R}}}_k|_{,N}\) in two different cases: \(k\gg He^N\) and \(k\ll He^N,\) in which the scale of \({{\mathscr {R}}}_k\) is sub- and super-horizon, respectively. Next, we bring the results into Eqs. (14) and (17) to obtain \(\theta _{k,N}\) and \(\varphi _{k,N}.\) Our main results are summarized in Tables 1 and 2.

We start our analytical approximation from the two parameters \(\varepsilon \) and \(\eta .\) First, from Fig. 2, in the SR stage when \(N<N_{\textrm{s}},\) \(\varepsilon \) is nearly invariant and \(\eta \) is very small, so we approximately have

$$\begin{aligned} \varepsilon \approx \varepsilon (N_{\textrm{s}})=\varepsilon _{\textrm{s}}, \quad \eta \approx 0. \end{aligned}$$

Second, in the USR stage when \(N_{\textrm{s}}<N<N_{\textrm{e}},\) as the inflaton potential is extremely flat, from Eqs. (2) and (4), we arrive at \(\eta \approx 3.\) Then, from Eqs. (1) and (2), \(\eta \) can be reexpressed as \(\eta =\varepsilon -{\varepsilon _{,N}}/(2\varepsilon ),\) so we approximately obtain

$$\begin{aligned} \varepsilon \approx \varepsilon _{\textrm{s}} e^{-6(N-N_{\textrm{s}})}, \quad \eta \approx 3. \end{aligned}$$

Third, when \(N>N_{\textrm{e}},\) the inflation dynamics significantly depends on the specific inflaton potential, so it is hard to have a general approximation method any more, and we will not consider it in this work.

Before exploring Eq. (5), we still need to make two assumptions. First, \(\eta \) and H are assumed to be constant. Second, \(\varepsilon \) is neglected when it is added to \(\eta .\) These assumptions are valid in both SR and USR stages, as can be evidently seen from Fig. 2. With these preparations, Eq. (5) can be exactly solved as a function of N as

$$\begin{aligned} {{\mathscr {R}}}_k&= e^{-\left(\frac{3}{2}-\eta \right)N} \nonumber \\&\quad \times \left[A J_{\frac{3}{2}-\eta }\left(\frac{k}{He^N}\right)+B J_{-\frac{3}{2}+\eta }\left(\frac{k}{He^N}\right)\right], \end{aligned}$$
(18)
$$\begin{aligned} {{\mathscr {R}}}_{k,N}&=e^{-\left(\frac{5}{2}-\eta \right)N}\frac{k}{H} \nonumber \\&\quad \times \left[-A J_{\frac{1}{2}-\eta }\left(\frac{k}{He^N}\right) +B J_{-\frac{1}{2}+\eta }\left(\frac{k}{He^N}\right)\right], \end{aligned}$$
(19)

where \(J_\alpha (x)\) is the Bessel function of the first kind, and the coefficients A and B are the functions of kH,  and \(\eta ,\) but are independent of N. Also, we should mention that both A and B are complex numbers, since \({{\mathscr {R}}}_k\) itself is complex.

3.1 \(k\gg He^N\)

We first consider the sub-horizon case with \(k\gg He^N.\) Under this circumstance, our analytical approximation is relatively simple, because \(J_{\alpha }(x)\sim 1/\sqrt{x}\) when \(x=k/(He^N)\gg 1,\) irrespective of the value of \(\alpha .\) As a result, from Eqs. (18) and (19), we easily obtain

$$\begin{aligned} |{{\mathscr {R}}}_k|&\approx c_1 e^{-(1-\eta )N}\sim {\left\{ \begin{array}{ll} e^{-N} &{} (\hbox {when}\ \eta \approx 0), \\ e^{2N} &{} (\hbox {when}\ \eta \approx 3), \end{array}\right. } \end{aligned}$$
(20)
$$\begin{aligned} |{{\mathscr {R}}}_{k,N}|&\approx c_2 e^{-(2-\eta )N}\sim {\left\{ \begin{array}{ll} e^{-2N} &{} (\hbox {when}\ \eta \approx 0), \\ e^{N} &{} (\hbox {when}\ \eta \approx 3), \end{array}\right. } \end{aligned}$$
(21)

where the coefficients \(c_1\) and \(c_2\) can be obtained from A and B,  but the explicit expressions are unnecessary here. We observe that the time-dependence of \({{\mathscr {R}}}_k\) and \({{\mathscr {R}}}_{k,N}\) is rather simple. More interpretation of these results will be presented in Sect. 4.

Table 1 All the possible results of \(|{{\mathscr {R}}}_k|,\) \(\theta _{k,N},\) \(|{{\mathscr {R}}}_k|_{,N},\) and \(\varphi _{k,N}\) in the sub-horizon case with \(k\gg He^N\) via our analytical approximation. The parameter \(\eta \approx 0\) or 3 corresponds to the SR or USR stage, respectively. From the behaviors of \(|{{\mathscr {R}}}_k|\) and \(\theta _{k,N},\) \({{\mathscr {R}}}_k\) either revolves towards the origin at a decelerated rate in the SR stage, or away from the origin at an accelerated rate in the USR stage
Table 2 Same as Table 1, but in the super-horizon case with \(k\ll He^N.\) The current situation is more complicated, due to the competition between the two terms in Eqs. (25) and (26), so all the possible four combinations of the leading-order terms are considered

Now, we discuss the SR and USR stages separately. In the SR case, \(\eta \approx 0,\) so from Eqs. (20) and (21), we have

$$\begin{aligned} |{{\mathscr {R}}}_k|\sim e^{-N}, \quad |{{\mathscr {R}}}_{k,N}|\sim e^{-2N}. \end{aligned}$$
(22)

Substituting \(|{{\mathscr {R}}}_k|\) into Eq. (14) and considering \(\varepsilon \approx \varepsilon _{\textrm{s}},\) we obtain the evolution of \(\theta _k\) as

$$\begin{aligned} \theta _{k,N}\sim \frac{1}{e^{2N}\varepsilon _{\textrm{s}} He^N e^{-2N}}\sim e^{-N}. \end{aligned}$$
(23)

Then, substituting Eqs. (22) and (23) into Eq. (15), we have

$$\begin{aligned} |{{\mathscr {R}}}_k|_{,N}\sim e^{-2N}. \end{aligned}$$

For the evolution of \(\varphi _k,\) substituting Eqs. (22) and (23) into Eq. (17), we obtain

$$\begin{aligned} \varphi _{k, N}\sim \frac{1}{e^{2N}\varepsilon _{\textrm{s}} He^N e^{-2N}}\sim e^{-N}. \end{aligned}$$
Fig. 3
figure 3

The power spectrum \({{\mathscr {P}}}_{{{\mathscr {R}}}}(k)\) in the USR inflation, with the model parameters in Sect. 2.1. On large scales (e.g., the CMB pivot scale \(k_*=0.05~\text {Mpc}^{-1}\) denoted with an asterisk), \({{\mathscr {P}}}_{{{\mathscr {R}}}}\) is nearly scale-invariant with an amplitude of \(2.10\times 10^{-9}\) [36]. On small scales, \({{\mathscr {P}}}_{{{\mathscr {R}}}}\) can be enhanced up to \({{\mathscr {O}}}(10^{-2})\) to produce abundant PBHs and SIGWs. Five typical scales with different \(N_{\textrm{out}}=50,\) 52.89, 55, 58.84, and 62 are marked as dots with different colors, corresponding to the nearly scale-invariant region, sharp dip, steep growth, peak, and falling stage, respectively. The minimum of \({{\mathscr {P}}}_{{\mathscr {R}}}\) in the sharp dip is merely \(1.61\times 10^{-14},\) but not zero

Fig. 4
figure 4

The evolutions of the primordial curvature perturbation \({{\mathscr {R}}}_k\) and the two arguments \(\theta _k\) and \(\varphi _k\) with \(N_{\textrm{out}}=50.\) Before \(N_{\textrm{out}},\) we have \(k\gg He^N\) and \(\eta \approx 0,\) and \(k^3|{{\mathscr {R}}}_k|^2/(2\pi ^2)\sim e^{-2N},\) \(\theta _{k,N}\sim e^{-N},\) and \(\varphi _{k,N}\sim e^{-N},\) so the slopes of the blue lines in ac are \(-2,\) \(-1,\) and \(-1.\) Similarly, when \(N_{\textrm{out}}<N<N_{\textrm{s}},\) we have \(k\ll He^N\) and \(\eta \approx 0,\) so the slopes of the red lines are 0, \(-3,\) and \(-1.\) When \(N_{\textrm{s}}<N<N_{\textrm{e}},\) we have \(k\ll He^N\) and \(\eta \approx 3,\) so the slopes of the purple lines are 0, 3, and \(-5.\) All these values match the numerical results (the dashed lines) very well. Moreover, the evolution of \({{\mathscr {R}}}_k\) from \(N=47\) to 52 (\(\varDelta N=0.01\)) is shown in the complex plane in d, in which \({{\mathscr {R}}}_k\) first revolves around the origin at a decelerated rate and finally stops near it. This means that, for the large-scale curvature perturbation that crosses the horizon much earlier before the USR stage, \(|{{\mathscr {R}}}_k|\) becomes almost constant once \(N>N_{\textrm{out}},\) and the influence from the USR inflation is negligible

Fig. 5
figure 5

Same as Fig. 4, but with \(N_{\textrm{out}}=52.89.\) Before the end of the USR stage, the evolutions of \(|{{\mathscr {R}}}_k|,\) \(\theta _{k},\) and \(\varphi _{k}\) basically resemble those in Fig. 4. However, around \(N_{\textrm{e}},\) \({{\mathscr {R}}}_k\) is no longer frozen but decreases dramatically as \(|{{\mathscr {R}}}_k|\sim e^{-3N},\) and \(\varphi _{k,N}\) is extremely small at the same time, so there appears a new linear evolution of \({{\mathscr {R}}}_k\) after the revolving process (not shown in the figure) towards but not through origin, as shown in d with \(N=56.8\) to 62.9 (\(\varDelta N=0.1\)). This linear evolution naturally explains the sharp dip in \({{\mathscr {P}}}_{{\mathscr {R}}},\) and the minimum of the dip is not exactly zero from both theoretical and numerical points of view

Fig. 6
figure 6

Same as Fig. 4, but with \(N_{\textrm{out}}=55.\) When \(N< N_{\textrm{out}},\) the evolutions of \(|{{\mathscr {R}}}_k|,\) \(\theta _{k},\) and \(\varphi _{k}\) are analogous to those in Fig. 4, as shown in the blue lines in ac. During the USR stage, we have \(k\ll He^N\) and \(\eta \approx 3,\) but there are two competitive terms in \(|{{\mathscr {R}}}_k|\) from Eq. (25). As a result, in the early phase, we have \(|{{\mathscr {R}}}_k|\sim {\mathrm{const.}},\) \(\theta _{k,N}\sim e^{3N},\) and \(\varphi _{k,N}\sim e^{-5N},\) so the slopes of the red lines are 0, 3,  and \(-5.\) Similarly, in the late phase, the slopes of the purple lines are 6, \(-3,\) and \(-3.\) Hence, \({{\mathscr {R}}}_k\) shows a linear evolution away from the origin and finally stops at a distant point in the complex plane, explaining the steep growth in \({{\mathscr {P}}}_{{\mathscr {R}}},\) as shown in d with \(N=52\) to 65 (\(\varDelta N= 0.01\)). Moreover, at \(N\approx 59.15,\) a sharp decrease appears in \(|{{\mathscr {R}}}_k|\) in a, accompanied by a spike in \(\theta _{k,N}\) in b, as a natural consequence from Eq. (14)

Fig. 7
figure 7

Same as Fig. 4, but with \(N_{\textrm{out}}=58.84.\) Here, our analytical approximation applies only to the situation when \(N<N_{\textrm{s}}\) with \(k\gg He^N.\) As shown in the blue lines in ac, the evolutions of \(|{{\mathscr {R}}}_k|,\) \(\theta _k,\) and \(\varphi _k\) are similar to the previous cases. From the dashed lines in a and c, the rapid increase of \(|{{\mathscr {R}}}_k|\) and decrease of \(\varphi _{k,N}\) almost maintain as before in Fig. 6a and c, so the linear evolution of \({{\mathscr {R}}}_k\) still exists. In d, \({{\mathscr {R}}}_k\) linearly moves away from the origin and eventually stops at a rather distant point, resulting in the peak in \({{\mathscr {P}}}_{{\mathscr {R}}}\) (the previous revolving process is not shown), with \(N=58.8\) to 62 \((\varDelta N=0.05)\)

Fig. 8
figure 8

Same as Fig. 4, but with \(N_{\textrm{out}}=62.\) Because \(N_{\textrm{out}}>N_{\textrm{e}},\) we only need to consider the simple case with \(k\gg He^N.\) When \(N<N_{\textrm{s}},\) we have \(\eta \approx 0,\) so the slopes of the blue lines in ac are \(-2,\) \(-1,\) and \(-1.\) Next, when \(N_{\textrm{s}}<N<N_{\textrm{e}},\) we have \(\eta \approx 3,\) so the slopes of the red lines are 4,  \(-1,\) and \(-1.\) Moreover, \({{\mathscr {R}}}_k\) first revolves towards the origin and then away from it, as can be seen in d with \(N=59\) to 64 \((\varDelta N=0.02).\) Although \(|{{\mathscr {R}}}_k|\) changes dramatically around the USR stage, the value of \(\varphi _{k,N}\) is still large enough, so \({{\mathscr {R}}}_k\) cannot show a linear evolution in the complex plane any more

Similarly, in the USR case, \(\eta \approx 3,\) so

$$\begin{aligned} |{{\mathscr {R}}}_k|\sim e^{2N}, \quad |{{\mathscr {R}}}_{k,N}|\sim e^N. \end{aligned}$$

Taking into account \(\varepsilon \approx \varepsilon _{\textrm{s}} e^{-6(N-N_{\textrm{s}})},\) we have

$$\begin{aligned} \theta _{k,N}\sim \frac{1}{e^{2N}\varepsilon _{\textrm{s}}e^{-6(N-N_{\textrm{s}})} He^N e^{4N}}\sim e^{-N}. \end{aligned}$$
(24)

Thus, \(\theta _{k,N}\) has the same time-dependence as that in Eq. (23). However, this happens only for the sub-horizon case. In the super-horizon case, \(\theta _{k,N}\) may have totally different results in the SR and USR cases, to be shown in Sect. 3.2. Furthermore, by the same means, we obtain

$$\begin{aligned} |{{\mathscr {R}}}_k|_{,N}\sim e^{N}, \end{aligned}$$

and

$$\begin{aligned} \varphi _{k, N}\sim \frac{1}{e^{2N}\varepsilon _{\textrm{s}}e^{-6(N-N_{\textrm{s}})} He^N e^{4N}}\sim e^{-N}. \end{aligned}$$

It is interesting to find that \(\varphi _{k,N}\) behaves the same in both SR and USR stages, but this is understandable. In the numerator of the second term in Eq. (17), \(2z^2He^N|{{\mathscr {R}}}_k||{{\mathscr {R}}}_k|_{,N}\) is actually constant no matter \(\eta \approx 0\) or 3. Thus, its derivative with respect to N vanishes, so \(\varphi _{k,N}\approx \theta _{k,N}\sim e^{-N}.\)

All the above results are summarized in Table 1, indicating that, when a scale is sub-horizon, \(|{{\mathscr {R}}}_k|\) is a monotonic function of N. Therefore, \({{\mathscr {R}}}_k\) either revolves towards the origin in the complex plane at a decelerated rate in the SR stage, or away from the origin at an accelerated rate in the USR stage. These conclusions will be compared with the numerical results in Sect. 4.Footnote 1

3.2 \(k\ll He^N\)

Next, we turn to the super-horizon case with \(k\ll He^N.\) The current situation is a little bit complicated than the sub-horizon case, and the essential reason is that, when \(x=k/(He^N)\ll 1,\) the asymptotic form of the Bessel function \(J_\alpha (x)\sim x^\alpha \) depends on \(\alpha .\) Consequently, from Eqs. (18) and (19), we obtain

$$\begin{aligned} |{{\mathscr {R}}}_k|&\approx \big |c_3 e^{-(3-2\eta )N}+c_4\big | \nonumber \\&={\left\{ \begin{array}{ll} \big |c_3e^{-3N}+c_4\big | &{} (\text{ when } \eta \approx 0), \\ \big |c_3e^{3N}+c_4\big | &{} (\text{ when } \eta \approx 3), \end{array}\right. } \end{aligned}$$
(25)
$$\begin{aligned} |{{\mathscr {R}}}_{k,N}|&\approx \big |c_5e^{-(3-2\eta )N}+c_{6}e^{-2N}\big | \nonumber \\&={\left\{ \begin{array}{ll} \big |c_5e^{-3N}+c_{6}e^{-2N}\big | &{} (\text{ when } \eta \approx 0), \\ \big |c_5e^{3N}+c_{6}e^{-2N}\big | &{} (\text{ when } \eta \approx 3), \end{array}\right. } \end{aligned}$$
(26)

where the coefficients \(c_3\) to \(c_6\) can also be obtained from A and B. Now, the situation becomes more intractable, as there are two competitive terms in both \(|{{\mathscr {R}}}_k|\) and \(|{{\mathscr {R}}}_k|_{,N}.\) Hence, we must be very careful to determine the leading-order one, and we discuss this issue in detail below.

We start from the SR stage with \(\varepsilon \approx \varepsilon _{\textrm{s}}\) and \(\eta \approx 0.\) If the leading-order terms are

$$\begin{aligned} |{{\mathscr {R}}}_k|\approx c_3 e^{-3N}, \quad |{{\mathscr {R}}}_{k,N}|\approx c_5e^{-3N}, \end{aligned}$$

from Eqs. (14), (15), and (17), by the same method in Sect. 3.1, we obtain

$$\begin{aligned} \theta _{k,N}\sim e^{3N}, \quad |{{\mathscr {R}}}_k|_{,N}\sim e^{-3N}, \quad \varphi _{k, N}\sim e^{3N}. \end{aligned}$$

Similarly, if the leading-order terms are

$$\begin{aligned} |{{\mathscr {R}}}_k|\approx c_4, \quad |{{\mathscr {R}}}_{k,N}|\approx c_6 e^{-2N}, \end{aligned}$$

the corresponding results read

$$\begin{aligned} \theta _{k,N}\sim e^{-3N}, \quad |{{\mathscr {R}}}_k|_{,N}\sim e^{-2N}, \quad \varphi _{k, N}\sim e^{-N}. \end{aligned}$$

Next, we move on to the USR stage with \(\varepsilon \approx \varepsilon _{\textrm{s}} e^{-6(N-N_{\textrm{s}})}\) and \(\eta \approx 3.\) If the leading-order terms are

$$\begin{aligned} |{{\mathscr {R}}}_k|\approx c_3 e^{3N}, \quad |{{\mathscr {R}}}_{k,N}|\approx c_5 e^{3N}, \end{aligned}$$

we obtain

$$\begin{aligned} \theta _{k,N}\sim e^{-3N}, \quad |{{\mathscr {R}}}_k|_{,N}\sim e^{3N}, \quad \varphi _{k, N}\sim e^{-3N}. \end{aligned}$$

Similarly, if the leading-order terms are

$$\begin{aligned} |{{\mathscr {R}}}_k|\approx c_4, \quad |{{\mathscr {R}}}_{k,N}|\approx c_6 e^{-2N}, \end{aligned}$$

we have

$$\begin{aligned} \theta _{k,N}\sim e^{3N}, \quad |{{\mathscr {R}}}_k|_{,N}\sim e^{-2N}, \quad \varphi _{k,N}\sim e^{-5N}. \end{aligned}$$

Again, all the above results are summarized in Table 2, according to the increasing rate of \(|{{\mathscr {R}}}_k|.\) Because of the competition between the two terms in Eqs. (25) and (26), there are altogether four different combinations at present. From the last two rows, we find that, during the USR stage, \(\varphi _{k,N}\) decreases dramatically, so \(\varphi _k\) will reach a constant very soon. Under these circumstances, \({{\mathscr {R}}}_k\) will have tiny angular velocity in the complex plane, so there will appear a linear evolution of \({{\mathscr {R}}}_k,\) to be shown in Figs. 56 and 7.

Last, we briefly comment that our analytical approximation mainly applies to the two limiting cases with \(k\gg He^N\) and \(k\ll He^N,\) but is unfortunately not suitable for the intermediate region around \(N\sim N_{\textrm{out}}\) with \(k\sim He^N,\) because at this time, it is difficult to determine the value of \([k/(He^N)]^2\) in front of \({{\mathscr {R}}}_k\) in Eq. (5).

4 Evolution of the primordial curvature perturbation

In this section, we utilize the results in Sect. 3 to provide a whole picture of the evolution of the primordial curvature perturbation \({{\mathscr {R}}}_k\) in the USR inflation. Meanwhile, we compare our analytical approximation with the numerical results. For such comparison, five typical scales are chosen, for which the numbers of e-folds when they cross the horizon are \(N_{\textrm{out}}=50,\) 52.89, 55, 58.84, and 62, respectively. The reason for choosing these numbers is that they correspond to five typical positions in the power spectrum \({{\mathscr {P}}}_{{\mathscr {R}}}(k)\): the nearly scale-invariant region, sharp dip, steep growth, peak, and falling stage, respectively, as shown in Fig. 3.

Below, the evolution of \({{\mathscr {R}}}_k\) with different \(N_{\textrm{out}}\) is systematically explored in order and compared with the numerical results (i.e., the lower right panels in Figs. 4567 and 8) in our previous work in Ref. [87].

First, we start from the case with \(N_{\textrm{out}}=50,\) in which the relevant scale crosses the horizon much earlier than \(N_{\textrm{s}},\) and our results are shown in Fig. 4. Before \(N_{\textrm{out}},\) we have \(k\gg He^N\) and \(\eta \approx 0.\) Therefore, \(|{{\mathscr {R}}}_k|\sim e^{-N}\) [or equivalently, \(k^3|{{\mathscr {R}}}_k|^2/(2\pi ^2)\sim e^{-2N}\)], \(\theta _{k,N}\sim e^{-N},\) and \(\varphi _{k,N}\sim e^{-N}\) (the first row in Table 1), as depicted in the blue lines in Fig. 4a–c. However, after \(N_{\textrm{out}},\) we should be cautious and distinguish the SR and USR stages. When \(N_{\textrm{out}}<N<N_{\textrm{s}}\) (i.e., the SR stage), we have \(k\ll He^N\) and \(\eta \approx 0,\) so \(|{{\mathscr {R}}}_k|\sim \mathrm{const.},\) \(\theta _{k,N}\sim e^{-3N},\) and \(\varphi _{k,N}\sim e^{-N}\) (the second row in Table 2). However, when \(N_{\textrm{s}}<N<N_{\textrm{e}}\) (i.e., the USR stage), we have \(k\ll He^N\) and \(\eta \approx 3,\) so \(|{{\mathscr {R}}}_k|\sim \mathrm{const.},\) \(\theta _{k,N}\sim e^{3N},\) and \(\varphi _{k,N}\sim e^{-5N}\) (the third row in Table 2). These behaviors can also be found in the red and purple lines in Fig. 4a–c. Altogether, we observe that our analytical approximation fits the numerical results (the dashed lines) perfectly in all these situations. Furthermore, the evolution of \({{\mathscr {R}}}_k\) from \(N=47\) to 52 is exhibited in the complex plane in Fig. 4d, with the difference of the number of e-folds between two adjacent dots set to be \(\varDelta N=0.01.\) From the results of \(|{{\mathscr {R}}}_k|\) and \(\theta _{k,N},\) it is easy to find that \({{\mathscr {R}}}_k\) revolves clockwise towards the origin and eventually stops near it, where \(|{{\mathscr {R}}}_k|\) tends to be constant. The interval between the dots monotonically decreases, meaning that \({{\mathscr {R}}}_{k}\) decelerates gradually. All these behaviors indicate that the USR stage has negligible effect on such large-scale curvature perturbation, and \(|{{\mathscr {R}}}_k|\) becomes almost frozen once \(N>N_{\textrm{out}},\) which is reflected in the nearly scale-invariant region in \({{\mathscr {P}}}_{{\mathscr {R}}}.\)

Second, we study the scale with \(N_{\textrm{out}}=52.89,\) and our results are shown in Fig. 5. In fact, most behaviors of \(|{{\mathscr {R}}}_k|,\) \(\theta _{k,N},\) and \(\varphi _{k,N}\) (e.g., the slopes of the blue, red, and purple lines in Fig. 5a–c) remain the same as those in the previous case. The unique distinction is that, from the first row in Table 2, around the end of the USR stage, \(|{{\mathscr {R}}}_k|\sim e^{-3N}\) decreases dramatically (the green line in Fig. 5a), and \(\varphi _{k,N}\) is extremely small. Hence, there appears a new linear evolution of \({{\mathscr {R}}}_k\) towards but not through the origin subsequent to the revolving process, as shown in Fig. 5d, with \(N = 56.8\) to 62.9. This interesting point naturally explains the sharp dip in \({{\mathscr {P}}}_{{\mathscr {R}}}.\) Moreover, we should emphasize that this dip actually cannot reach zero, consistent with the numerical result (the dashed line in Fig. 5a). If there were \({{\mathscr {R}}}_k\) that eventually stops at the origin, we would have \(|{{\mathscr {R}}}_k|=0\) and \(\theta _{k,N}=0\) simultaneously. However, from Eq. (14), this is strictly forbidden, so the minimum of \({{\mathscr {P}}}_{{\mathscr {R}}}\) can never be zero, contrary to the naive inference in Ref. [92].

Third, we consider the scale with \(N_{\textrm{out}}=55,\) which is before but near \(N_{\textrm{s}}.\) When \(N<N_{\textrm{out}},\) the situation is analogous to the former two cases, as shown in the blue lines in Fig. 6a–c. However, when \(N_{\textrm{s}}<N<N_{\textrm{e}},\) we have \(k\ll He^N\) and \(\eta \approx 3,\) so we have to face the two competitive terms in Eq. (25) and determine the leading-order one in advance. As a result, in the early phase of the USR stage, we follow the third row in Table 2, so \(|{{\mathscr {R}}}_k|\sim \mathrm{const.},\) \(\theta _{k,N}\sim e^{3N},\) and \(\varphi _{k,N}\sim e^{-5N},\) corresponding to the red lines. In contrast, in the late phase of the USR stage, we arrive at the last row in Table 2, so \(|{{\mathscr {R}}}_k|\sim e^{3N},\) \(\theta _{k,N}\sim e^{-3N},\) and \(\varphi _{k,N}\sim e^{-3N},\) as shown with the purple lines. As \(|{{\mathscr {R}}}_k|\) is significantly amplified and \(\varphi _{k,N}\) is tiny simultaneously, \({{\mathscr {R}}}_k\) shows a linear evolution away from the origin and finally terminates at a remote point in the complex plane, inducing the steep growth in \(\mathcal{P}_{{\mathscr {R}}},\) as shown Fig. 6d with \(N=52\) to 65. In addition, the competition between the two terms in Eq. (25) causes a sharp decrease of \(|{{\mathscr {R}}}_k|\) at \(N\approx 59.15.\) From Eq. (14), this decrease also explains a spike in \(\theta _{k,N},\) as can be seen in Fig. 6a and b (such divergent behaviors are beyond our analytical approximation).

Fourth, the peak of the power spectrum \({{\mathscr {P}}}_{{\mathscr {R}}}\) corresponds to the scale with \(N_{\textrm{out}}=58.84,\) and our analytical approximation is partly suitable to this special case. Now, \(N_{\textrm{out}}\) is almost in the middle between \(N_{\textrm{s}}=56.81\) and \(N_{\textrm{e}}=60.93,\) but the approximation is inapplicable when \(k\sim He^N.\) Therefore, we can only find the similar behaviors of \(|{{\mathscr {R}}}_k|,\) \(\theta _{k,N},\) and \(\varphi _{k,N}\) to those in the previous three cases when \(N<N_{\textrm{s}}\) (not \(N_{\textrm{out}},\) since it is larger than \(N_{\textrm{s}}\) now). During the USR stage, the numerical method is still needed. From the dashed lines in Fig. 7a and c, the rapid increase of \(|{{\mathscr {R}}}_k|\) and decrease of \(\varphi _{k,N}\) almost maintain as before in Fig. 6a and c, so the linear evolution of \({{\mathscr {R}}}_k\) still exists, resulting in the peak in \({{\mathscr {P}}}_{{\mathscr {R}}}\) eventually, as shown in Fig. 7d with \(N=58.8\) to 62.

Last, we finish our discussion with \(N_{\textrm{out}}=62.\) Because \(N_{\textrm{out}}>N_{\textrm{e}}\) at present, we only need to take into account the simple case with \(k\gg He^N.\) Therefore, in the SR stage, we have \(\eta \approx 0,\) so \(|{{\mathscr {R}}}_k|\sim e^{-N},\) \(\theta _{k,N}\sim e^{-N},\) and \(\varphi _{k,N}\sim e^{-N}\) (the first row in Table 1), as shown in the blue lines in Fig. 8a–c. Next, during the USR stage, we have \(\eta \approx 3,\) so \(|{{\mathscr {R}}}_k|\sim e^{2N},\) \(\theta _{k,N}\sim e^{-N},\) and \(\varphi _{k,N}\sim e^{-N}\) (the second row in Table 1), as shown in the red lines. Since \(|{{\mathscr {R}}}_k|\) first decreases and then increases in the current situation, \({{\mathscr {R}}}_k\) first revolves towards the origin and then away from it, as can be seen in Fig. 8d with \(N=59\) to 64. Besides, there is no more linear evolution of \({{\mathscr {R}}}_k,\) which can also be understood from the behavior of \(\varphi _k,\) since \(\varphi _{k,N}\) is still large enough in Fig. 8c.

In summary, we compare our analytical approximation of the evolution of \({{\mathscr {R}}}_k\) to the numerical calculation in our previous work in Ref. [87]. Except the epoch around \(N_{\textrm{out}}\) with \(k\sim He^N,\) the theoretical and numerical results match each other very well. We admit that there is indeed some situations where our approximation is insufficient. These intricacies mainly root from the competition between the two terms in Eqs. (25) and (26), because it is not always easy to determine the leading-order one between them, so numerical calculation is still necessary there. Last, we claim that there are revolving processes of \({{\mathscr {R}}}_k\) in all the five typical cases, but the subsequent linear evolutions only exist for the three intermediate scales, and the relevant criterion should be that \(|{{\mathscr {R}}}_k|\) changes dramatically and \(\varphi _{k,N}\) is tiny simultaneously.

5 Conclusion

The USR inflation is receiving increasing research interest in recent years, in which the primordial curvature perturbation \({{\mathscr {R}}}_k\) behaves completely different in contrast to that in the SR case. The most remarkable character is that \({{\mathscr {R}}}_k\) can still significantly increase after the horizon-exit and thus greatly enhance the power spectrum \({{\mathscr {P}}}_{{\mathscr {R}}}\) on small scales. As a result, PBHs at certain mass windows can be produced in abundance, serving as an effective candidate of DM. Meanwhile, SIGWs at certain frequencies can also be intense enough to be discovered by current and future gravitational wave detectors.

In this paper, we provide an analytical approximation to systematically investigate the evolution of \({{\mathscr {R}}}_k\) and the relevant \({{\mathscr {P}}}_{{\mathscr {R}}}\) in the USR inflation, which is usually a highly numerical and time-consuming task. The asymptotic solutions of the moduli and arguments of \({{\mathscr {R}}}_k\) and \({{\mathscr {R}}}_{k,N}\) are obtained in both sub- and super-horizon limits with \(k\gg He^N\) and \(k\ll He^N,\) perfectly matching the numerical results in our previous work in Ref. [87]. Moreover, five typical scales with different \(N_{\textrm{out}}\) are studied in order, corresponding to five different positions in \({{\mathscr {P}}}_{{\mathscr {R}}}.\) In all, we hope to construct a framework to facilitate the analytical calculation of \({\mathcal {R}}_k.\) Our basic conclusions are summarized as follows.

First, we obtain the approximate forms of the parameters \(\varepsilon \) and \(\eta \) in both SR and USR stages. By this means, all the analytical solutions of \(|{{\mathscr {R}}}_k|,\) \(|{{\mathscr {R}}}_k|_{,N},\) \(\theta _{k,N},\) and \(\varphi _{k,N}\) can be expressed in simple exponential forms. However, due to the different asymptotic expansion of the Bessel function, the solutions with \(k\ll He^N\) are more complicated than those with \(k\gg He^N,\) because we have to determine the leading-order term of \({{\mathscr {R}}}_k\) in advance. All the possibilities are summarized in Tables 1 and 2.

Second, as for the evolution of \({{\mathscr {R}}}_k\) for the five typical scales with different \(N_{\textrm{out}},\) the numerical results validate our analytical approximation where it is applicable, as shown in Figs. 4567 and 8. The evolution of \({{\mathscr {R}}}_k\) can basically be divided into two types: the revolving and linear processes. The former exists in all the five cases. At early times when \(k\gg He^N,\) in the first four cases with \(N_{\textrm{out}}<N_{\textrm{e}},\) we have \(|{{\mathscr {R}}}_k|\sim e^{-N}\) and \(\theta _{k,N}\sim e^{-N},\) so \({\mathcal {R}}_k\) presents a revolving motion towards the origin in the complex plane with decreasing angular velocity, as shown in Figs. 4d, 5d, 6d, and 7d. However, in the last case with \(N_{\textrm{out}}>N_{\textrm{e}},\) the USR stage causes \(|{{\mathscr {R}}}_k|\sim e^{2N}\) and \(\theta _{k,N}\sim e^{-N},\) so \({{\mathscr {R}}}_k\) revolves away from the origin, as shown in Fig. 8d.

Third, for the three intermediate scales with \(N_{\textrm{out}}=52.89,\) 55, and 58.84,  besides the revolving motion, a linear evolution of \({{\mathscr {R}}}_k\) follows afterward. In these situations, \(|{{\mathscr {R}}}_k|\) changes violently, but at the same time \(\varphi _{k,N}\) is extremely small, so \({{\mathscr {R}}}_k\) nearly evolves along a straight line in the complex plane. When \(N_{\textrm{out}}=52.89,\) \({{\mathscr {R}}}_k\) moves towards but not through the origin in Fig. 5d, inducing the sharp dip in \({{\mathscr {P}}}_{{\mathscr {R}}}\) but with a nonvanishing minimum. On the contrary, when \(N_{\textrm{out}}=55\) and 58.84,  \({{\mathscr {R}}}_k\) departs away from the origin in Figs. 6d and 7d, explaining the steep growth and the peak in \({{\mathscr {P}}}_{{\mathscr {R}}}.\) However, if \(N_{\textrm{out}}\) is too small or too big, \(\varphi _{k,N}\) remains large enough during the USR stage, so there is only revolving process of \({{\mathscr {R}}}_k\) without the subsequent linear evolution.

Altogether, by exploring the evolution of \({{\mathscr {R}}}_k\) via our analytical approximation, we wish to provide a whole picture and thorough understanding of the primordial curvature perturbation and the power spectrum in the USR inflation from a theoretical perspective. Our work will be helpful to the model building of the USR inflation and contribute our knowledge of PBH and gravitational wave physics.