Common Lagrangian in introductory courses of classical mechanics depends of $q$
and $\dot{q}$. We refer to them as \emph{first-order Lagrangian} because the
highest derivative term is $\dot{q}$, i.e.\ the first derivative of $q(t)$ with
respect to time. If in the Lagrangian appears also $\ddot{q}$ it will be a
\emph{second-order Lagrangian} and so on. \\

In the following, the mathematical machinery develop in Section~\ref{section:
lagrangin and hamiltonian formalism}, will be used to study two apparently
similar Lagrangian and introduce the core problem of HD systems.

\subsection{First-order vs Second-order Lagrangian}
The simplest non-trivial Lagrangians are those in polynomial form like
\begin{equation*}
  L_1(q, \dot{q}) = \frac{\dot{q}^2}{2} - \alpha(q) \qquad
  L_2(q, \ddot{q}) = \frac{\ddot{q}^2}{2} - \alpha(q)
\end{equation*}
where $\alpha$ is a smooth function of $q$ only~\cite{Chen13}. No Greek indices
appear on $q^{(i)}$ because we limit ourself to one-dimensional systems. $L_1$
is a first-order Lagrangian while $L_2$ it's a second-order one.

\paragraph{First-order Lagrangian} In order to get the equation of motion we
used the Euler-Lagrange equation~\eqref{eq:euler-lagrange}

\begin{equation} \label{eq: first-order_motion_eq_lagr}
  \frac{\partial L_1}{\partial q} -
  \frac{d}{dt}\frac{\partial L_1}{\partial \dot{q}} +
  \underbrace{
  \frac{d^2}{dt^2}\frac{\partial L_1}{\partial \ddot{q}} - \ldots }_0 = 0
  \qquad \Rightarrow \qquad
  \ddot{q} = - \frac{d\alpha(q)}{dq}
\end{equation}

If $q$ is the position coordinate of a unit mass point, the function $\alpha(q)$
can be interpreted as the potential energy and~\eqref{eq:
first-order_motion_eq_lagr} is nothing less than Newton's second law. To get the
Hamiltonian description apply the recipe in Section~\ref{subsection:
hamiltonian_formalism}

\begin{equation*}
  \begin{cases}
    Q_1 := q \\
    P_1 := \frac{\delta L_1}{\delta \dot{q}} = \dot{q}
  \end{cases}
  \qquad
  \tilde{H}_1(q, \dot{q}) :=
  \frac{\delta L_1}{\delta \dot{q}} \dot{q} - L_1(q, \dot{q})
\end{equation*}

\begin{equation} \label{eq: first-order_motion_eq_ham}
  H_1(Q_1, P_1) = \frac{P_1^2}{2} + \alpha(Q_1) \qquad
  \begin{cases}
    \dot{Q_1} =   \frac{\partial H_1}{\partial P_1} = P_1 \\
    \dot{P_1} = - \frac{\partial H_1}{\partial Q_1} = - \frac{d\alpha(Q_1)}{dQ_1}
  \end{cases}
\end{equation}

Since $L_1$ does not explicitly depend on time (only through $q$ and $\dot{q}$),
the Hamiltonian $H_1$ can be interpreted as the total energy of the system (sum
of kinetic and potential energy). Combining the two equation in~\eqref{eq:
first-order_motion_eq_ham} we restore~\eqref{eq: first-order_motion_eq_lagr}
proving the equivalence between Lagrangian and Hamiltonian mechanics for simple
systems describe by $L_1$.


\paragraph{Second-order Lagrangian} Proceeding as before, equation of motion
follow from Euler-Lagrange equation~\eqref{eq:euler-lagrange}

\begin{equation} \label{eq: second-order_motion_eq_lagr}
  \frac{\partial L_2}{\partial q} -
  \frac{d}{dt}\frac{\partial L_2}{\partial \dot{q}} +
  \frac{d^2}{dt^2}\frac{\partial L_2}{\partial \ddot{q}} -
  \underbrace{
    \frac{d^3}{dt^3}\frac{\partial L_2}{\partial q^{(3)}} + \ldots
  }_0 = 0
  \quad \Rightarrow \quad
  q^{(4)} = - \frac{d\alpha(q)}{dq}
\end{equation}

Equation~\eqref{eq: second-order_motion_eq_lagr} never appear in classical
mechanics, but can it describe some kind of system? The study of the Hamiltonian
can give us some clues about the total energy of the system.

