\section{Introduction}
\label{sec:Introduction}
Complex oscillations are ubiquitous in neuronal data, where cells respond to sufficiently strong input stimuli not only by emitting action potentials or spikes but also by displaying more complicated electrical oscillations such as \textit{bursting}, which corresponds to an alternation between epochs of quiescent (subthreshold) activity and groups of consecutive spikes referred to as bursts.
A large repertoire of bursting patterns are observed in neuronal experimental recordings~\cite{alving1968spontaneous,adams1985generation,butera1995analysis,cook1984electrical,jian2004novel,pinault1992voltage} as well as in other excitable cells~\cite{atwater1980nature}. Mathematical models of such complex neuronal oscillatory behaviour often possess multiple timescales, which can be explicit or not, and many different types of experimentally-observed bursting oscillations have been classified over the past few decades~\cite{rinzel1987formal,izhikevich2000neural,golubitsky01} using tools from bifurcation theory and singularity theory; see below.

Biologists and modellers are interested in understanding possible changes (upon input) between different neuronal activity regimes, in particular how neurons may transition from spiking to bursting activity. One scenario of such transition involves special solutions called \textit{torus canards}; they have been evidenced and studied in both biophysical~\cite{kramer2008new} and phenomenological~\cite{benes2011elementary} models. According to this scenario, the membrane potential response to stimuli evolves from repetitive spiking (regular large-amplitude oscillations) to amplitude-modulated spiking, where a second frequency emerges in the solution profile, to then bursting activity. Torus canards are organising this transition, which occurs in a very narrow range of input parameter values.

Mapping the transition between bursting and tonic firing, which may carry
different information content~\cite{latorre2006,welday2011}, is a key question in neuroscience and it has
partially been addressed in single cell recordings\cite{ranck1973studies,
mason1988passive, connors1990intrinsic}, and investigated in a
FitzHugh–Nagumo–Rinzel type model \cite{wojcik2011}. As hinted at above, a first classification of bursting behaviour was proposed by Rinzel~\cite{rinzel1987formal} and later extended by Izhikevich
\cite{izhikevich2000neural}. These classifications hinge upon the slow-fast structure
of the models supporting bursting, and upon the bifurcation structure of certain
associated singular systems, which we now introduce.

In this paper, we will study \textit{elliptic bursting}~\cite{azad2010,ermentrout2005,lajoie2011,su2003effects}, which occurs
minimally in systems of $2$-fast $1$-slow variables. These models are written in
the \emph{fast-time parametrisation}
\begin{align}
  x' & =  f (x,y,\mu,\varepsilon), && y' = g(x,y,\mu,\varepsilon), && \mu' = \varepsilon h(x,y,\mu,\varepsilon), && t \in \mathbb{R}^+,
\label{eq:fastTime} \\
\intertext{or in an equivalent slow-time parametrisation}
  \varepsilon \dot x & =  f (x,y,\mu,\varepsilon), && \varepsilon \dot y = g(x,y,\mu,\varepsilon), && \dot \mu = h(x,y,\mu,\varepsilon), && \tau \in \mathbb{R}^+,
\label{eq:slowTime}
\end{align}
where $x,y$ are fast variables, $\mu$ is the slow variable, $0 < \varepsilon \ll 1$ is a
small parameter indicating the time scale separation, and in which prime and overdot
denote differentiations with respect to the fast time $t$ and to the slow time $\tau
= \varepsilon t$, respectively. The functions $f,g,h$ are supposed to be sufficiently
regular.

As in many bursting patterns, the timescale difference in elliptic bursting systems
manifests itself as a composition of fast oscillatory components coupled to a slow
process driving the fast variables alternatively from quasi-stationary to
quasi-periodic dynamics~\cite{rinzel1987formal}. Intuitively, this means that one
uses both scalings presented above to understand bursting solutions.
Despite being equivalent to each other for
$\varepsilon\neq 0$, systems~\cref{eq:fastTime} and~\cref{eq:slowTime} converge to two
non-equivalent singular limits for $\varepsilon=0$:
\begin{align}
  x' & = f(x,y,\mu,0),  & & y' = g(x,y,\mu,0) & & \mu' = 0,                 & & t \in \mathbb{R}^+,  \label{eq:fastSubs}\\
   0 & = f(x,y,\mu,0),  & & 0  = g(x,y,\mu,0) & & \dot \mu = h(x,y,\mu,0),  & & \tau \in \mathbb{R}^+, \label{eq:slowSubs}
\end{align}
which are referred to as \textit{fast} and \textit{slow
subsystems}, respectively.

An important object in slow-fast dynamical systems is the \textit{critical manifold}
$S_0$
\[
  S_0:=\big\{(x,y,\mu)\in \mathbb{R}^3 \colon 0=f(x,y,\mu,0), \quad 0=g(x,y,\mu,0) \big\}.
\]
Note that $S_0$, which is a $1$-dimensional manifold in this case, appears both as
the set of equilibria of the fast subsystem (where the slow variable $\mu$ is frozen
and considered as a bifurcation parameter) and as the phase space of the slow
subsystem.

According to Izhikevich's classification, elliptic bursting is characterized by the
following pair of fast-subsystem bifurcations: a subcritical Hopf bifurcation (HB),
which triggers the onset of burst, and a saddle node (SN) of limit cycles, which
terminates the burst. This type of bursting is one example where the transition to and from bursting is associated with torus canards (TCs), which emerge in a exponentially-small parameter intervals~\cite{kramer2008new,benes2011elementary,burke2012showcase,ju2018bottom,vo2017generic,roberts2015averaging}.

We review these concepts in a neuronal population Wilson--Cowan-type model~\cite{wilson1972excitatory} which supports
elliptic bursting behaviour; see~\cite{burke2012showcase,izhikevich2000neural}. In the fast scaling, the model reads
\begin{equation}\label{eq:slowSS_meanField}
\begin{aligned}
x'= & -x+N_x(J_{xx}x+J_{xy}y+\mu) ,\\
y'= & \delta(-y+N_y(J_{yx}x+J_{yy}y+\rho)),\\
\mu^\prime= &\varepsilon(k-x-y),
\end{aligned}
\end{equation}
where $x$ and $y$ model a fast excitatory and inhibitory population, respectively,
and where $\mu$ denotes a slow, excitatory, adaptation current. In the model
above, $N_x,N_y$ are the sigmoidal functions,
\begin{equation}
N_x = \frac{\lambda_x}{2}\left[+\mathrm{erf}\left(\frac{g_1x}{\sqrt{2(1+g_1^2\sigma_x^2)}}\right)\right],\quad
N_y = \frac{\lambda_y}{2}\left[1+\mathrm{erf}\left(\frac{g_2y}{\sqrt{2(1+g_2^2\sigma_y^2)}}\right)\right],
\end{equation}
and $J_{xx}$, $J_{xy}$, $J_{yx}$, $J_{yy}$, $\delta$, $\varepsilon$,
$\rho$, $g_1$, $g_2$, $v$, $\lambda_x$, $\lambda_y$, $k$, are control parameters; see also~\cite{touboul2015noise}.

In~\cref{fig:ellipticClassicalScn} (a), we plot a bifurcation diagram of the fast subsystem, onto
which we superimpose an elliptic bursting cycle, whose time profile is shown in the panel (b) of the figure. The critical manifold $S_0$ is a cubic-like curve, and the
fast subsystem bifurcations relevant to elliptic bursting are the right-most HB and
the SN of limit cycles.

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{Ellburst.pdf}}
\caption{Elliptic bursting cycle from the Wilson-Cowan-type system~\cref{eq:slowSS_meanField}. (a) Projection onto the $(z,x)$ plane superimposed onto the bifurcation diagram of the fast subsystem. (b) Time profile of the $x$ variable over a few periods of the bursting cycle. Parameter values are: $J_{xx}=12$, $J_{xy}=-4$, $J_{yx}=13$, $J_{yy}=-9$, $\delta=0.05$, $\varepsilon=0.001$, $\rho=-1.1$, $g_1=0.19$, $g_2=0.4$, $\sigma_x=\sigma_y=1.001$, $\lambda_x=2$, $\lambda_y=8.5$, $k=2.56$.}
\label{fig:ellipticClassicalScn}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{TorusCanardsWC}
\caption{Torus canard dynamics in the Wilson-Cowan type model~\cref{eq:slowSS_meanField}. Left and central panels show torus canard orbits in the $(z,x)$ plane superimposed onto the fast subsystem bifurcation diagram; the latter are zoomed views of the former. Right panels show the associated time series for $x$. Rows (ai), (bi), (ci) and (di) correspond in the same order to a headless TC, a TC with head, a headless mixed-type TC and a mixed-type TC with head. Fixed parameter values are as in~\cref{fig:ellipticClassicalScn}. Panels (ai)-(bi): $\varepsilon=0.0001$, $k=2.536083$ (a) and $k=2.536084$ (b). Panels (ci)-(di): $\varepsilon=0.001$, $k=3.60118$ (c) and $k=3.6011$ (d).}
\label{fig:TorCanWC}
\end{figure}

As parameter $k$ varies, we observe transitions to different regimes. When
$k$ decreases, we observe a transition to tonic
spiking. This transition is organised by TC cycles (see~\cref{fig:TorCanWC}(a)-(b)).
One geometrical property of TCs is that their fast oscillations do not stop
near the SN point of the fast subsystem, but instead continue along the branch of
unstable limit cycles. This feature justifies the term \textit{canard} because the
fast oscillations follow a repelling object of the fast
subsystem~\cite{BCDD1981,DR96,krupa01}. The name TC reflects the fact that this
transition is accompanied by a torus bifurcation in the full system (not shown), that
is, when $\varepsilon\neq 0$ and $\varepsilon\ll 1$ . As other canard solutions, torus canards
come both with head (\cref{fig:TorCanWC}(b)) and without head (\cref{fig:TorCanWC}(a)): the latter winds around an invariant torus of
the full system ($0 < \varepsilon\ll 1$).

In addition to the aforementioned one, the Wilson--Cowan-type system generates
another type of torus canard that seems to be less studied. We refer to it as \textit{mixed-type torus canard} for its link to the notion of canards of mixed type reported
in~\cite{desroches2012canards}. Like the classical TCs, these cycles are related to
the presence of a Neimark Sacker bifurcation in the full system. The name mixed type is used
because these solutions contain an attracting slow segment followed by a repelling
fast segment, hence they differ from classical TCs; mixed-type TCs appear near a
subcritical HB as shown in~\cref{fig:TorCanWC}(c)-(d). This is different from the classical TCs, which appear near a SN of cycles of the
fast subsystem.

In the present paper, we are motivated by studying minimal elliptic bursters that can
display all the solution types mentioned above. Classical TCs have been studied in
various neuronal models~\cite{kramer2008new,benes2011elementary,burke2012showcase}.
Mixed-type TCs, have received less attention, albeit they have been observed and
described in a recent study by Ju, Neiman and Shilnikov~\cite{ju2018bottom}, where
the main focus is on the full-system torus bifurcation.

The canonical form proposed by Izhikevich~\cite{izhikevich2000subcritical} for Bautin bursters, was adapted in~\cite{ju2018bottom} to the framework of elliptic bursting systems. Dynamics relevant to elliptic bursting, in particular torus canards, was therein studied with this adapted canonical form. We will refer to this adapted form as the \textit{canonical elliptic burster}. This model is considered to be a minimal
system for general elliptic bursting dynamics but, as we will show below, it does not
support mixed-type TCs. We first study this model and in particular we find that one boundary of the elliptic bursting regime is degenerate, owing to the presence of a continuum family of
equilibria. At such boundary, the critical manifold $S_0$ has a non-transversal
intersection with the slow nullcline. To overcome this limitation, we introduce a minimal perturbation of the canonical elliptic burster, which restores transversality at the aforementioned boundary. Consequently, the boundary in this new model is organised around mixed-type TCs. And away from this boundary, the new model behaves just like the canonical elliptic burster. Therefore, this new model appears to be a minimal torus canard system capturing transitions to and from the elliptic bursting regime.

In order to analyse all possible dynamics of the canonical elliptic burster as well as our minimal model, we mostly use singular flows and their concatenation, which define singular cycles. Full system's solutions, both of TC and of elliptic bursting type, stay close to such singular cycles for $\varepsilon \neq 0$ small enough. While studying the different singular limits to construct such cycles has been a classical tool in planar canard systems~\cite{krupa01} as well as in systems with mixed-mode oscillations~\cite{brons2006}, it has not been used in bursting systems so far. Indeed, singular cycles perturbing to bursting cycles have not been studied at all, let alone in the TC regime; this is one of the main contributions of the present work.
To do so, we use an additional subsystem obtained in the
singular limit $\varepsilon=0$, the so-called \textit{average slow subsystem}. This subsystem
approximates the slow flow on the manifold of limit cycles of the fast subsystem by averaging out the fast oscillations.  As a result,
the average slow flow retains a slowly-varying observable from the drifting fast
oscillations of the burst, and it allows to trace segments of singular solutions in
an additional $\varepsilon=0$ slow limit. We then can realise all possible singular cycles by suitably concatenating segments from fast, slow and average slow subsystems, and we can understand how to recover singular equivalents of both classical and mixed-type TCs.

The manuscript is organized as follows. In \cref{sec:singular}, we introduce the concept of singular orbits and explain how they help to understand all possible dynamics of the full system, using the classical van der Pol model. Then in \cref{sec:canonicalforms}, we present and analyse both the canonical elliptic burster and our proposed generalized form; in particular, we focus on their singular orbits and relate them to full system's dynamics. Finally, in the conclusion section, we summarise our findings and propose a number of perspectives and questions for future work.

\section{Singular orbits}\label{sec:singular}
We review here the concept of singular
orbit~\cite{benoit1990,brons2006,krupa01}, which will be relevant for the analysis of
elliptic bursters presented below. The main idea behind singular orbits is to
extract information on solutions to the slow-fast dynamical system~\cref{eq:fastTime}, with
$0 < \varepsilon \ll 1$, by constructing orbits for $\varepsilon=0$. The latter are therefore
candidate orbits, that are hopefully descriptive of the perturbed orbits with $\varepsilon
\neq 0$. As we have seen, there exist two non-equivalent singular
limits in an $m$-slow $n$-fast dynamical system, namely the
fast subsystem~\cref{eq:fastSubs} and the slow subsystem~\cref{eq:slowSubs}, respectively. It is therefore natural to construct a singular orbit by concatenating
orbit segments of the fast- and slow-subsystem, as the following definition suggests
\begin{definition}[Singular orbit] \label{def:singOrbit}.
  Let $0=t_0 < t_1 < \cdots < t_l = T$
  be a partition of $[0,T] \subset \mathbb{R}$, for some $T>0$. A singular orbit
  starting at $p$ is a continuous curve $\gamma \colon [0,T] \to \mathbb{R}^{n+m}$
  \[
    \gamma(t) =
    \begin{cases}
      \gamma_1(t)  & t \in [t_0,t_1], \\
      \gamma_2(t)  & t \in [t_1,t_2], \\
      \ldots       &                 \\
      \gamma_{m}(t)  & t \in [t_{m-1},t_m]
    \end{cases}
    \qquad
    \gamma_{1}(t) = \varphi_1^t(p),
    \qquad
    \gamma_{i}(t) = \varphi_i^t(\gamma_{i-1}(t_{i-1})),
    \qquad
    i = 2,\ldots,l,
  \]
  where, for each $i=1,\ldots,l$, $(t,p) \mapsto \varphi_i^t(p)$ is the flow associated to either
  \cref{eq:fastSubs} or \cref{eq:slowSubs}.
\end{definition}

In the analysis below, we will also be interested in periodic solutions to the full
problem, and hence, in closed singular orbits. We exemplify the construction of open
and closed singular orbits using the classical Van der Pol system~\cite{BCDD1981,DR96,krupa01}, known to be a prototypical planar canard system
\begin{equation}\label{eq:vandp}
  x' = x - \frac{x^3}{3} - \mu, \qquad \mu' = \varepsilon(x-\mu).
\end{equation}
In \cref{fig:vandP} we show examples of singular trajectories. The top row sketches
slow and fast flows associated to~\cref{eq:vandp}, using single and double arrows,
respectively. The sketches are made for different values of $k$, for which the slow nullcline
intersects the critical manifold to the left of $p_0$ (a1), at $p_0$ (b1), and to the
right of $p_0$ (c1).

When the intersection is to the right of $p_0$, an equilibrium $p$ of the full
problem exists to the left stable branch of $S_0$. From panels (a2)-(a4) we see that
\emph{singular orbits are not unique}: concatenating slow (green) and fast (red) segments
and respecting the orientation on curve patches, one can construct infinitely many
singular orbits starting at
$p_0$, and terminating at $p$. In (a2) we see that an initial slow segment starting
from $p_0$ can be concatenated to any of the purple
segments, giving rise to orbits that jump to the left attracting branch of $S_0$, and
terminate at $p$. Of all such orbits, we highlight one that leaves the repelling
branch of $S_0$ at the origin (red segment). In (a3) we present a singular orbit that
starts at $p_0$ but then jumps at the origin to the right attracting branch of $S_0$,
highlighting in purple other potential concatenations. In (a4) we present a singular
orbit that starts at $p_0$, from where it jumps to the right attracting branch of
$S_0$ at $p_1$; in this case, the singular orbit must follow $S_0$ up to the fold at
$p_2$, from where it jumps towards $p_3$: it is impossible to concatenate to a fast
segment before reaching the fold at $p_2$, while respecting \cref{def:singOrbit}.

When the intersection of the $w$-nullcline is exactly at $p_0$, one can construct
singular cycles (panels (b1)--(b4)). Singular cycles may start at $p_0$, follow the
repelling branch and then jump to the left (b2) or to the right (b3) attracting
branch of $S_0$. A singular cycle that jumps from the starting point $p_0$ to the
right attracting branch of $S_0$ can also be uniquely constructed. When $\varepsilon >0$,
the singular solutions give rise to headless canard cycles (b2), canard cycles with
head (b3), and relaxation a unique relaxation oscillation (b4).

Finally, when the intersection of the $w$-nullcline is to the right of $p_0$, one can
also construct a unique singular cycle corresponding to a relaxation oscillation.

\begin{figure}
  \includegraphics[width=\textwidth]{SingularVDP}
  \caption{Examples of singular orbits in the Van der Pol system~\cref{eq:vandp} for
    different values of $k$ for which the slow nullcline intersects the critical
    manifold to the left of $p_0$ (a1--a4), at $p_0$ (b1--b4), and to the
    right of $p_0$ (c1--c2). (a1), (b1), (c1): sketches of slow (single arrow) and
    fast (double arrows) flows associated to~\cref{eq:vandp}, which give indication
    of how one can construct singular orbits. (a2)--(a4): Singular orbits starting at
    $p_0$ can follow the repelling branch of the critical manifold $S_0$ and then
    jump to the left (a2) or to the right (a3) attracting branch; we indicate slow
    segments in green, and fast segments in red; singular orbits are non unique, and
    we sketch in purple alternative fast segments; in (a4) we show an orbit starting
    at $p_0$ and jumping directly to the right attracting branch $S_0$; this singular
    orbit must follow $S_0$ up to the fold at $p_2$, from where it jumps towards
    $p_3$: it is impossible to concatenate to a fast segment before reaching the fold
    at $p_2$, while respecting \cref{def:singOrbit}, and indeed no purple segment is
    indicated. (b2)--(b4): Singular cycles can be constructed when the slow nullcline
    intersect $S_0$ at $p_0$, and similar considerations to (a2)--(a4) apply here;
    when $\varepsilon >0$, these singular orbits perturb to headless canard cycles (a2),
    canard cycles with head (a3), and relaxation oscillations, respectively. (c2):
    when the $w$-nullcline intersects $S_0$ to the right of $p_0$ one can construct a
    unique singular cycle corresponding to a relaxation oscillation.
}
  \label{fig:vandP}
\end{figure}

\FloatBarrier

\section{Canonical forms for elliptic bursting}\label{sec:canonicalforms}
We now proceed to study the two canonical systems for elliptic bursting of interest
to us. The first one is the canonical elliptic burster in~\cite{ju2018bottom}, adapted  from~\cite[Thm 2.1]{izhikevich2000subcritical},
which we repeat here for convenience
\begin{equation}\label{eq:normalFormClassicalRep}
 \begin{aligned}
   x'   &=  -y+x\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2  \big ),\\
   y'   &=   x+y\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2 \big ),\\
 \mu' & = \varepsilon (k-x^2-y^2).
 \end{aligned}
 \end{equation}
 A polar coordinate transformation $(x,y) = (r \cos \theta, r \sin \theta)$ leads to
 the system
\begin{equation}\label{eq:canonicalStandard}
\begin{aligned}
r' &=  r(\mu+2r^2-r^4),\\
\theta' &= 1,\\
\mu'  &=  \varepsilon (k-r^2).
\end{aligned}
\end{equation}
We will consider $3$ subsystems for the canonical elliptic burster: in addition to the fast and
slow subsystems of \cref{eq:normalFormClassicalRep}, written as
\begin{equation}\label{eq:normalFormClassicalRepFast}
 \begin{aligned}
   x'   &=  -y+x\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2  \big ),\\
   y'   &=   x+y\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2 \big ),\\
 \mu' & = 0,
 \end{aligned}
 \end{equation}
 and