\begin{equation*}
  \begin{cases}
    Q_1 := q \\
    P_1 := \frac{\delta L_2}{\delta \dot{q}}
         = \frac{\partial L_2}{\partial \dot{q}} -
           \frac{d}{dt} \left( \frac{\partial L_2}{\partial \ddot{q}} \right)
         = 0 - \frac{d}{dt} \left( \ddot{q} \right) = -q^{(3)}
  \end{cases}
  \begin{cases}
    Q_2 := \dot{q} \\
    P_2 := \frac{\delta L_2}{\delta \ddot{q}}
         = \frac{\partial L_2}{\partial \ddot{q}} = \ddot{q}
  \end{cases}
\end{equation*}

Therefore
\begin{equation}\label{eq:canonical_coordinates_second_order_lagrangian}
  \begin{cases}
    q        = Q_1 \\
    \dot{q}  = Q_2 \\
    \ddot{q} = P_2
  \end{cases}
  \qquad
  \begin{cases}
    \frac{\delta L_2}{\delta \dot{q}}  = P_1 \\
    \frac{\delta L_2}{\delta \ddot{q}} = P_2
  \end{cases}
\end{equation}

\begin{equation*}
  \tilde{H}_2(q, \dot{q}) :=
    \frac{\delta L_2}{\delta \dot{q}} \dot{q} +
    \frac{\delta L_2}{\delta \ddot{q}} \ddot{q} -
    L_2(q, \ddot{q})
    = \frac{\delta L_2}{\delta \dot{q}} \dot{q} +
    \frac{\delta L_2}{\delta \ddot{q}} \ddot{q} -
    \frac{\ddot{q}^2}{2} + \alpha(q)
\end{equation*}

Substituting relations~\eqref{eq:canonical_coordinates_second_order_lagrangian}
in the expression for $\tilde{H}_2$ we get the Hamiltonian

\begin{align} \label{eq: second-order_motion_eq_ham}
  H_2(Q_1, Q_2, P_1, P_2)
    &= P_1 Q_2 + P_2 P_2 - \frac{P_2^2}{2} + \alpha(Q_1) \notag \\
    &= P_1 Q_2 + \frac{P_2^2}{2} + \alpha(Q_1)
\end{align}

There is a significant difference between the two spectra of $H_1$ and $H_2$:
the first is bounded from below while the latter is not (and both are not
bounded from above).


\subsection{Linear Ostrogradskian instability}~\label{subsection:
linear_ostrogradskian_instability}
Consider an \emph{isolated} system describe by the Hamiltonian $H_2$; if the
energy is conserve even though the spectra is not bounded the energy stay
constant. Things start to going wrong when ones consider \emph{interacting}
system (e.g. $H_1$-system interacting with $H_2$-system). The phase space of
$H_2$-system infinitely extend where the Hamiltonian is negative. In the example
this is due to the \emph{linear} term $P_1$. When two systems interact the
$H_1$-system tends to occupy higher and higher energies states while the
$H_2$-system fall lower and lower in negative energies; this behavior entails
the conservation of energy~\cite{Kallosh08, Eliezer89}.
This is the so called \emph{linear Ostrogradskian instability}. \\

One can encounter the same kind of instability in quantum mechanics when try to
canonically quantized higher order Hamiltonian (e.g. $H_2$). This process lead
to \emph{negative norm states} (or \emph{negative energy states}) which are
often called \emph{``ghosts''}. Systems (even classical ones) involving this
type of ghosts are called \emph{ghost-like} or \emph{ghost-ridden} systems. \\

One may wonder when these ghosts are summoned; the answer is contains in the
Ostrogradsky theorem

\begin{theorem}[Ostrogradsky]\label{th:ostrogradsky_classical}
  If in the Lagrangian~\eqref{eq:general_lagrangian} $n \geq 2$ and the
  canonical momentum $\bm{P}_n$ does not vanish, the
  Hamiltonian~\eqref{eq:general_hamiltonian} may acquire an arbitrary real
  value.
\end{theorem}
\begin{proof}
  The proof is given for a one-dimensional system (i.e. $\mu=1$) but it is
  generalize straight to higher dimensions. Consider the
  Hamiltonian~\eqref{eq:general_hamiltonian} in terms of conjugate coordinates
  \begin{equation} \label{eq:general_hamiltonian_1dim}
    H = P_n h + P_{n-1} Q_n + \cdots + P_1 Q_2 - L(Q_1, Q_2, \ldots, h)
  \end{equation}
  The function $h$ express $q^{(n)}$ in terms of conjugate coordinates under
  the assumption of regular Lagrangian, i.e.
  \begin{equation*}
    q^{(n)} = h(Q_1, \ldots, Q_n, P_n)
  \end{equation*}
  Observe that $h$ does not depend on $P_1, \ldots, P_{n-1}$ so the only
  contribution of these momenta to the Hamiltonian is in the linear factor
  \begin{equation*}
    P_{n-1} Q_n + \cdots + P_1 Q_2
  \end{equation*}
  which is linear in momenta and hence can acquire any real value and so does
  the Hamiltonian~\eqref{eq:general_hamiltonian_1dim}.
\end{proof}

For many years, due to consequences of Theorem~\ref{th:ostrogradsky_classical},
HD theories are consider intrinsically sick and not worthy of further study.
However a deeper investigation reveals that Ostrogradsky instability can be cure
and HD theories should not be discard in the first instance. One way to
exorcising Ostrogradsky's ghost is through imposition of
constraints~\cite{Chen13}.


\subsection{Removing Ostrogradskian instability with constraints}

In Hamiltonian mechanics, \emph{constraints} are relation between coordinates
and conjugate momenta. Constraints are classify in two group: \emph{first class
constraints} and \emph{second class constraints}. Second class constraints can
be thought as ``physical'' (e.g.\ train on a railway, particle on a plane) so
that the solutions of equations of motion are different with or without them.
First-class are instead related to a gauge freedom: the solutions of equations
of motion differ by a function of time and hence represent the same system.
Moreover the Poisson Bracket of a first class constraint with all the other
constraints vanishes on the constraints surface in the phase space (i.e.
$\{\phi, \tilde{\phi}\} \approx 0$, see later for the notation used).


Another independent classification can be made for constraints (and not confused
with the previous one): \emph{primary constraints} and \emph{secondary
constrains}. Primary constraints are relation between $Q$s and $P$s which are
given at the beginning.

\begin{equation} \label{eq:constraint}
  \phi_1(Q, P) = 0
\end{equation}

As the name suggests, secondary constraints are derived from the primary ones
imposing the conservations of the constraints during the evolution of the
system. This generates a series of constraints where the constraints relations
are called \emph{consistency relations}.

\begin{equation} \label{eq:consistency_relations}
  \phi_1 \approx 0 \quad \Rightarrow \quad
  \left\{ \phi_1 , H \right\} =: \phi_2 \approx 0 \quad \Rightarrow \quad
  \left\{ \phi_2 , H \right\} =: \phi_3 \approx 0 \quad \Rightarrow \quad
  \ldots
\end{equation}

The weak equality symbol ``$\approx$'' has been used to highlight the fact that
these relations vanish only on the hypersurface where all constraints are
satisfied. When we will be said that a constrain $\phi_i$ \emph{vanishes} it is
to be intended in this weak sense.

\subsubsection{Lagrangian and Hamiltonian with $m$ auxiliary variables}
Constraints in Hamiltonian mechanics can be imposed with the help of auxiliary
variables $\lambda_i$ in the Lagrangian:

\begin{equation*}
  L = L(
  q, \dot{q}, \ddot{q}, \ldots, q^{(n)},
  \lambda_1, \lambda_2, \ldots, \lambda_m
  )
\end{equation*}

So the Lagrangian, besides the equations of motion~\eqref{eq:euler-lagrange},
need to satisfy the constraints equation
\begin{equation}
  \frac{\partial L}{\partial \lambda_i} = 0 \qquad i=1, 2, \ldots, m
\end{equation}

Canonical coordinates have to be given also for $\lambda_i$ and the following
choice can be done
\begin{equation} \label{eq:def_canonical_coordinates_lambda}
  \Lambda_{i}:= \lambda_{i}
  \quad \leftrightarrow \quad
  \Pi_{i} := \frac{\delta L}{\delta \dot{\lambda_{i}}} = 0
  \qquad i = 1, 2, \ldots, m