\begin{equation}\label{eq:normalFormClassicalRepSlow}
 \begin{aligned}
   0   &=  -y+x\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2  \big ),\\
   0   &=   x+y\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2 \big ),\\
   \dot \mu  &=  (k-r^2),
 \end{aligned}
 \end{equation}
 respectively, we will consider the averaged slow subsystem~\cite{butera1997transient,roberts2015averaging}, which is
 obtained by averaging the slow equation $\mu$ in the slow-time parametrisation of
 \cref{eq:normalFormClassicalRep} over one period of the fast limit
 cycles of the fast subsystem~\cref{eq:normalFormClassicalRepFast}, and approximating
 $\mu$ by its average $\langle \mu \rangle$. The fast-subsystem limit
 cycles, which we denote $(x_p(t;\mu),y_p(t;\mu))$, have period $2\pi$ for every
 $\mu>0$, hence we obtain
 \begin{equation}
   \dot{\langle \mu \rangle} = \frac{1}{2\pi} \int_0^{2\pi} \big( k - x^2_p(t;\langle
   \mu \rangle) - y^2_p(t;\langle \mu \rangle) \big)\, dt.
 \end{equation}

The fast, slow, and averaged slow subsystems discussed above have a compact
representation in polar coordinates, given by
\begin{align}
  r' &=  r(\mu+2r^2-r^4), && \theta' = 1, && \mu'  =  0 , \label{eq:BautinFS}\\
  0  &=  r,               &&              && \dot \mu  = k-r^2, \label{eq:BautinSS}\\
  0  &=  \mu +2r^2 - r^4, &&              && \dot \mu  = k-r^2, \label{eq:BautinASS}
\end{align}
respectively, where we have replaced $\langle \mu \rangle$ by $\mu$ in the last equation.
Note that there is no $\theta$ dynamics in none of the slow and average slow subsystems.

As stated in the introduction, in addition to the canonical elliptic burster, we
consider an extended model, with a perturbative term in the slow equation for $\mu$
\begin{equation}\label{eq:leidenator}
 \begin{aligned}
   x'   &=  -y+x\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2  \big ),\\
   y'   &=   x+y\big ( \mu+2(x^2+y^2)-(x^2+y^2)^2 \big ),\\
   \mu'  &=  (k-x^2-y^2 - \alpha \mu),
 \end{aligned}
 \end{equation}
 where $\alpha >0$. We will henceforth refer to the model above, or to its polar
 representation
\begin{equation}\label{eq:LeidenStandard}
\begin{aligned}
r' &=  r(\mu+2r^2-r^4),\\
\theta' &= 1,\\
\mu'  &=  \varepsilon (k-r^2 - \alpha \mu),
\end{aligned}
\end{equation}
 as the Leidenator model. A stochastic version of this model was simulated in a study on
 epilepsy~\cite[Supplementary material]{hebbink2017}. Reasoning as for the canonical
 elliptic burster, we arrive at the following slow, fast, and averaged slow
 subsystems for the Leidenator
\begin{align}
  r' &=  r(\mu+2r^2-r^4), && \theta' = 1, && \mu'  =  0 , \\
  0  &=  r,               &&              && \dot \mu  = k-r^2 - \alpha \mu, \\
  0  &=  \mu +2r^2 - r^4, &&              && \dot \mu  = k-r^2 - \alpha \mu,
\end{align}
respectively. As expected, the canonical and the Leidenator bursters differ only in
the slow subsystems.

In the following, we will analyse for each $k$ interval both the standard canonical
form and the Leidenator model in terms of singular and perturbed (non-singular)
dynamics. We will focus on selected singular orbits in each case, which are the most
relevant for the bursting.

\begin{figure}[!b]
\centerline{\includegraphics[scale=0.29]{preLeidenator.png}}
\caption{Sketches of the singular orbits generated by the canonical burster as $k$
  varies. Orbit segments of the slow/average slow subsystems
  \cref{eq:BautinSS,eq:BautinASS} are indicated in green, and segments of the fast
  subsystem \cref{eq:BautinFS} are indicated in red. We denote by SN and HB
  bifurcations of the fast subsystem. Singular orbits are constructed by
  concatenating red and green segments, with concatenation points indicated by dots.
  Purple fast segments indicate alternative fast orbit segments. Solid (dashed)
  curves indicate branches of stable (unstable) attractors of the fast subsystem.
}
\label{fig:preLeidenator}
\end{figure}

\subsection{Canonical elliptic burster}
We now consider the fast subsystem \cref{eq:BautinFS} of the canonical elliptic
burster. The critical manifold is given by the horizontal line $r=0$. The fast
subsystem has 2 main bifurcations, a subcritical Hopf, HB, and a saddle-node of
cycles, SN, as shown in \cref{fig:preLeidenator}. The fast subsystem is therefore
bistable, and the branch of periodic orbits emanating from HB are circles of radius
$r = \sqrt{\mu}$, and lie in the $(r,\mu)$ plane on the curve $0 = \mu + 2r^2 -r^4$.

We will study the dynamics of the full system when $k$ is varied. As we shall see
below, the asymptotic dynamics of the full system is either a stable limit cycle, or
a bursting cycle, and the transition between these regimes occurs around a value $k =
\sqrt{r_{SN}}$, where $r_{SN}$ is the fold point of $0=\mu + 2r^2 - r^4$.
Henceforth, we will pose $k_{SN} =\sqrt{r_{SN}}$. We will also show below that the bursting regime
terminates around $k=0$, where the dynamics is degenerate.

\subsubsection*{Case 1: $\boldsymbol{k>k_{SN}}$}
\paragraph{Singular dynamics} The slow average subsystem has a unique stable
equilibrium at $p_3$ as shown in~\cref{fig:preLeidenator}(a). All singular orbits terminate at $p_3$ in this region of parameter space. We show
the construction of a typical orbit, which starts on the critical manifold at $p_0$,
follows the critical manifold past HB, up to an arbitrary large value $p_1$. From
$p_1$ the orbit concatenates with a fast segment to $p_2$; the final segment from
$p_2$ to $p_3$ is on the slow averaged flow. There is a one-parameter family of
orbits connecting $p_0$ to $p_3$; this family is parametrised by $p_1$.

\paragraph{Nonsingular dynamics}
The full system ($\varepsilon \neq 0$) selects a unique orbit starting at $p_0$; this orbit
approaches an $\varepsilon$-perturbation of $p_3$, that is, a tonic firing solution of the
full system. This orbit lifts off from an $\varepsilon$-neighbourhood of the critical manifold, at a
point uniquely determined using the so-called \textit{way-in-way-out
function}~\cite{diener1991} once $k$, $\varepsilon$ and $p_0$ are fixed.

\subsubsection*{Case 2: $\boldsymbol{k=k_{\text{SN}}}$}

\paragraph{Singular dynamics}
In this case, the $\mu$-nullcline crosses the fold at SN, hence the slow average system does not have an equilibrium, but a turning
point at SN and, as for the van der Pol system, this generates singular canard cycles,
as presented in the three panels \cref{fig:preLeidenator}(b1)--(b3).

In the first scenario, which is illustrated in~\cref{fig:preLeidenator}(b1), we
construct a one-parameter family of singular orbits parametrized by $p_1$,
consisting of one average slow segment, which goes through SN, and
one fast segment. SN is a canard point for the $(r,\mu)$ averaged slow flow
system~\cref{eq:BautinASS}. Each member of this one-parameter family is uniquely defined based
on the location of $p_1$, which can vary between SN and HB (excluded), thereby generating
singular canard orbits in~\cref{eq:BautinASS}. We call these orbits singular TC
orbits because they correspond to TC orbits in the full system.

In the second scenario, in~\cref{fig:preLeidenator}(b2), $p_1$ occurs at HB, or at an
arbitrary point to the right of HB on the critical manifold. These orbits are the
so-called \textit{singular maximal TCs}, because their repelling segment on
the branch of periodic orbits of the fast subsystem has maximal length.

In the third scenario we observe a
two-parameter family of singular TCs with head (with 2 fast segments) and parametrized by $p_1$ and $p_3$;
see \cref{fig:preLeidenator}(b3). The location of $p_1$, which is arbitrarily chosen
between SN and HB, determines the canard's head size. As in the second scenario,
$p_3$ can be located anywhere on the unstable critical manifold.

\paragraph{Nonsingular dynamics} When $\varepsilon \neq 0$, we expect to find TC cycles
in the full system~\cref{eq:normalFormClassicalRep}; these orbits will be close to
singular orbits from the three scenarios mentioned above, and parametrised by $k$ in
an exponentially small region of parameter space. This one-parameter family of TC
cycles starts with headless TCs realised from~\cref{fig:preLeidenator}(b1). Out of
the singular maximum TC cycles of~\cref{fig:preLeidenator}(b2), we expect that only
one is realised, namely the one for which $p_1$ is at HB. This is best understood by
looking at the third scenario~\cref{fig:preLeidenator}(b3): computing the
way-in-way-out function, we see that $p_2$ and $p_3$ must be equidistant from HB.
This also implies that, when $\varepsilon \neq 0$, the TC cycles with realised in the full
system are a one-parameter family, depending solely $p_2$, as opposed to their
singular counterpart which appear in a 2-parameter family parametrised by $p_2$ and
$p_3$.

\subsubsection*{Case 3: $\boldsymbol{0<k<k_{\text{SN}}}$}

\paragraph{Singular dynamics}
In this case the equilibrium of the averaged slow system~\cref{eq:BautinASS} is on
the unstable branch of limit cycles, therefore singular orbits do not display canard
segments, as they connect SN to a point $p_1$ on the critical manifold. One can
construct a one-parameter family of singular cycles, parametrised by the value of
$p_2$ at which the singular orbit is concatenated to a fast segment, as shown
in~\cref{fig:preLeidenator}(c). We call these orbits singular elliptic bursting
cycles.

\paragraph{Nonsingular dynamics}
For each $k,\varepsilon>0$, we expect that only elliptic bursting orbit persists, selected by the way-in-out
function as above.

\subsubsection*{Case 4: $\boldsymbol{k=0}$}

\paragraph{Singular dynamics}
In this case, the slow subsystem has a continuum of trivial equilibria on the
critical manifold, as illustrated in~\cref{fig:preLeidenator}(d). All singular orbits
starting at a point $p_0$ on the stable slow average branch terminate at $p_1$, the
projection of SN on $S_0$. One can construct other singular orbits, terminating at
any of the equilibria on $S_0$, by varying $p_0$.

\paragraph{Non-singular dynamics}
$S_0$ is now also a line of equilibria for the full
system~\cref{eq:normalFormClassicalRepFast}. The non-singular behaviour is
qualitatively the same of the singular one: the system comes to a rest on $S_0$,
after a jump from the branch of stable limit cycles.