\end{equation}
so the primary constraints in~\eqref{eq:constraint} simply become
$\phi_{1, i}: \Pi_i = 0$. To the Hamiltonian~\eqref{eq:Ham_in_q} are now added
the auxiliary variables terms and it becomes
\begin{equation} \label{eq:Ham_constraints_in_q}
  \tilde{H} :=
  \sum_{j=1}^{n} \frac{\delta L}{\delta q^{(j)}} q^{(j)} +
  \sum_{i=1}^{m} \frac{\delta L}{\delta \dot{\lambda_i}} \dot{\lambda_i} -
  L(q, \ldots, q^{(n)}, \lambda_1, \ldots, \lambda_m)
\end{equation}

The $\dot{\lambda_i}$ can be obtain by the consistency
relations~\eqref{eq:consistency_relations} but there is no need to explicitly
calculate them. $\dot{\lambda_i}$ can be written as functions of canonical
coordinates: $\dot{\lambda_i} = u_i(Q_1, \ldots, Q_n,  P_1, \ldots, P_n)$.
Expressing~\eqref{eq:Ham_constraints_in_q} using canonical
coordinates we obtain the Hamiltonian
\begin{align} \label{eq:Ham_constraints}
  H =\ & P_n h + P_{n-1} Q_n + \cdots + P_1 Q_2 \nonumber \\
       & + \phi_m u_m + \phi_{m-1} u_{m-1} + \cdots + \phi_1 u_1 \nonumber \\
       & - L ( Q_1, Q_2, \ldots, h, \Lambda_1, \ldots, \Lambda_m)
\end{align}


After calculating the secondary constraints
with~\eqref{eq:consistency_relations}

\begin{equation*}
  \phi_{2, i} := \left\{ \phi_{1,i}, H \right\} =
  \frac{\partial \phi_{1,i}}{\partial \Pi_i}
  \frac{\partial H}{\partial \Lambda_i} =
  \frac{\partial L}{\partial \lambda_i} \Big|_{\lambda_i = \Lambda_i}
  \approx 0
\end{equation*}
two cases are possible:

\begin{enumerate}
  \item $\{\phi_{1,i}, \phi_{2,j}\} \not\approx 0$: $\phi_{1,i}$, $\phi_{2,j}$ are
    both second-class constraints thus there are no further constraints
    derivable from consistency relations.
  \item $\{\phi_{1,i}, \phi_{2,j}\} \approx 0$: $\phi_{1,i}$, $\phi_{2,j}$ are
    not second-class constraints then further constraints can be found employing
    the chain~\eqref{eq:consistency_relations} until either \emph{case 1} is
    reached or the new constraint can be expressed using those found previously.
\end{enumerate}

After all constraints have been found and the gauge symmetries are fixed, one
can use these equations to write the auxiliary canonical variables in terms of
the ``usual'' canonical variables
\begin{equation} \label{eq:constraints_in_terms_of_canonical_vars}
  \begin{cases}
    \Lambda_i =\ f_i(Q_1, \ldots, Q_n, P_n) \\
    \Pi_i =\ 0
  \end{cases}
  \qquad i = 1, 2, \ldots, m \\
\end{equation}
and then substitute in the Hamiltonian~\eqref{eq:Ham_constraints}
\begin{equation*}
  H =\ P_n h + P_{n-1} Q_n + \cdots + P_1 Q_2
     - L ( Q_1, Q_2, \ldots, h, f_1, \ldots, f_m)
\end{equation*}

To recap, we started form a phase space of dimension $2n$ (n$Q$ + n$P$), then we
enlarged it by adding $2m$ auxiliary variables (m$\Lambda$ + m$\Pi$) so the
total phase space dimension is $2(n+m)$. Then we associated with every pair of
canonical auxiliary variables two constraints (\emph{case 1}) or more
(\emph{case 2}) have been found. During the last substitution the dimension of
phase space is bring down to it's original dimensions (\emph{case 1}: $2(n+m) -
2m = 2n$) or further reduced (\emph{case 2}: e.g. $2(n+m) - 3m < 2n$).

\subsubsection{Example of non-stable system with constraints}
We introduce here a widely studied HD Lagrangian: the
\emph{Pais-Uhlenbeck oscillator} (PU).
\begin{equation} \label{eq:lagrangian_PU}
  L_{PU} = \frac{1}{2} \left[
    \ddot{q}^2 - (\omega_1^2 + \omega_2^2) \dot{q}^2 + \omega_1^2 \omega_2^2 q^2
  \right]