\subsubsection*{Case 5: $\boldsymbol{k<0}$}
\paragraph{Singular and nonsingular dynamics}
There is no equilibrium of the slow subsystems. Both the singular and the nonsingular
dynamics display a drift in $\mu$ towards $-\infty$, while $r(t)$ tends monotonically
to $0$.

\subsection{The Leidenator model}
The canonical bursting model captures well the transition from tonic firing to
elliptic bursting via TCs, observed in most known elliptic bursters, by decreasing
$k$. However, as $k$ decreases further, the end of the elliptic bursting regime
occurs via a continuum of equilibria, which is non generic, because the fast nullcine
(the critical manifold) intersects non transversally the $\mu$-nullcline (the two
coincide). In addition, in numerical simulations of the Wilson-Cowan-type model~\cref{eq:slowSS_meanField}, we have seen that the end of the bursting
regime occurs via mixed-type TC cycles, like the ones presented
in~\cref{fig:TorCanWC}(ci)-(di). The canonical elliptic burster can not capture
solutions of this type.

We introduce the Leidenator model in order to overcome these problems. In this model, the slow nullcline is the parabola $\mu = (k - r^2) / \alpha$ and it intersects transversally the critical manifold for all $k$ values. For sufficiently small $\alpha>0$, one therefore obtains an isolated equilibrium of the slow flow; see~\cref{fig:preLeidenator}. In addition, for $k>0$, this equilibrium corresponds to a so-called \textit{buffer point}~\cite{diener1991}, which affects the non-singular dynamics as will be explained below. Note that one recovers the canonical elliptic burster from the Leidenator for $\alpha=0$.

\subsubsection*{Cases 1-3: $\boldsymbol{k>0}$}
As seen in \cref{fig:captureLeidenator}(a)-(c), the Leidenator behaves like the canonical elliptic bursting model for $k>0$. The only difference is the existence of $p_b$, which determines an upper bound for the points on $S_0$ where the right-most fast segment can be placed. For $\varepsilon \neq 0$, the buffer point at $p_b$ marks the maximal delay that any full-system trajectory can have.
\begin{figure}
\centerline{\includegraphics[scale=0.27]{leidenator1.png}}
\caption{Sketches of the singular orbits generated by the Leidenator model. Graphical conventions as in \cref{fig:preLeidenator}. A buffer point $p_b$ arises and is indicated by a blue dot. Panels (d1) and (d2) show new configurations, which give rise to singular mixed-type TC with head and without head, respectively.}
\label{fig:captureLeidenator}
\end{figure}
\subsubsection*{Case 4: $\boldsymbol{k=0}$}
\begin{figure}
\centerline{\includegraphics[scale=0.45]{torcan.png}}
\caption{TC cycles computed in the Leidenator model. From left to right: headless TC, maximal TC and TC with head. Top panels display the time profile of the $x$-variable ($x=r\cos\theta$). Bottom panels show the projection of the solution onto the $(\mu,r)$ phase plane together with the fast subsystem bifurcation diagram (the dots correspond to the HB and the SN points). Solutions of these types are also supported by the canonical elliptic burster. Parameter values: $\varepsilon=0.001$, $\alpha=0.2$. The value of $k$ is different in each panel with numerical values within $10^{-12}$ of  $0.79984990182$.}
\label{fig:torcanleiden}
\end{figure}
\paragraph{Singular dynamics}
This case differs substantially from the canonical elliptic burster. Indeed, $p_b$ coincides with HB and there is no equilibrium of the average slow flow (point $p$ from \cref{fig:preLeidenator}(c) is no longer present). This provides a new way to concatenate singular segments, namely slow and average slow segments, so as to form new types of singular canard cycles. We refer to these cycles as singular mixed-type TC cycles. Similar to all singular canard cycles, we come both with and without head, which are shown in~\cref{fig:captureLeidenator}(d1) and (d2), respectively. At the transition between these two subfamilies, there is a unique singular maximal mixed-type TC cycle.

\begin{figure}
\centerline{\includegraphics[scale=0.45]{mixedtypetorcan.png}}
\caption{Mixed-type TC cycles computed in the Leidenator model. From left to right: headless mixed-type TC, maximal mixed-type TC and mixed-type TC with head. Top panels display the time profile of the $x$-variable ($x=r\cos\theta$). Bottom panels show the projection of the solution onto the $(\mu,r)$ phase plane together with the fast subsystem bifurcation diagram (the dots correspond to the HB and the SN points). These solutions are not supported by the canonical elliptic burster. Parameter values: $\varepsilon=0.05$, $\alpha=0.2$. The value of $k$ is different in each panel with numerical values within $10^{-14}$ of  $0.0015128438002$.}
\label{fig:mixtyptorcanleiden}
\end{figure}

\paragraph{Non-singular dynamics}
When $\varepsilon \neq 0$, we expect to find mixed-type TC cycles
in the full system~\cref{eq:leidenator}; these orbits will be parametrised by $k$ in
an exponentially small region of parameter space and will span canard cycles with and without head.

\subsubsection*{Case 5: $\boldsymbol{k<0}$}
This case is identical to Case 5 in the canonical elliptic burster.

\subsection{Numerical results}
We have used numerical continuation with the software package \textsc{auto}~\cite{auto} in order to compute canard cycles in the $(r,\mu)$ system and therefore obtain TC cycles in the Leidenator model via the polar coordinate transform. The results are included in \cref{fig:torcanleiden}, which shows TC cycles corresponding to Case 2, as well as in~\cref{fig:mixtyptorcanleiden}, which shows mixed-type TC cycles corresponding to Case 4.

\section{Conclusion}
\label{sec:conclusion}
In this paper, we analysed the dynamical repertoire of elliptic bursting systems by using two minimal systems, namely the canonical elliptic burster~\cref{eq:normalFormClassicalRep} and the Leidenator~\cref{eq:leidenator}. The presented results focuses on singular orbits in order to show that the canonical elliptic burster cannot produce all dynamical patterns that an elliptic burster may display. More precisely, it cannot possess mixed-type TCs due to the presence of a continuum of equilibria which correspond to one boundary of the elliptic bursting regime where mixed-type TC should appear.  For this reason, we extended the canonical elliptic burster into the Leidenator, which manage to overcome this limitation and possess generic transitions involving torus canards (classical and mixed type) at both ends of the bursting regime.

Our approach is based on the singular limits of slow-fast systems of elliptic bursting type, and it employs three limiting systems obtained for $\varepsilon = 0$: the fast, slow and average slow subsystems. We define a singular orbit as a compatible concatenation of orbit segments solutions to these three subsystems. The idea behind the singular orbits is to extract information about the full system with $0<\varepsilon \ll 1$ by constructing piecewise-defined orbits valid for $\varepsilon = 0$. The orbits of the full system remain close to the singular orbits. Therefore, the singular orbits contain information about the full system within a simplified setting.

In parameter space, the transition to and from elliptic bursting is organized by TCs, which appear in exponentially-small parameter intervals. Two families of TCs are then of interest: classical and mixed-type. In both families, canard solutions are divided into headless canards, maximal canard and canards with head. We relate such TC solutions, as well as the bursting and tonic spiking solutions, to appropriately constructed singular orbits. As we do so, the fact that the canonical elliptic burster~\cref{eq:normalFormClassicalRep} cannot produce mixed-type TCs becomes clear and show the limitation of this model, as mixed-type TCs are observed in elliptic bursters, e.g. in the Wilson-Cowan-type model given by~\cref{eq:slowSS_meanField}. This is due to the presence of a continuum of equilibria resulting from a nontransversal intersection between the (slow) $\mu$ nullcline and the critical manifold for $k=0$, as shown in~\cref{fig:preLeidenator}(d1). We then propose the Leidenator system, in which we introduce a perturbative term in the slow flow that minimally suffices to overcome this problem. This term makes the $\mu$ nullcline (a parabola) keep a transversal intersection with the critical manifold for any $k$ value, as shown in~\cref{fig:captureLeidenator}. Consequently, this guarantees that, for proper parameter choices, the full system has a unique equilibrium when $k=0$. At the level of singular systems, it enables the construction of singular mixed-type TCs when $k=0$, which is not possible in the case of the canonical elliptic burster.

Using various $\varepsilon=0$ systems in order to construct singular cycles, and hence infer the possible dynamics of a slow-fast system, is a standard tool that has been regularly used in problems with canard cycles~\cite{krupa01} or mixed-mode oscillations~\cite{brons2006}. However, it has hardly been used in the context of bursting systems, where a third singular flow must be considered. The present paper shows the advantage of this approach in order to derive valuable information in a given family of bursting systems, namely elliptic bursters. We anticipate that it could give useful in order type of bursters, starting by the other two classes introduced by Rinzel~\cite{rinzel1987formal} which are square-wave and parabolic bursters. This is an interesting question for future work. Also, the question of rigorous existence of mixed-type TCs is very natural and we will consider it in a follow-up study. Finally,
the occurrence of mixed-type TCs in biophysical neuron model, and their physiological relevance, is also of interest, possibly in link with resonator neurons that display subthreshold oscillations, as well as the possibility to find such structures in a spatially-extended context.

\bibliographystyle{siamplain}


\section{Introduction}

The stochastic gravitational wave background is
of great interest of cosmologists and astrophysicists
since the detection of its power spectrum reveals to us
high-energy events and their unknown physics in the early universe.
The power spectrum for $f\lesssim 10^{-15}\,{\rm Hz}$ is constrained
by experiments of the cosmic microwave background (CMB),
whereas that in the high-frequency range is constrained by
experiments of pulsar-timing arrays and space/ground-based interferometers.
Recently, the NANOGrav collaboration reported that they found a strong preference
for a stochastic common-spectrum process in their
12.5-yr data set~\cite{Arzoumanian:2020vkk}
although definitive evidence of quadrupolar interpulsar
correlations~\cite{Hellings:1983fr} is lacked.
If the common-spectrum process is due to supermassive black hole binaries,
it would be the first evidence of their formation and coalescence.
Other possible explanations for the NANOGrav result have been proposed by recent studies~\cite{Ellis:2020ena, Blasi:2020mfx, Buchmuller:2020lbh, Samanta:2020cdk, Nakai:2020oit, Addazi:2020zcj, Ratzinger:2020koh, Neronov:2020qrl, Abe:2020sqb,Vaskonen:2020lbd, DeLuca:2020agl, Kohri:2020qqd, Domenech:2020ers, Bhattacharya:2020lhc, Namba:2020kij, Kitajima:2020rpm, Bian:2020bps, Vagnozzi:2020gtf, Li:2020cjj}
from the view point of
cosmic strings~\cite{Ellis:2020ena, Blasi:2020mfx, Buchmuller:2020lbh, Samanta:2020cdk},
dark phase transition~\cite{Nakai:2020oit, Addazi:2020zcj, Ratzinger:2020koh},
QCD phase transition~\cite{Neronov:2020qrl, Abe:2020sqb}, and
second-order gravitational waves from
primordial black holes~\cite{Vaskonen:2020lbd, DeLuca:2020agl, Kohri:2020qqd, Domenech:2020ers, Bhattacharya:2020lhc}
or resonant amplification of other fields~\cite{Namba:2020kij,Kitajima:2020rpm}.

Among a variety of possibilities, in this paper we discuss whether it is possible
to generate gravitational waves of inflationary origin~\cite{Guzzetti:2016mkm}
that have amplitudes large enough to be
detected in the nanohertz range with pulsar timing arrays while evading
other observational constraints.
Taking the tensor amplitude and spectral index simply as parameters,
one can derive the region of the parameter space consistent with
the NANOGrav result at $f\sim 1/$yr and the Planck limit at $f\lesssim 10^{-15}\,$Hz.
However, if one assumes that the gravitational wave energy density $\Omega_{\rm GW}(f)$
has a single power-law form with $n_t\sim 1$ up to the high frequency end of the spectrum,
Big Bang nucleosynthesis (BBN) then excludes the entire parameter space~\cite{Vagnozzi:2020gtf}.
To explore the inflationary explanation for the NANOGrav result,
we therefore have to resolve the following issues.
First of all, it is difficult to produce such an extremely blue gravitational wave
spectrum directly from inflation.
A blue inflationary gravitational wave spectrum implies the violation of the
null energy condition (NEC), which usually results in some kinds of instabilities.
Even if stable violation of the NEC were possible,
generating an extremely blue spectrum with $n_t\sim 1$ would be even more difficult
because the Hubble expansion rate that determines the tensor
amplitude is supposed to be nearly constant during inflation.
Furthermore, the spectrum must be less steep for $f\gtrsim 1/$yr
to evade the BBN constraint~\cite{Kuroyanagi:2014nba}.

At the level of the background solution, one can
obtain NEC-violating inflation (super-inflation, $\dot H>0$) with a phantom scalar
field~\cite{Piao:2003ty,Piao:2004tq,Piao:2004jg,Piao:2006jz,Piao:2007sv}
or a k-essence field~\cite{Baldi:2005gk}.
Unfortunately, at the level of linear perturbations
all of these earlier models are unstable. It was, however, pointed out that
stable violation of the NEC is possible in a consistent effective field
theory~\cite{Creminelli:2006xe}. Indeed, galileon and Horndeski
theories~\cite{Horndeski:1974wa,Deffayet:2011gz,Kobayashi:2011nu,Kobayashi:2019hrl} admit
stable NEC-violating inflation~\cite{Deffayet:2010qz,Kobayashi:2010cm}.
(See also~\cite{Anisimov:2005ne} for a higher-derivative scalar-field theory
admitting a similar super-accelerating phase.)
In the same context, a more radical scenario has been proposed
in which the universe starts expanding from Minkowski
(galilean genesis)~\cite{Creminelli:2010ba}.
Variants of the galilean genesis scenario are found in~\cite{Liu:2011ns,Liu:2012ww,Creminelli:2012my,Hinterbichler:2012fr,Hinterbichler:2012yn,Liu:2013xt,Elder:2013gya,Pirtskhalava:2014esa,Nishi:2015pta,Kobayashi:2015gga,Cai:2016gjd,Nishi:2016ljg,Volkova:2019jlj,Ilyas:2020zcb}.
Some of these introduce an NEC-violating phase connected smoothly to
a subsequent quasi-de Sitter inflationary phase~\cite{Pirtskhalava:2014esa,Kobayashi:2015gga}.
See~\cite{Rubakov:2014jja} for a review on NEC violation and
its cosmological consequences.
An NEC-violating inflationary solution has also been considered
in different contexts such as loop quantum gravity~\cite{Copeland:2007qt,Copeland:2008kz}.

In the above NEC-violating examples of the early universe scenarios,
gravitational waves originated from
quantum fluctuations of the metric typically have a blue spectrum.
This is because the amplitude of the metric perturbations is basically given by
$h_{ij}\sim H/M_{\rm Pl}$, where $H$ is the Hubble parameter at the horizon crossing time
(which is increasing when the NEC is violated)
and $M_{\rm Pl}$ is the Planck mass.
Note, however, that its increasing rate, $\dot H/H^2$, is small
in quasi-de Sitter inflation and consequently the tensor spectral index
$n_t$ is much smaller than one. (See, e.g.,~\cite{Mishima:2019vlh,Capurri:2020qgz}
for recent analyses of a blue tensor tilt from slow-roll inflation.)
Therefore, the expansion of the early universe must be accelerating
and significantly away from de Sitter
to produce gravitational waves with $n_t\sim 1$.
(There is a possibility that the primordial tensor spectrum
is blue-tilted due to the time variation of the tensor propagation speed~\cite{Cai:2015yza,Cai:2016ldn},
but its effect can be absorbed entirely into the time variation of
the Hubble parameter in the Einstein frame~\cite{Creminelli:2014wna}.
A non-standard symmetry-breaking
pattern~\cite{Gruzinov:2004ty,Endlich:2012pz,Cannone:2014uqa,Bartolo:2015qvr,Ricciardone:2016lym,Fujita:2018ehq}
and a non-Bunch-Davis initial state~\cite{Ashoorioon:2014nta} can also yield
a blue tensor tilt.)

The purpose of the present paper is to propose an early universe scenario
that can generate large primordial gravitational waves
in the nanohertz range and at the same time is consistent with other observations.
Our scenario is composed of NEC-violating super-inflation
and a subsequent (relatively short)
phase of NEC-preserving slow-roll inflation.
We construct an explicit model within the cubic Horndeski theory.
Our scenario is closely related to, but different in many aspects from,
those in~\cite{Pirtskhalava:2014esa,Kobayashi:2015gga}.

This paper is organized as follows.
In the next section, we give a brief sketch of the scenario.
We then present in more detail the explicit Lagrangian for our model
and study the background cosmological evolution in Sec.~\ref{sec3}.
The gravitational wave spectrum generated in our model is presented in Sec.~\ref{sec4}.
The mechanism for producing scalar perturbations is discussed in Sec.~\ref{sec5}.
Section~\ref{sec6} is devoted to conclusions.

\section{Sketch of the Idea}\label{sec2}

Let us begin with introducing our basic idea
for generating an extremely blue tensor spectrum
in the early universe.
Our scenario is composed of the following three phases:
an NEC-violating phase, subsequent relatively short inflation, and (standard) reheating.
Each phase is characterized by the different dynamics of a scalar field in the potential,
as shown schematically in Fig.~\ref{fig:potential.pdf}.

  \begin{figure}[tb]
    \begin{center}
            \includegraphics[keepaspectratio=true,height=55mm]{potential.pdf}
       \caption{The shape of the potential and the dynamics of the scalar field
       in our scenario.
  	}
       \label{fig:potential.pdf}
  	\end{center}
  \end{figure}

The scale factor and the Hubble parameter during the NEC-violating phase we consider
are given, respectively, by
\begin{align}
    a\propto (-t)^{-p},\quad H=\frac{p}{-t},\label{eq:Hubble-NECV}
\end{align}
where $t\,(<0)$ is the proper time and $p\,(>0)$ is a constant.
Clearly, the NEC is violated: $\dot H= p/(-t)^2>0$.
Note that in this phase the universe undergoes an
accelerated expansion, $\ddot a/a=p(1+p)/(-t)^2>0$, so that
the horizon and flatness problems are resolved in the same way
as in standard quasi-de Sitter inflation.
The amplitude of tensor modes is basically determined by
the Hubble parameter at horizon crossing.
In this NEC-violating phase where $H$ is increasing with time as in Eq.~\eqref{eq:Hubble-NECV},
the tensor spectral index is given by
\begin{align}
    n_t=\frac{2}{p+1}, \label{eq:tensor-index}
\end{align}
and hence $0<n_t<2$ for $p>0$.

The cosmic expansion with $\dot H>0$ is caused by
a scalar field having (effectively) negative kinetic energy
and climbing up the potential.
It is likely that the curvature perturbation exhibits ghost/gradient
instabilities around such a background, which indeed occurs in
$P(\phi,(\partial\phi)^2)$ theory. However, by adding a galileon-type interaction term
$\sim (\partial\phi)^2\Box\phi$ to the Lagrangian,
one can stabilize such an NEC-violating background,
as is demonstrated in Refs.~\cite{Creminelli:2010ba,Deffayet:2010qz,Kobayashi:2010cm}.
We also use this idea to stabilize the NEC-violating phase.

The first phase described above ends when the scalar field goes over
the top of the potential. The Hubble parameter takes the maximum value $H_*$
at that moment.
The scalar field then rolls down the potential and
quasi-de Sitter inflation occurs. This inflationary phase is
rather standard except that its duration is relatively short.
The tensor modes generated during this phase have a nearly flat spectrum.
Finally, the scalar field starts to oscillate around the bottom of the potential
and thereby the reheating process proceeds in a standard way.

In the next section we present a concrete model realizing the entire
cosmological evolution of this scenario.

\section{Model}\label{sec3}

In this section, we seek for the Lagrangian that admits
the aforementioned early universe scenario within
the following subclass of the Horndeski theory:
\begin{align}
    {\cal L}=G_2(\phi,X)-G_3(\phi,X)\Box\phi+\frac{M_{\rm Pl}^2}{2}R,\label{cubic-Horndeski}
\end{align}
where $X=-g^{\mu\nu}\partial_\mu\phi\partial_\nu\phi/2$.
The general background and perturbation equations in this theory
are summarized in the appendix.

\subsection{The Lagrangian for the NEC-violating phase}\label{subsec3a}

Before investigating the full Lagrangian
that is capable of describing the entire cosmological evolution
sketched in the previous section,
let us focus on the Lagrangian just for the NEC-violating phase.
The Lagrangian we consider for the moment is given by
\begin{align}
{\cal L}=\frac{M_{\rm Pl}^2}{2}R +\beta X-V(\phi)-\frac{\alpha M_{\rm Pl}}{U(\phi)}X\Box\phi,
\label{Lag:NEC-V}
\end{align}
where
\begin{align}
   V(\phi)=U(\phi)=U_0e^{2\phi/M_{\rm Pl}},
   \quad U_0={\rm const}>0,
   \label{Pot:U=V}
\end{align}
and $\alpha$ and $\beta$ are dimensionless parameters.
In the Horndeski language, this is obtained by taking
$G_2=\beta X-V(\phi)$ and $G_3=\alphaM_{\rm Pl} X/U(\phi)$ in Eq.~\eqref{cubic-Horndeski}.