\end{equation}

According to Ostrogradsky Theorem~\ref{th:ostrogradsky_classical} this system
shows a problematic instability. Then we can try to cure that by introducing a
constraint such as $\ddot{q}^2 - \dot{q}^2 = 0$ and see if the Hamiltonian is
still two-side unbounded. Using the auxiliary variable $\lambda$ one can write
the constraint version of~\eqref{eq:lagrangian_PU}:
\begin{equation} \label{eq:lagrangian_PUC}
  L_{PUC} = L_{PU} + \frac{\lambda}{2} \left(\ddot{q}^2 - \dot{q}^2\right)
\end{equation}
and from it we derive the Hamiltonian following the given prescription.
\begin{align*}
  & \begin{cases}
    Q_1 := q \\
    P_1 := \frac{\delta L_{PUC}}{\delta \dot{q}} =
        - (\lambda + 1) q^{(3)} - \dot{\lambda}\ddot{q}
        - (\lambda + \omega_1^2 + \omega_2^2)\dot{q}
  \end{cases} \\
  & \begin{cases}
    Q_2 := \dot{q} \\
    P_2 := \frac{\delta L_{PUC}}{\delta \ddot{q}} =
        (\lambda + 1) \ddot{q}
        \quad  \Rightarrow \quad \ddot{q} = \frac{P_2}{(\Lambda + 1)}
  \end{cases} \\
  & \begin{cases}
    \Lambda := \lambda \\
    \Pi := \frac{\delta L_{PUC}}{\delta \lambda} = 0
        \quad  \Rightarrow \quad \phi_1: \Pi \approx 0
  \end{cases}
\end{align*}
\begin{equation} \label{eq:ham_PUC_with_aux_vars}
   H_{PUC} =\ P_1 Q_2 + \frac{P_2^2}{2(\Lambda + 1)} + \phi_1 u_1
  + \frac{Q_2^2}{2} \left(\Lambda + \omega_1^2 + \omega_2^2 \right)
  - \frac{Q_1^2}{2} \omega_1^2 \omega_2^2
\end{equation}

\begin{equation} \label{eq:secondary_constraint_PU}
  \phi_2 :=
  \left\{ \phi_1, H_{PUC} \right\} =
  \frac{1}{2} \left[ \frac{P_2}{\Lambda + 1} - Q_2 \right]
  \left[ \frac{P_2}{\Lambda + 1} + Q_2 \right]
  \approx 0
\end{equation}

The secondary constraint equation~\eqref{eq:secondary_constraint_PU}
vanish if one of the factors in square brackets vanish. The two square bracket
represent two hypersurfaces and one solution is chosen instead of the other
based on initial conditions. In the following, we choose the bracket with the
minus sign. Moreover note that $\phi_1$ and $\phi_2$ are already second-class
constraints, indeed
\begin{equation*}
  \begin{cases}
    \phi_1: \Pi = 0 \\
    \phi_2: \frac{P_2}{\Lambda + 1} - Q_2 \approx 0
  \end{cases}
  \quad \Rightarrow \quad
  \left\{ \phi_1, \phi_2 \right\} = -
  \frac{\partial\phi_1}{\partial\Pi} \frac{\partial\phi_2}{\partial\Lambda} =
  \frac{P_2}{{(\Lambda + 1)}^2} \not\approx 0
\end{equation*}
so no further constraints are needed (\emph{case 1}). The equation for
$\phi_1$ (trivial) and for $\phi_2$ can be locally
inverted~\eqref{eq:constraints_in_terms_of_canonical_vars} obtaining that
$\Pi = 0$ and $\Lambda = P_2/Q_2 - 1$. Substituting them
in~\eqref{eq:ham_PUC_with_aux_vars} we get the final Hamiltonian of the system
but unfortunately the Ostrogradsky instability persists.
\begin{equation} \label{eq:ham_PUC_instable}
   H_{PUC} =\ P_1 Q_2 + P_2 Q_2
  + \frac{Q_2^2}{2} \left(\omega_1^2 + \omega_2^2 - 1 \right)
  - \frac{Q_1^2}{2} \omega_1^2 \omega_2^2
\end{equation}


\subsubsection{Example of stable system with constraints}
Now consider the following constrained Pais-Uhlenbeck oscillator~\cite{Chen13}
\begin{equation}
  L_{PUC} = L_{PU} +
  4 \omega_1^2\omega_2^2 q^2 \lambda (1+ \lambda) +
  2 \sqrt{2} \omega_1\omega_2 \lambda q \ddot{q}
\end{equation}

Analogously to the previous example, canonical variable and hamiltonian
$H_{PUC}$ can be derived in similar fashion
\begin{align*}
  & \begin{cases}
    Q_1 := q \\
    P_1 := - q^{(3)} -
    \left( 2\sqrt{2} \omega_1\omega_2\lambda + \omega_1^2 + \omega_1^2 \right)
    \dot{q} - 2\sqrt{2}\omega_1\omega_2\dot{\lambda} q
  \end{cases} \\
  & \begin{cases}
    Q_2 := \dot{q} \\
    P_2 := \ddot{q} + 2\sqrt{2} \omega_1\omega_2\lambda q
    \quad \Rightarrow \quad
    \ddot{q} = P_2 - 2\sqrt{2} \omega_1\omega_2\Lambda Q_1
  \end{cases} \\
  & \begin{cases}
    \Lambda := \lambda \\
    \Pi := 0
    \quad  \Rightarrow \quad \phi_1: \Pi \approx 0
  \end{cases}
\end{align*}
\begin{align} \label{eq:ham_PUC_with_aux_vars_stable}
  H_{PUC} =\ & P_1 Q_2 + \frac{P_2^2}{2} -
  \omega_1^2\omega_2^2\frac{Q_1^2}{2} +
  \left( \omega_1^2 + \omega_2^2 \right) \frac{Q_2^2}{2} \\
  & - 4 \omega_1^2\omega_2^2\Lambda Q_1^2 -
  2\sqrt{2}\omega_1\omega_2\Lambda P_2 Q_1 + \phi_1 u_1
\end{align}

Secondary constraints relations carry useful information up to the forth order
\footnote{
  To obtain $\phi_4$, information from $\phi_1$ and $\phi_2$ have been used in
  order to write $\phi_4$ in terms of $\Lambda$.
}
\begin{align*}
  & \phi_1 : \Pi = 0 \\
  & \phi_2 := \{\phi_1, H_{PUC}\} =
    P_2 + \sqrt{2} \omega_1\omega_2 Q_1 \approx 0 \\
  & \phi_3 := \{\phi_2, H_{PUC}\} =
    P_1 + \left(\omega_1^2 + \omega_2^2 - \sqrt{2}\omega_1\omega_2 \right) Q_2
    \approx 0 \\
  & \phi_4 := \{\phi_3, H_{PUC}\} =
  \omega_1 \omega_2 (3+8\Lambda) - \sqrt{2}
  \left( \omega_1^2 + \omega_2^2 \right) \left(1 + 2\Lambda \right) \approx 0
\end{align*}

These consistency relations can be inverted so that $P_2$, $Q_2$ and $\Lambda$
are express in terms of $P_1$ and $Q_1$. Substituting them
into~\eqref{eq:ham_PUC_with_aux_vars_stable} obtaining

\begin{equation} \label{eq:ham_PUC_stable}
  H_{PUC} =\ \frac{\omega_1^2\omega_2^2}{2} Q_1^2 +
  \frac{\omega_1\omega_2}{\sqrt{2} {\left(\sqrt{2} \omega_1\omega_2 -
  \omega_1^2 - \omega_2^2 \right)}^2} P_1^2
\end{equation}

The Hamiltonian~\eqref{eq:ham_PUC_stable} is bounded from below (and not from
above, just like an harmonic oscillator) and hence this system does not suffer
of Ostrogradsky instability unlike the previous example. \\

The main difference between the two is the role played by constraints. In the
first example the dimensions of phase space of the unconstrained and constrained
Hamiltonian are both four (i.e.\ four variable appear in the reduced
Hamiltonian~\eqref{eq:ham_PUC_instable}: $Q_1$, $P_1$, $Q_2$ and $P_2$). In the
second example the constraints reduce the dimensionality of the phase space from
the original four to two (indeed just $Q_1$ and $P_1$ appear
in~\eqref{eq:ham_PUC_stable}). As shown in~\cite{Chen13} a reduction of the
phase space, due to imposed constraints, can cure the typical Ostrogradsky
instability of HD theories.