It is straightforward to see
from the field equations~\eqref{app:Fridmann}--\eqref{app:Scalar}
that the Lagrangian~\eqref{Lag:NEC-V}
with~\eqref{Pot:U=V} admits the desired NEC-violating solution of the form
\begin{align}
H=\frac{p}{-t},\quad \phi = M_{\rm Pl} \ln \left[\frac{M_{\rm Pl}}{\sqrt{q U_0}(-t)}\right],
\label{Soln:NEC-V}
\end{align}
where $p>0$ and $q>0$ are dimensionless parameters characterizing the solution.
These parameters are determined solely from $\alpha$ and $\beta$:
$p=p(\alpha,\beta)$, $q=q(\alpha,\beta)$.
Compact expressions are obtained if one instead writes $\alpha$ and $\beta$
in terms of $p$ and $q$:
\begin{align}
\alpha = \frac{2p}{q}-\frac{2}{q^2(1+3p)},
\quad
\beta =-2p(2+3p)+\frac{2}{q}.\label{Eq:abpq}
\end{align}
Figure~\ref{fig:stable.pdf} shows constant $\alpha$ and constant $\beta$ curves
in the $p$-$q$ plane. For given $(\alpha,\beta)$, one finds two
sets of $(p,q)$ satisfying $p,q>0$.

  \begin{figure}[tb]
    \begin{center}
            \includegraphics[keepaspectratio=true,height=95mm]{stable.pdf}
       \caption{Constant $\alpha$ curves (solid lines) and constant $\beta$ curves
       (dashed lines) in the $p$-$q$ plane. The two stability conditions are
       satisfied in the shaded region.
  	}
       \label{fig:stable.pdf}
    \end{center}
  \end{figure}

In order for the NEC-violating background~\eqref{Soln:NEC-V} to be stable,
it is required that
\begin{align}
{\cal F}_S&=-1+(-1-2p+3p^2)q>0,\label{Stab:1}
\\
{\cal G}_S&=2 p (1+3 p)^3 q^2-3 (1+3 p)^2 q+3>0.\label{Stab:2}
\end{align}
The two conditions are satisfied
in the shaded region in Fig.~\ref{fig:stable.pdf}.
A stable NEC-violating solution can be found for $\alpha>0$ and $\beta<-10$.
Note that all solutions with $p\le  1$ are unstable.
As we have mentioned above, we have two sets of $(p,q)$ for given $(\alpha,\beta)$.
It turns out that at least one of them leads to an unstable background.
(Since this unstable branch typically has $p\gg 1$ and $q\ll 1$,
it is difficult to judge the stability by seeing Fig.~\ref{fig:stable.pdf}.)

\subsection{The full Lagrangian}\label{subsec3b}

We have thus seen that the Lagrangian~\eqref{Lag:NEC-V} with~\eqref{Pot:U=V}
admits the NEC-violating solution~\eqref{Soln:NEC-V},
and there is a parameter region where the solution is stable
against scalar perturbations.
In this subsection, we modify the Lagrangian~\eqref{Lag:NEC-V}
so that it can be used to describe the subsequent inflation and reheating stages
as well as the NEC-violating phase.

The full Lagrangian is given by
\begin{align}
    {\cal L}=\frac{M_{\rm Pl}^2}{2}R+B(\phi)X-V(\phi)-\frac{\alphaM_{\rm Pl}}{U(\phi)}X\Box\phi,
    \label{Lag:full}
\end{align}
where we now take
\begin{align}
    V(\phi)&=\frac{U(\phi)\cdot m^2\phi^2}{2U(\phi)+m^2\phi^2},\label{Pot:full}
    \\
    B(\phi)&=\frac{1+\beta}{2}+\frac{1-\beta}{2}\tanh[c(\phi-\phi_*)/M_{\rm Pl}],\label{Func:kinatic}
\end{align}
with $U(\phi)=U_0e^{2\phi/M_{\rm Pl}}$.
The potential $V(\phi)$ is taken so that $V\simeq U$ for negatively large $\phi$
and $V\simeq m^2\phi^2/2$ for small $\phi$, as shown in Fig.~\ref{fig:potential.pdf}.
Any potential with a similar form is possible as well, but we
use this particular form for simplicity.
The potential peaks at $\phi=\phi_*\,(<0)$, where
\begin{align}
    \phi_*=-\frac{3}{2}M_{\rm Pl} W_0(z),
    \quad
    z=\frac{2^{4/3}}{3}\left(\frac{U_0}{M_{\rm Pl}^2m^2}\right)^{1/3},
    \label{Def:phi-ast}
\end{align}
with $W_0$ being the Lambert $W$ function (the principal branch of the inverse functions of $z=We^{W}$).
The coefficient of the kinetic term, $B(\phi)$, is taken so that
$B\simeq \beta \,(<0)$ for $\phi<\phi_*$ and $B\simeq 1$ for $\phi>\phi_*$.
Here, $c$ is some constant of ${\cal O}(1)$--${\cal O}(10)$
whose explicit value is not important.
The sign flip of the kinetic term soon after the end of the NEC-violating phase
is necessary for the subsequent inflation and reheating stages to be
free from ghost instability. This idea is analogous to what is supposed to happen
at the end of k-inflation~\cite{ArmendarizPicon:1999rj}.

The Lagrangian~\eqref{Lag:full} reduces to
Eq.~\eqref{Lag:NEC-V} with~\eqref{Pot:U=V} for negatively large $\phi$,
giving rise to the NEC-violating solution~\eqref{Soln:NEC-V}.
The sign of the kinetic term is flipped when $\phi$ goes over the top of the potential,
and inflation then starts
in the potential $V\simeq m^2\phi^2/2$
provided that $|\phi_*|\gtrsim M_{\rm Pl}$. After inflation,
reheating occurs in a usual way at the bottom of the $\phi^2$ potential.
We thus expect that the scenario outlined in the previous section is
realized by the Lagrangian~\eqref{Lag:full}.
This will be demonstrated explicitly by numerical calculations
in the next subsection.

We have four key quantities characterizing our model:
$p$ in the expansion law in the NEC-violating phase,
the maximum value of the Hubble parameter reached at the end of
the NEC-violating phase (denoted by $H_*$),
the inflationary Hubble scale (hereafter denoted by $H_{\rm inf}$),
and the duration of inflation (hereafter denoted by $\Delta{\cal N}_{\rm inf}$).
These quantities are determined by the four parameters in the Lagrangian~\eqref{Lag:full}
($\alpha$, $\beta$, $U_0$, and $m$).\footnote{The parameter $c$
plays essentially no role in the dynamics.}
As we have already seen, $p$ is determined by $\alpha$ and $\beta$
through the relation~\eqref{Eq:abpq}.
Assuming that the transition from the NEC-violating phase
to inflation is sudden, one can estimate as
\begin{align}
    H_*&\sim \frac{p\sqrt{q U_0}}{M_{\rm Pl}}e^{\phi_*/M_{\rm Pl}}=p\sqrt{\frac{q}{2}}m
    \left(\frac{|\phi_*|}{M_{\rm Pl}}\right)^{3/2},
    \label{Estimate:H-ast}
    \\
    H_{\rm inf}&\sim\frac{m|\phi_*|}{\sqrt{6}M_{\rm Pl}},
    \label{Estimate:H-inf}
\end{align}
where $\phi_*=\phi_*(U_0,m)$ was already defined in Eq.~\eqref{Def:phi-ast}.
In deriving the estimate~\eqref{Estimate:H-ast}, we used
$H_*\sim p/(-t_*)$ and $\phi_*\sim M_{\rm Pl} \ln[M_{\rm Pl}/\sqrt{qU_0}(-t_*)]$,
where $t_*$ is the time at the end of the NEC-violating phase.
Finally, a simple textbook calculation yields
\begin{align}
    \Delta{\cal N}_{\rm inf}\sim \frac{1}{4}\left(\frac{|\phi_*|}{M_{\rm Pl}}\right)^{2}.
    \label{Estimate:N-inf}
\end{align}

As seen from Eq.~\eqref{Estimate:N-inf},
we can adjust the number of
e-folds of inflation by shifting the top of the potential $\phi_*$.
To demonstrate how $\phi_*$ depends on $U_0$,
it is convenient to introduce a new parameter $\phi_0$ defined as
$\phi_0:=(1/2)M_{\rm Pl}\ln(U_0/M_{\rm Pl}^4)$
(i.e., $U_0 = M_{\rm Pl}^4 e^{2 \phi_0 / M_{\rm Pl}}$).
The shift $\phi_0\to\phi_0+c$ with $c>0$ ($c<0$)
causes left (right) shifting of $U(\phi)$.
The Lambert $W$ function $W_0(z)$ can be approximated roughly by
$W_0(z) \approx \ln{(z/\ln{z})}$ or more crudely $W_0(z) \sim \ln{z}$ for large $z>0$.
Using the latter approximation here in Eq.~\eqref{Def:phi-ast}, we get
$\phi_* \sim -\phi_0 + {\rm const.}$
We typically take $\phi_0/M_{\rm Pl} \sim \mathcal{O}(10)$ to sustain
a sufficient duration of inflation.
One might suspect that such a value of $\phi_0$ yielded extremely large $U_0$
and the universe experienced super-Planck scale $H$. However,
since $M_{\rm Pl}^2H^2\lesssim V(\phi_*) = (1/2) m^2 \phi_*^3 / (\phi_*-M_{\rm Pl})$,
it is not the case as long as $\left| \phi_* \right| \ll M_{\rm Pl}^2/m$.
Hereafter we use $\phi_0$ instead of $U_0$.

\subsection{Numerical results}

  \begin{figure}[tb]
    \begin{center}
            \includegraphics[keepaspectratio=true,height=45mm]{attractor.pdf}
       \caption{Typical phase-space trajectories. The Friedmann constraint
       does not have real roots for $H$ in the shaded region.
  	}
       \label{fig:attractor.pdf}
    \end{center}
  \end{figure}

Let us first affirm that the NEC-violating solution~\eqref{Soln:NEC-V}
is an attractor. Figure~\ref{fig:attractor.pdf} shows the behavior
of the solutions for $(\alpha,\beta)= (1.09,-12.4)$ with different initial conditions,
obtained by solving the dynamical equations~\eqref{app:Evolution}
and~\eqref{app:Scalar} numerically.\footnote{The initial value for $H$ is determined from
the Friedmann constraint~\eqref{app:Fridmann}. We check that the constraint~\eqref{app:Fridmann}
is satisfied at each time step.}
We see that the trajectories converge to either of the two curves,
$\dot\phi\propto \sqrt{q}e^{\phi/M_{\rm Pl}}$, with different $q$.
In the upper left region the solutions converge to
the stable NEC-violating solution with $(p,q)\simeq (1.2,2)$, while in the lower right region
the solutions converge to the unstable background with $(p,q)\simeq (6.1,8.4\times 10^{-3})$.
The two regions of initial conditions are separated by the ``no-go'' region
(the shaded region in Fig.~\ref{fig:attractor.pdf})
where the Friedmann constraint equation does not have real roots for $H$.
In the following analysis, we assume that the initial conditions for $(\phi,\dot\phi)$
are in the upper left region so that they lead to a stable NEC-violating universe.

We then investigate the entire evolution numerically.
In Figure~\ref{fig:time-evolution-of-vars},
we plot the time evolution
of $\phi$, $H$, $\mathcal{F}_S$, and $c_s^2 := \mathcal{F}_S/ \mathcal{G}_S$
for the parameter choice
$p=3/2$ ($n_t=0.8$), $q=10$, $m^2=10^{-6}M_{\rm Pl}^2$, and $\phi_0=12M_{\rm Pl}$.
In the NEC-violating phase ($0<\ln a < 26.0$), we have positive $\mathcal{F}_S$.
Note that $H_*$ is as high as $\mathcal{O}(0.1)\times M_{\rm Pl}$.
In the short decelerating phase ($26.0 < \ln a < 27.3$),
$\mathcal{F}_S$ becomes negative,
which does not contradict the stability condition~\eqref{Stab:1}
since the system is no longer on the attractor~\eqref{Soln:NEC-V}.
This point will be discussed shortly.
In the inflationary phase ($27.3 < \ln a \lesssim 67$), $\mathcal{F}_S$ is approximately given by
$- {\rm d}  \ln H / {\rm d}  \ln a$, which is positive in an NEC-preserving universe.
At $\ln a =67.7$, $\phi$ passes the origin and starts to oscillate,
leading to the expansion of a matter-dominated universe.

  \begin{figure}[tb]
    \begin{center}
        \includegraphics[keepaspectratio=true,height=100mm]{lna_vars.pdf}
       \caption{
       Time evolution of $\phi$, $H$, $\mathcal{F}_S$, and $c_s^2 (= \mathcal{F}_S/ \mathcal{G}_S)$.
       The model parameters are chosen as $p=3/2$ ($n_t=0.8$), $q=10$, $m^2=10^{-6}M_{\rm Pl}^2$,
       and $\phi_0=12M_{\rm Pl}$.
       In the plots,
       the dashed lines denote the negative value having the same absolute value.
  	}
       \label{fig:time-evolution-of-vars}
    \end{center}
  \end{figure}

Our numerical result implies a potential drawback of the present model
that it exhibits a gradient instability in the short intermediate stage
between the NEC-violating and inflationary phases.
This is similar to what happens at the transition from
the galilean genesis phase to inflation~\cite{Pirtskhalava:2014esa},
and in fact the appearance of gradient instabilities is generic to
non-singular cosmological solutions in galileon/Horndeski
theories~\cite{Libanov:2016kfc,Kobayashi:2016xpl}.
In the present case, however, the instability could in principle be
avoidable by more elaborate choices of the Lagrangian and model parameters,
because the background~\eqref{eq:Hubble-NECV} is
singular in the $t\to-\infty$ limit\footnote{{Since
$\int^t_{-\infty}a{\rm d} t<\infty$ for $p>1$, this NEC-violating universe
is past incomplete. In this sense, the present model is
as bad as inflation.}}
 and hence the ``no-go''
theorems in~\cite{Libanov:2016kfc,Kobayashi:2016xpl} do not apply.
The instability issue can also be resolved
by introducing ``beyond-Horndeski'' or higher-order
operators~\cite{Creminelli:2016zwa,Cai:2016thi,Cai:2017tku,Kobayashi:2015gga}.

\section{Gravitational Wave Spectrum}\label{sec4}

Since we do not consider the scalar field derivatively coupled to
the curvature through the so-called
Horndeski $G_4(\phi,X)$ and $G_5(\phi, X)$ terms,
the tensor perturbations $h_{ij}$ evolve following the
same equation as in general relativity.
Their equation of motion is
\begin{align}
    h_{ij}'' + 2 (a'/a) h_{ij}'-\partial^2 h_{ij} = 0,
\end{align}
where a prime denotes the derivative with respect to
the conformal time $\eta$.
We move to Fourier space, expand $h_{ij}$ into the basis
of the polarization tensors, and write
\begin{align}
    h_{ij}(\eta,\Vec{x})=\frac{1}{(2\pi)^3}\sum_{A=+,\times}\int{\rm d}^3k \,
    \tilde h_A(\eta,\Vec{k})e^A_{ij}(\Vec{k})e^{i\Vec{k}\cdot\Vec{x}},
\end{align}
where $e^A_{ij}$ is transverse, traceless, and normalized as $e^A_{ij}e^{A'}_{ij}=\delta^{A}_{A'}$.
Following the usual quantization procedure,
we obtain the canonically normalized positive frequency mode
in the NEC-violating phase as
\begin{align}
    \tilde h_A=\frac{\sqrt{\pi}}{ a M_{\rm Pl} } \sqrt{-\eta} H^{(1)}_{\nu} (-k\eta),
    \label{posfreqmode}
\end{align}
where $H^{(1)}_{\nu}$ is the Hankel function and $\nu := (3p+1)/2(p+1)$.
The initial conditions for the tensor perturbations
are set by Eq.~\eqref{posfreqmode}.
In the superhorizon regime, $|k\eta|\ll 1$, we have
\begin{align}
    k^3|\tilde h_A|^2\simeq \frac{2^{2\nu}\Gamma^2(\nu)}{\pi M_{\rm Pl}^2}
    \left(\frac{2H_c}{2\nu-1}\right)^{2\nu-1}k^{3-2\nu}=
    {\rm const},
\end{align}
where we introduced the constant
$H_c:=Ha^{-1/p}$.
From this, we get the primordial power spectrum of gravitational waves as
\begin{align}
    \mathcal{P}_{{\rm GW}}
    =2\cdot \frac{k^3}{2\pi^2}|\tilde h_A|^2
                    = \frac{2^{2\nu} \Gamma(\nu)^2}{\pi^3 M_{\rm Pl}^2}
                    \! \left( \frac{2 H_c}{2\nu-1} \right)^{\! 2\nu-1} \!\! k^{3-2\nu} ,
                    \label{pgw}
\end{align}
Note that for $\nu=3/2$
this reproduces the usual primordial power spectrum in de Sitter inflation.
From Eq.~\eqref{pgw} we see that
the modes which exit the horizon in the NEC-violating phase
have a primordial spectrum with a blue spectral tilt~\eqref{eq:tensor-index},
whereas those which exit the horizon in the subsequent inflationary phase have
an almost flat primordial spectrum with $n_t =-2\epsilon\, (<0)$,
where $\epsilon$ is the usual small slow-roll parameter,
$\epsilon:=-\dot H/H^2$.
After re-entering the horizon their amplitudes decay proportionally to $1/a$.

The present energy density spectrum of gravitational waves is given by
\begin{align}
    \Omega_{\rm GW}(k) =\frac{k^5 \left(
    |\tilde{h}_{+} |^2 + |\tilde{h}_{\times} |^2\right) }{24 \pi^2 a^2 H_0^2},
    \label{Def:GWspectrum}
\end{align}
and below we show $\Omega_{\rm GW}(f)$ plots, where $k$ is the comoving wavenumber,
$f$ is the physical frequency $f = k/2\pi a$,
and $H_0$ is the Hubble constant.
Each spectrum in our model is typically composed by the three parts;
an extremely blue-tilted spectrum in the low frequency range;
{wiggles} in the middle frequency range
(due to the intermediate decelerating stage);
and a slightly red-tilted spectrum in the high frequency range.
We plot the energy density spectra
for different values of $\phi_0$ in Fig.~\ref{fig:Energy_phi0}.
A larger value of $\phi_0$ leads to a larger number of e-folds, a lower peak frequency,
and as a result larger $\Omega_{\rm GW}(f)$ on CMB scales.
Note in passing that here and hereafter we assume that
the Hubble parameter at which the reheating process becomes so efficient
that the energy is converted from matter to radiation rapidly, $H_{\rm R}$,
is given by $H_{\rm R}=10^{-8}M_{\rm Pl}$.

  \begin{figure}[tb]
    \begin{center}
            \includegraphics[keepaspectratio=true,height=55mm]{f_Omega_phi0.pdf}
       \caption{Gravitational wave spectra in the cases of
       $\phi_0/M_{\rm Pl}=12$ (solid), $10$ (left dashed), $8$ (left dotted),
       $6$ (right dashed), and $4$ (right dotted).
       The other parameters are the same as in Fig.~\ref{fig:time-evolution-of-vars}.
  	}
       \label{fig:Energy_phi0}
  	\end{center}
  \end{figure}

Let us study the effects of changing
the other parameters $m$, $q$, and $p$.
Since we have $|\phi_*|/M_{\rm Pl} \sim \mathcal{O}(1-10)$,
$m$ is the only parameter which determines
the order of an energy scale of the slow-rolling phase,
although $\phi_*$ also has an implicit logarithmic dependence on $m$ and $\phi_0$.
We plot the energy density spectra for different $m$ in Fig.~\ref{fig:Energy_m}.
  \begin{figure}[tb]
    \begin{center}
        \includegraphics[keepaspectratio=true,height=55mm]{f_Omega_m.pdf}
       \caption{Gravitational wave spectra in the cases of
       $m^2/M_{\rm Pl}^2=10^{-6}$ (solid), $10^{-8}$ (dashed), and $10^{-10}$ (dotted).
       In each spectrum, we set $\phi_0/M_{\rm Pl} = 12+\ln{(10^3 m/M_{\rm Pl})}$ to have nearly constant $\phi_*$.
       The other parameters are the same as in Fig.~\ref{fig:time-evolution-of-vars}.
  	}
       \label{fig:Energy_m}
    \end{center}
  \end{figure}
From Eq.~\eqref{Estimate:H-ast}, once we set $qm^2$ constant,
we get similar values of $H_*$ even with different $m$.
We show energy density spectra in the case of different $m$ but the same $qm^2$
in Fig.~\ref{fig:Energy_qm}.
  \begin{figure}[tb]
    \begin{center}
        \includegraphics[keepaspectratio=true,height=55mm]{f_Omega_qm.pdf}
       \caption{Gravitational wave spectra in the cases of
       $m^2=10^{-5}$ (dashed), $10^{-6}$ (solid), and $10^{-7}$ (dotted).
       The parameter $q$ are taken as $q m^2=10^{-5}$
       in order to produce similar blue-tilted spectra.
       The other parameters are the same as in Fig.~\ref{fig:time-evolution-of-vars}.
  	}
       \label{fig:Energy_qm}
    \end{center}
  \end{figure}

To compare our results with the NANOGrav observation~\cite{Arzoumanian:2020vkk},
we parametrize the characteristic strain of gravitational waves
around the reference frequency $f_{\rm yr}=1/{\rm yr}$ as
\begin{align}
    h_{c}(f) = A_{\rm GW} \left( \frac{f}{f_{\rm yr}} \right)^{\alpha},
\end{align}
which is related to $\Omega_{\rm GW}$ as $\Omega_{\rm GW} = 2\pi^2 f^2 h_c^2 / (3H_0^2)$.
The relation between $\alpha$ and the slope of the cross-power spectral density $\gamma=3-2\alpha$
is also useful for direct comparison to the observational implication.
Since $\Omega_{\rm GW}(f) \propto f^{n_t}$ is observed
for the modes which re-enter the horizon in the radiation-dominated era,
we get $\alpha = (n_t-2)/2 = -{p}/{(p+1)}$ from Eq.~\eqref{eq:tensor-index},
if the extremely blue-tilted gravitational waves from the NEC-violating phase
give rise to the spectrum around $f_{\rm yr}$.
In Fig.~\ref{fig:Energy_nt95}, we plot the
energy density spectra for different values of $p$.
We also show the sensitivity curves for present and future experiments
of pulsar timing arrays and space/ground-based interferometers.
For the plotted spectra,
their spectral slopes and amplitudes are
$(\gamma, \log_{10} A_{\rm GW})= (4.2,-14.525)$,
$(4.1 , -14.542)$, and $(4.05 , -14.535)$, for $n_t=0.8$, $0.9$, and $0.95$, respectively,
which are within the 2-$\sigma$ posterior contour
for the five frequency power-law model
of a spatially-uncorrelated common-spectrum process~\cite{Arzoumanian:2020vkk}.

  \begin{figure*}[tb]
    \begin{center}
            \includegraphics[keepaspectratio=true,height=100mm]{f_Omega_nt95.pdf}
       \caption{Energy density spectra of gravitational waves in the cases of $n_t=0.8$ (solid), $0.9$ (dashed), and $0.95$ (dot-dashed),
       which correspond to $p=3/2$, $11/9$, and $21/19$, respectively.
       The other parameters are the same as in Fig.~\ref{fig:time-evolution-of-vars}.
       We give the power-law-integrated sensitivity curves~\cite{Schmitz:2020syl} for present and future experiments of pulsar-timing arrays
       and space/ground-based interferometers.
       For comparison with the lower-frequency constraint, we also plot as the gray line
       the spectrum with $n_t=0$ and $r=0.07$,
       the upper limit from the {\it BICEP2/Keck Array} and {\it Planck} experiments~\cite{Ade:2018gkx},
       between the frequencies corresponding to multipoles $\ell=2$ and $10$.
  	}
       \label{fig:Energy_nt95}
    \end{center}
  \end{figure*}

Finally, let us discuss the BBN constraint~\cite{Allen:1996vm,Maggiore:1999vm,Smith:2006nka,Boyle:2007zx,Kuroyanagi:2014nba} in our model.
We calculate the total energy density of gravitational waves for the case with $n_t=0.8$ in Fig.~\ref{fig:Energy_nt95}
and obtain its contribution to the effective number of relativistic species $\Delta N_{\rm eff}$ as
\begin{align}
     \Delta N_{\rm eff} \simeq \int {\rm d}
     \ln{\!f} \; \frac{\Omega_{\rm GW}h^2}{5.6 \times 10^{-6}} \simeq  1.3 \times 10^{-3}.
\end{align}
This is much smaller than the upper bound from the BBN constraint,
$\Delta N_{\rm eff} < \mathcal{O}(0.1)$~\cite{Cooke:2017cwo,Hsyu:2020uqb}.
(The similar bound
also comes from the large scale observations~\cite{Aghanim:2018eyx,Aiola:2020azj}.)
Since we have the nearly flat part of the spectrum in the high-frequency range,
we can avoid the problem pointed out by~\cite{Vagnozzi:2020gtf},
in which a single power-law, blue-tilted spectrum violates the BBN constraint
if it is consistent with both the CMB constraint and the NANOGrav result.

\section{Scalar Perturbations}\label{sec5}

Large-scale curvature perturbations must be generated
in the early universe to account for the observed nearly scale-invariant
spectrum of fluctuations. In our scenario, however, it cannot be due to
the fluctuations of the $\phi$ field, because in the NEC-violating phase
we have ${\cal G}_S=\,$const and ${\cal F}_S=\,$const, which means that
the curvature perturbation $\zeta$ has a highly blue spectrum
with $n_s-1=n_t>0$ and hence is irrelevant on CMB scales.

To generate nearly scale-invariant curvature perturbations on CMB scales,
we need another scalar field. Since the background in the NEC-violating phase
is significantly away from de Sitter,
a minimally coupled scalar field does not give rise to a scale-invariant spectrum
of fluctuations.
We therefore extend the mechanism proposed originally
in the context of the galilean genesis scenario~\cite{Creminelli:2010ba}
to a more general background, and
consider a non-minimally coupled scalar field,
\begin{align}
    S_\sigma=\int{\rm d}^4x\sqrt{-\widetilde g}\left[
    -\frac{1}{2}\widetilde g^{\mu\nu}\partial_\mu\sigma\partial_\nu\sigma-{\cal V}(\sigma)
    \right],
\end{align}
with
\begin{align}
    \widetilde g_{\mu\nu}=e^{2\phi/M_{\rm Pl}}g_{\mu\nu}.\label{confomal-metric}
\end{align}
The conformally related metric~\eqref{confomal-metric} reads more explicitly
\begin{align}
    \widetilde g_{\mu\nu}{\rm d} x^\mu{\rm d} x^\nu =
    \widetilde a^2\left(-{\rm d}\eta^2+{\rm d}\Vec{x}^2\right),
\end{align}
where
\begin{align}
    \widetilde a =e^{\phi/M_{\rm Pl}}a \propto \frac{1}{(-\eta)},
\end{align}
in the NEC-violating phase, irrespective of the value of $p$.
This means that the non-minimally coupled scalar field $\sigma$
effectively lives in de Sitter, generating a nearly scale-invariant spectrum.
We therefore expect that fluctuations of the $\sigma$ field are the origin
of the observed spectrum of scalar perturbations.
To obtain a more detailed prediction, we need to investigate the
process of converting the fluctuations of $\sigma$ to adiabatic modes,
which is model-dependent and deserves further study.

\section{Conclusions}\label{sec6}

In this paper,
we have explored the possibility to generate nanohertz gravitational waves
consistent with both of the CMB and BBN constraints.
We investigated the scenario which results in an extremely blue-tilted spectrum appropriate for
the CMB constraint at $f \lesssim 10^{-15} \;{\rm Hz}$ and the NANOGrav result at $f \sim 1/{\rm yr}$ and a (nearly) flat spectrum for $f\gtrsim 1$/yr
to evade the BBN constraint.
To produce an extremely blue spectrum,
we considered an NEC-violating phase in the early universe
undergoing the (super-)accelerated expansion with $a\propto (-t)^p$ ($p>0$),
in which the tensor spectral index $n_t>0$
is related to the power-law index $p$ as in Eq.~\eqref{eq:tensor-index}.
To obtain a flat spectrum in a high frequency range,
we examined the scenario in which slow-roll inflation follows after the NEC-violating phase and ends with standard reheating,
as sketched out in the potential of Fig.~\ref{fig:potential.pdf}.

In Sec.~\ref{sec3},
we proposed a successful model in the cubic Horndeski theory.
First, we gave the minimal action~\eqref{Lag:NEC-V} which can cause a stable NEC-violating phase.
We obtained the relation between the model parameters $(\alpha,\beta)$ and the solution parameters $(p,q)$ in Eq.~\eqref{Eq:abpq},
and the stability condition in terms of $(p,q)$ in Eqs.~\eqref{Stab:1} and \eqref{Stab:2}.
From these conditions, we found that a solution with $p \le 1$ is always unstable within our model,
while stable solutions with $p>1$ are indeed possible.
We expect that a stable model with $p \le 1$ can be realized by adding other higher-order terms in the action,
but the analysis would be more complicated.
Next, we assume that the potential
and the coefficient of the kinetic term are given by Eq.~\eqref{Func:kinatic},
which changes the field dynamics for $\phi > \phi_*$,
connecting the NEC-violating phase to the subsequent inflationary phase.
The phenomenological parameters $(H_\ast, H_{\rm inf}, \Delta{\cal N}_{\rm inf})$ are estimated in
Eqs.~\eqref{Estimate:H-ast}, \eqref{Estimate:H-inf}, and \eqref{Estimate:N-inf},
enabling us to design the gravitational wave spectrum.
In Fig.~\ref{fig:time-evolution-of-vars},
we numerically gave the typical time evolution of the background
and showed that a gradient instability occurs in the intermediate decelerating phase.
We have discussed that this issue can be resolved by higher-order operators such as
in beyond-Horndeski theories~\cite{Creminelli:2016zwa,Cai:2016thi,Cai:2017tku,Kobayashi:2015gga}.

In Sec.~\ref{sec4},
we have presented the gravitational wave spectrum~\eqref{Def:GWspectrum} in our model with different parameter choices.
As is expected, those spectra are composed of
the blue-tilted part and the nearly flat part,
of which amplitudes can be
controlled almost independently of each other.
We gave gravitational wave spectra with $n_t = 0.8$, $0.9$, and $0.95$
for $f\lesssim 1/$yr
in Fig.~\ref{fig:Energy_nt95}, which are consistent with the recent NANOGrav result~\cite{Arzoumanian:2020vkk}.
We found that thanks to the flat part of the spectrum we are no longer bothered by the BBN bound and
this flat part is possibly detected by future space/ground-based observations.
In Sec.~\ref{sec5}, we discussed how the nearly scale-invariant
scalar perturbations observed by CMB experiments are produced
by introducing another non-minimally coupled scalar field~\eqref{Lag:full}
which effectively lives in de Sitter even in the NEC-violating phase.

In summary, we have obtained an explicit early universe model
that can yield an extremely blue gravitational wave spectrum peaking
at $f\sim 1/$yr explored by pulsar timing arrays, while being consistent
with other observational constraints.

\acknowledgments
We thank Kohei Kamada and Sachiko Kuroyanagi for useful discussions.
HWHT was supported by the Grant-in-Aid for JSPS Research Fellow No.~JP19J11640.
TK was partially supported by JSPS KAKENHI Grant Nos.~JP20H04745 and~JP20K03936.