# Extended Chern-Simons model for a vector multiplet

D.S<EMAIL_ADDRESS>$\,{}^{1,2}$, S.L.
<EMAIL_ADDRESS>$\,{}^{1}$, and O.D.
<EMAIL_ADDRESS>$\,{}^{1}$

(1 Physics Faculty, Tomsk State University, Tomsk 634050, Russia
2 _Lebedev Institute of Physics, Leninsky ave. 53, Moscow 119991, Russia_)

###### Abstract

We consider a gauge theory of vector fields in $3d$ Minkowski space. At the
free level, the dynamical variables are subjected to the extended Chern-Simons
(ECS) equations with higher derivatives. If the color index takes $n$ values,
the third-order model admits a $2n$-parameter series of second-rank conserved
tensors, which includes the canonical energy-momentum. Even though the
canonical energy is unbounded, the other representatives in the series can
have bounded from below $00$-component. The theory admits consistent self-
interactions with the Yang-Mills gauge symmetry. The Lagrangian couplings
preserve the unbounded from below energy-momentum tensor, and they do not lead
to a stable non-linear theory. The non-Lagrangian couplings are consistent
with the existence of conserved tensor with a bounded from below
$00$-component. These models are stable at the non-linear level. The dynamics
of interacting theory admits a constraint Hamiltonian form. The Hamiltonian
density is given by the $00$-component of the conserved tensor. In the case of
stable interactions, the Poisson bracket and Hamiltonian do not follow from
the canonical Ostrogradski construction. The particular attention is paid to
the ”triply massless” ECS theory, which demonstrates instability already at
the free level. It is shown that the introduction of extra scalar field,
serving as Higgs, can stabilize dynamics in the vicinity of the local minimum
of energy. The equations of motion of stable model are non-Lagrangian, but
they admit the Hamiltonian form of dynamics with a bounded from below
Hamiltonian.

##  1 Introduction

The higher derivative theories are notorious for better convergency properties
at classical and quantum level, and a wider symmetry. In most of the
instances, these advantages come with a price of dynamics instability, being
the typical problem for the models in this class. At the classical level the
solutions to the equations of motion demonstrate runaway behavior
(”collapse”). At the quantum level the ghost poles appear in the propagator,
so the unitarity of dynamics is an issue. These peculiarities follow from a
single fact: the canonical energy is unbounded in every non-singular higher
derivative theory. For a review of the problem, we cite recent articles [2, 1,
3] and references therein. The canonical energy of singular models can be
bounded. The examples include $f(R)$-theories of gravity [4, 5, 6, 7]. These
models do not demonstrate instability. The stability problem for constrained
higher derivative theories is discussed in [8]. In the majority of interesting
(constrained or unconstrained) models, the instability has the special form:
the classical dynamics is regular (no ”collapse”), but the canonical energy is
unbounded from below. The Pais-Uhlenbeck oscillator [9], Podolsky [10] and
Lee-Wick [11, 12] electrodynamics, ECS theory [13], and conformal gravity [14]
are examples. For the current studies we mention [16, 17, 18, 15] and
references therein. In all these models the quantum theory is expected to be
well-defined [2, 19], but the application of the canonical quantization
procedure based on the Ostrogradski procedure444The canonical Hamiltonian
formalism for the non-singular Higher-derivative theories has been first
proposed in [20]. Its generalization for singular theories has been developed
in [21]: see also [22]. For recent studies see [23, 24]. leads to the model
with an unbounded from below spectrum of energy. That is why the stability and
unitarity of quantum dynamics is the most important issue.

Recently, it has been recognized that the higher derivative dynamics can be
stable at the classical and quantum levels even if the canonical energy of the
model is unbounded. The articles [25, 26] explain the stability of Pais-
Uhlenbeck theory form the viewpoint of existence of the Hamiltonian form of
dynamics with a bounded (non-canonical) Hamiltonian. The quantization of the
model with an alternative Hamiltonian leads to quantum theory with a well-
defined vacuum state with the lowest energy. The papers [27, 28] use a special
PT-symmetry for construction of stable quantum mechanics in the higher
derivative oscillator model. The articles [2, 1, 3, 29] tell us that the non-
linear higher derivative models may have a well-defined classical dynamics
without ”collapsing” trajectories with the runaway behavior. In the last
articles, the existence of stable classical dynamics serves as a necessary
prerequisite for construction of a well-defined quantum theory. For the
studies of field theories, we refer to articles [30, 31]. All the mentioned in
this paragraph models have one common feature: they admit an alternative
(different from the canonical energy) bounded from below conserved quantities.
The quantum stability is achieved if the model admits Hamiltonian form of
dynamics, with the bounded conserved quantity being the Hamiltonian. This
means that the stable higher derivative theory is characterized by two
ingredients: the bounded conserved quantity, and the Poisson bracket on its
phase space that brings the equations of motion to the Hamiltonian form with a
bounded from below Hamiltonian.

In the article [32], the stability is studied in the class of models of
derived type. At the free level, the wave operator of the theory is given by a
polynomial (characteristic polynomial) in another formally self-adjoint
operator of lower order. It has been shown that each derived theory admits a
series of conserved tensors, which includes the canonical energy-momentum
[33]. The number of entries in the series grows with the order of
characteristic polynomial. Even though the canonical energy is unbounded, the
other conserved tensors in the series can be bounded from below [33]. The
bounded conserved quantity stabilizes the classical dynamics of the theory.
The quantum stability is explained by the existence of the Lagrange anchor
that relates a bounded conserved quantity with the invariance of the model
with respect to the time translations.555The concept of the Lagrange anchor
has been proposed in [34] in the context of quantization of non-Lagrangian
field theories. Later, it has been recognized that the Lagrange anchor
connects symmetries and conserved quantities [35]. In the first-order
formalism, the Lagrange anchor defines the Poisson bracket [36]. The
Hamiltonian is given by the conserved quantity in the series, being expressed
in terms of phase-space variables. The linear higher derivative theory is
stable if a bounded conserved quantity, being related to the time-translation
symmetry, is admitted by the model [32]. The articles [37, 38] consider a
problem of consistent deformations of symmetries and conserved quantities that
preserve the stability of dynamics in the class of derived type models. In all
the cases interaction vertices do not come from the least action principle
with higher derivatives, but the equations of motion admit the Hamiltonian
form of dynamics. The main difficulty of the cited above procedures is that
they do not automatically preserve gauge invariance. This restricts their
applications in gauge systems.

The ECS model [13] is the simplest gauge theory of derived type. In the
current studies, it is often serves as a prototype of the class of gauge
theories with higher derivatives, see eg. [39, 40, 41]. The stability of the
ECS model has been first studied in [33]. It has been observed that the theory
of order $p$ admits an $p-1$-parameter series of second-rank conserved
tensors, which includes the canonical energy-momentum. The canonical energy of
the model is always unbounded from below. The other tensors in the series may
have a bounded from below 00-component. The stability conditions for the
theory are determined by the structure of the roots of characteristic
polynomial [33]. The theory is stable if all the nonzero roots of
characteristic polynomial are real and different and zero root has
multiplicity one or two. The stability of the ECS theory is confirmed in work
[42]. The ECS model has been shown to be a multi-Hamiltonian at the free level
in [43]. The series of Hamiltonians includes the canonical (Ostrogradski)
Hamiltonian, which is unbounded, and the other representatives, being bounded
or unbounded depending on the model parameters. If the Hamiltonian is bounded,
it ensures the stability of the model at both classical and quantum levels.
The model admits inclusion of stable non-Lagrangian interactions with scalar,
fermionic and gravitational fields that preserve a selected representative in
the series of conserved quantities of free model [43, 44, 38]. However, the
gauge symmetry is abelian in the sector of vector field in all these examples.

The concept of consistency of interaction between the gauge fields has been
developed in [45]. The consistent couplings between gauge vectors has been
studied in [46, 47]. It has been shown that the Yang-Mills vertex is unique in
the covariant setting. The Lagrangian self-interactions of the gauge vector
multiplet subjected to the ECS equations at free level has been reconsidered
in article [42]. The most general consistent non-linear theory is proven to
have the Yang-Mills gauge symmetry, and the Lagrangian is given by the
covariantization of the free ECS Lagrangian. The interacting model
demonstrates the Ostrogradski instability at the non-linear level in all the
instances, even though the free theory is stable. This result implies a no-go
theorem for stable Lagrangian interactions in the ECS model. The conclusion is
not surprising because the Lagrangian couplings preserve the conservation law
of canonical energy, being unbounded already at the free level. The inclusion
of non-Lagrangian interactions can solve the issue of the dynamics stability
at the interacting level, because such couplings preserve bounded conserved
quantities. However, the problem of construction of stable consistent non-
Lagrangian (self-)couplings (stable or unstable) has not been studied in the
ECS model in the literature before.

In the current work, we present a class of stable self-interactions in the
theory of vector multiplet. The free fields are subjected to the ECS equations
of third order. The interaction is (in general) non-Lagrangian. The non-linear
equations of motion are consistent with the Yang-Mills gauge symmetry. A
selected second-rank conserved tensor of the free model is preserved by the
coupling. Depending on the values of coupling constants, it can be bounded or
unbounded. The Lagrangian interaction vertex of the article [42] is unstable.
The non-Lagrangian coupling can be consistent with the existence of bounded
conserved tensor. The equations of motion admit the Hamiltonian form of
dynamics. On shell, the Hamiltonian density is given by the $00$-component of
the conserved tensor. For Lagrangian interactions, the canonical formalism
with the unbounded Ostrogradski Hamiltonian is reproduced. The bounded
Hamiltonians do not follow from the Ostrograski procedure, so we have non-
canonical Hamiltonian formalism in this case. In all the instances, the
Poisson bracket is a non-degenerate tensor, so the model admits a Hamiltonian
action principle. The quantization of the first-order action with the bounded
Hamiltonian leads to a stable quantum theory with a well-defined vacuum state.

In the case of resonance (multiple roots of charlatanistic polynomial), the
dynamics of the theory is unstable already at the free level. The inclusion of
self-interaction does not improve stability properties of dynamics of this
model because the conserved quantities have one and the same structure in the
free and non-linear models. In the current article, a model with the third-
order zero resonance root is of interest. The free field is subjected to the
”triply massless” Chern-Simons (CS) equations. The wave operator of free
equations is given by a cube of the CS operator. To stabilize the dynamics at
the non-linear level, we apply the ”Higgs-like” mechanism described in [48].
Introducing an extra scalar, we generate a coupling such that the energy of
the interacting theory gets a local minimum for a solution with nonzero values
of dynamical variables. The motions of small fluctuations in the vicinity of
this solution are stable because the non-linear theory has no resonance. The
dynamics of fluctuations admits the Hamiltonian form, with the Hamiltonian
being given by the positive definite quadratic form of the fields. This means
that the dynamics of the theory with a resonance can be stabilized by the
inclusion of appropriate interaction.

The article is organized as the following. In Section 2, we consider the
third-order ECS model for a vector multiplet. The particular attention is paid
to the structure of symmetries and conserved quantities of the theory, and the
stability of dynamics. In Section 3, we propose stable self-interactions
between the multiplet of vector fields with the Yang-Mills gauge symmetry. The
general non-model in the class model preserve a selected conserved tensor of
the free theory, which can be bounded or unbounded from below depending on the
values of coupling constants. The stable interactions are inevitably non-
Lagrangian. In Section 4, we construct the constrained Hamiltonian formalism
for the non-linear model. The density of Hamiltonian is given by the
00-component of the conserved tensor, being expressed in terms of the phase-
space variables. In Section 5, we consider the issue of stability of ”triply
massless” ECS theory. We propose the class of consistent couplings with the
scalar field stabilizing the dynamics in the vicinity of equilibrium position.
The equations of motion are non-Lagrangian, but they admit the constrained
Hamiltonian form with a bounded Hamiltonian. The concluding section summarizes
the results of the article.

##  2 Higher-derivative Chern-Simons model

We consider the most general third-order gauge theory of the vector multiplet
$A{}^{a}=A{}^{a}{}_{\mu}(x)dx^{\mu}$, $\mu=0,1,2$, $a=1,2,\ldots,n$ on three-
dimensional Minkowski space. The action functional of the model reads

$S[A(x)]=\frac{1}{2}\int
A^{a}{}_{\mu}(\alpha_{1}F^{a\mu}+\alpha_{2}G^{a\mu}+\alpha_{3}K^{a\mu})d^{3}x\,.$
(1)

The real numbers $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ are model
parameters. Without loss of generality, we assume that the coefficient at the
higher derivative term is nonzero, $\alpha_{3}\neq 0$. The vectors
$F{}^{a}{}_{\mu}$, $G{}^{a}{}_{\mu}$, $K{}^{a}{}_{\mu}$ denote the
(generalized) field strengthes of the potential $A{}^{a}{}_{\mu}$,

$F{}^{a}{}_{\mu}=\varepsilon_{\mu\nu\rho}\partial^{\nu}A{}^{a\rho}\,,\qquad
G{}^{a}{}_{\mu}=\varepsilon_{\mu\nu\rho}\partial^{\nu}F{}^{a\rho}\,,\qquad
K{}^{a}{}_{\mu}=\varepsilon_{\mu\nu\rho}\partial^{\nu}G{}^{a\rho}\,,$ (2)

with the Levi-Civita symbol $\epsilon_{\mu\nu\rho},\epsilon_{012}=1$ being
antisymmetric with respect to indices. All the tensor indices are raised and
lowered by the Minkowski metric. We use the mostly positive convention for the
signature of metric throughout the paper. The summation is implied over
repeated at one level isotopic indices $a=1,\ldots,n$ unless otherwise stated.
For $\alpha_{3}=\alpha_{1}=1$, $\alpha_{2}=0$, $n=1$, the action functional
(1) has been first proposed in the article [13]. We call the model (1) the ECS
theory for a vector multiplet.

The least action principle for the functional (1) brings us to the following
Lagrange equations for the field $A_{\mu}$:

$\frac{\delta S}{\delta
A^{a\mu}}\equiv\alpha_{1}F{}^{a}{}_{\mu}+\alpha_{2}G{}^{a}{}_{\mu}+\alpha_{3}K{}^{a}{}_{\mu}=0\,,\qquad
a=1,\ldots,n\,.$ (3)

These equations have the derived form [33] because the wave operator is a
polynomial in the CS operator $\ast d$. The symbol $\ast$ denotes the Hodge
star operator, and $d$ is the de Rham differential. The structure of the
Poincare group representation, being described by the equation (3), is
determined by the roots of the characteristic polynomial

$M(\alpha;z)=\alpha_{1}z+\alpha_{2}z^{2}+\alpha_{3}z^{3}\,.$ (4)

The polynomial $M(\alpha;z)$ follows from (3) after the formal replacement of
the CS operator $\ast d$ by the complex-valued variable $z$. The following
cases are distinguished in [33].

(1a) $\alpha_{1}\neq 0\,,\alpha{}_{2}{}^{2}-4\alpha_{1}\alpha_{3}>0\,$. The
characteristic polynomial has a zero root, and two different nonzero real
roots. Equation (3) describes two massive spin-1 self-dual fields.666A theory
of massive spin-1 field subjected to the self-duality condition has been
proposed in [49]. For the theory of representations of $3d$ Poincare group, we
cite [50, 51, 52]. A zero root corresponds to the CS mode, being pure gauge.

(1b) $\alpha_{1}\neq 0\,,\alpha_{2}^{2}-4\alpha_{1}\alpha_{3}<0$. The
characteristic polynomial has a zero root and two complex roots. The case is
similar to (1a), but the masses of vector modes are complex. The Poincare
group representation is non-unitary.

(2a) $\alpha_{2}\neq 0,\alpha_{1}=0$. The characteristic polynomial has a
multiplicity two zero root and simple nonzero root. The set of
subrepresentations includes a massless spin-1 field and a massive spin-1 mode
subjected to the self-duality condition.

(2b) $\alpha_{2}\neq 0,\alpha_{2}^{2}-4\alpha_{1}\alpha_{3}=0$. The
characteristic polynomial has a multiplicity two nonzero root, and a simple
zero root. The subrepresentations describe ”doubly massive” massive mode, and
a gauge CS mode. The Poincare group representation is non-unitary.

(3) $\alpha_{1}=\alpha_{2}=0$. The characteristic polynomial has a
multiplicity three zero root. Equation (3) describes the ”triply” massless
extended CS theory. The Poincare group representation is non-unitary.

As we see, the ECS model describes the unitary representation of the Poincare
group if all the nonzero roots of characteristic equation are different, and
the zero root has multiplicity one or two. The dynamical degrees of freedom
include the massive spin 1 vector subjected to the self-duality condition
and/or spin 1 massless field, which meets the $3d$ Maxwell equations. In all
the instances, the model has two local physical degrees of freedom (four
physical polarizations).

The action functional (1) is preserved by the $2n$-parameter series of
infinitesimal transformations

$\delta_{\xi;\beta}A{}^{a}{}_{\mu}=-\varepsilon_{\mu\nu\rho}\xi{}^{\rho}(\beta{}^{a}{}_{1}F{}^{a}{}^{\rho}+\beta{}^{a}{}_{2}G{}^{a}{}^{\rho})\,,\qquad
a=1,\ldots,n,$ (5)

(no summation in $a$). The transformation parameters are the constant vector
$\xi{}^{\mu}$ and real numbers $\beta{}^{a}{}_{k}$, $a=0,\ldots,n$, $k=1,2$.
The series (5) includes the space-time translations with the independent
parameters for the individual vector $A{}^{a}{}_{\mu}$ in the multiplet, and
higher order transformations, whose value is determined by
$\beta{}^{a}{}_{2}$. The Noether theorem associates symmetries (5) with the
$2n$-parameter series of second-rank conserved tensors, which has the form

$\begin{array}[]{c}\displaystyle\Theta^{\mu\nu}(\beta;\alpha)=\sum_{a=1}^{n}(\beta{}^{a}{}_{1}\Theta{}^{a\mu\nu}{}_{1}(\alpha)+\beta{}^{a}{}_{2}\Theta{}^{a\mu\nu}{}_{2}(\alpha))\,,\end{array}$
(6)

where

$\begin{array}[]{c}\displaystyle\Theta{}^{a\mu\nu}{}_{1}(\alpha)=\alpha_{3}(G{}^{a\mu}F{}^{a\nu}+G{}^{a\nu}F{}^{a\mu}-g^{\mu\nu}G{}^{a}{}_{\rho}F{}^{a\rho})+\alpha_{2}(F{}^{a\mu}F{}^{a\nu}-\frac{1}{2}g^{\mu\nu}F{}^{a}{}_{\rho}F{}^{a\rho})\,;\end{array}$
(7)
$\begin{array}[]{c}\displaystyle\Theta{}^{a\mu\nu}{}_{2}(\alpha)=\alpha_{3}(G{}^{a\mu}G^{a\nu}-\frac{1}{2}g^{\mu\nu}G{}^{a}{}_{\rho}G{}^{a\rho})-\alpha_{1}(F{}^{a\mu}F{}^{a\nu}-\frac{1}{2}g^{\mu\nu}F{}^{a}{}_{\rho}F{}^{a\rho})\end{array}$
(8)

(no summation in $a$). The quantities $\Theta{}^{a\mu\nu}{}_{1}$,
$a=1,\ldots,n$ (7) represent the canonical energy-momentum tensors of
individual fields in the vector multiplet $A{}^{a}{}_{\mu}$. The tensors
$\Theta{}^{a\mu\nu}{}_{2}$, $a=1,\ldots,n$ (8) are other conserved quantities.
The total number of independent conserved tensors in the free theory is $2n$
because each field in the multiplet admits two different symmetries.

The 00-component of the conserved tensor (6) reads

$\begin{array}[]{c}\displaystyle\Theta^{00}(\beta;\alpha)=\frac{1}{2}\sum_{a=1}^{n}\Big{\\{}\beta_{2}\alpha_{3}(G{}^{a0}G{}^{a0}+G{}^{ai}G{}^{ai})+2\beta_{1}\alpha_{3}(G{}^{a0}F{}^{a0}+G{}^{ai}F{}^{ai})+\\\\[14.22636pt]
\displaystyle+(\beta_{1}\alpha_{2}-\beta_{2}\alpha_{1})(F{}^{a0}F{}^{a0}+F{}^{ai}F{}^{ai})\Big{\\}}\,.\end{array}$
(9)

The summation over repeated at one level index $i=1,2$ is implied. The
quadratic form (9) can be reduced to the principal axes as follows

$\begin{array}[]{c}\displaystyle\Theta^{00}(\beta;\alpha)=\sum_{a=1}^{n}\Big{\\{}\frac{\alpha_{3}(X{}^{a0}X{}^{a0}+X{}^{ai}X{}^{ai})}{2\beta{}^{a}{}_{2}}+\frac{C(\beta^{a};\alpha)(F{}^{a0}F{}^{a0}+F{}^{ai}F{}^{ai})}{2\beta{}^{a}{}_{2}}\Big{\\}}\,.\end{array}$
(10)

Here,
$X{}^{a}_{\mu}=\beta{}^{a}{}_{1}F{}^{a}{}_{\mu}+\beta{}^{a}{}_{2}G{}_{\mu}$,
and the notation is used

$C(\beta{}^{a};\alpha)=-(\beta{}^{a}{}_{2}){}^{2}\alpha{}_{1}+\beta{}^{a}{}_{2}\beta{}^{a}{}_{1}\alpha_{2}-(\beta{}^{a}{}_{1}){}^{2}\alpha_{3}\,.$
(11)

In Section 4, we see that $X{}^{a0}X{}^{a0}+X{}^{ai}X{}^{ai}$ and
$F{}^{a0}F{}^{a0}+F{}^{ai}F{}^{ai}$ depend on different initial data.777See
the detailed discussion in [43]. In this case, the $00$-component (9) is
bounded from below if the coefficients at squares are positive,

$\phantom{\frac{1}{2}}\beta{}^{a}{}_{2}\alpha{}_{3}>0,\qquad
C(\beta{}^{a};\alpha)>0\,,\qquad a=1,\ldots,n\,.\phantom{\frac{1}{2}}$ (12)

We ignore the case of semi-positive quadratic form because the degenerate
conserved quantities do not ensure the stability of dynamics. Relations (12)
are consistent if and only if the model parameters
$\alpha_{1},\alpha_{2},\alpha_{3}$ meet conditions (1a), (2a) of
classification on pages 4-5. In the cases (1b), (2b), (3) of this
classification, conditions (12) are inconsistent. This means that the free
model (1) is stable if the wave equation (3) describes unitary representation
of the Poincare group. This result is quite natural because the theories,
which correspond to the unitary representations of the Poincare group, should
have a stable dynamics.

##  3 Stable interactions

In this section, we present an example of stable self-interactions in the
model (1). The interactions are non-Lagrangian. The dynamics of the non-linear
theory is determined by the equations of motion. The interactions are
associated with the deformations of free equations of motion that preserve the
gauge symmetries and gauge identities of free model. The interaction is
consistent if the number of physical degrees of freedom is preserved by
coupling. For details of the concept of consistency of interaction in the
class of not necessary Lagrangian theories, we refer article [53].

We start construction of non-linear theory by assuming that the dynamical
fields take values in the Lie algebra of a semisimple Lie group with the
generators $t{}^{a},a=1,\ldots,n$,

$\mathcal{A}{}_{\mu}=A{}^{a}{}_{\mu}(x)t^{a}dx^{\mu}\,,\qquad[t{}^{a},t{}^{b}]=f^{abc}t^{c}\,,\qquad\text{tr}(t{}^{a}t{}^{b})=\delta^{ab}\,.$
(13)

The covariant analogs of the (generalized) field strength vectors (2) are
defined as follows:

$\mathcal{F}{}_{\mu}=\epsilon_{\mu\nu\rho}(\partial^{\nu}\mathcal{A}{}^{\rho}+\frac{1}{2}[\mathcal{A}{}^{\nu},\mathcal{A}{}^{\rho}])\,,\qquad\mathcal{G}_{\mu}=\epsilon_{\mu\nu\rho}D^{\nu}\mathcal{F}{}^{\rho}\,,\qquad\mathcal{K}_{\mu}=\epsilon_{\mu\nu\rho}D^{\nu}\mathcal{G}{}^{\rho}.$
(14)

The vectors $\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$,
$\mathcal{K}{}_{\mu}$ lie in the Lie algebra of a semisimple Lie group, and
$D$ stands for the covariant derivative

$D_{\mu}=\partial_{\mu}+[\mathcal{A}_{\mu},\,\cdot\,]\,.$ (15)

The vector $\mathcal{F}{}_{\mu}$ represents a dual of the standard Yang-Mills
strength tensor. The generalized field strengthes $\mathcal{G}{}_{\mu}$,
$\mathcal{K}{}_{\mu}$ are other covariant quantities that involve second and
third derivatives of $\mathcal{A}{}_{\mu}$.

We consider the interactions that are polynomial in the gauge invariants (14)
and does not involve the highest derivative in the non-linear part. In this
setting, the most general self-consistent non-linear theory is determined by
the following equations of motion:

$\mathbb{T}{}^{\mu}=\alpha{}_{3}\mathcal{K}{}^{\mu}+\alpha{}_{2}\mathcal{G}{}^{\mu}+\alpha{}_{1}\mathcal{F}{}^{\mu}-\frac{\alpha{}_{3}{}^{2}}{2C(\beta;\alpha)}\epsilon^{\mu\nu\rho}[\beta{}_{1}\mathcal{F}{}_{\nu}+\beta{}_{2}\mathcal{G}{}_{\nu},\beta{}_{1}\mathcal{F}{}_{\rho}+\beta{}_{2}\mathcal{G}{}_{\rho}]=0\,.$
(16)

We prove the uniqueness of this coupling in Appendix A. The model parameters
are the real numbers $\alpha_{k},k=1,2,3$, $\beta_{l}$, $l=1,2$, and
$f^{abc},a,b,c=1,\ldots,n$. The constants $\alpha{}_{1},\alpha{}_{2}$
determine the free limit of the equations (16). The numbers $\beta{}_{1}$,
$\beta{}_{2}$ distinguish admissible couplings with one and the same gauge
group. Throughout this section and below, we assume that $C(\beta;\alpha)\neq
0$. The equations of motion (16) come from the least action principle if
$\beta_{1}=1,\beta_{2}=0$. In this case, the action functional reads

$S[\mathcal{A}(x)]=\frac{1}{2}\text{tr}\int\mathcal{F}{}_{\mu}(\alpha_{1}\mathcal{A}{}_{\mu}+\alpha_{2}\mathcal{F}{}^{\mu}+\alpha_{3}\mathcal{G}{}^{\mu})d^{3}x\,.$
(17)

This action functional was first derived in [42].888The higher-derivative
Yang-Mills theory having a similar structure of the Lagrangian has been known
long before [54]. The same paper tells us that the (17) is the most general
form of consistent self-interactions in the gauge theory of vector fields.
This means that the most general consistent Lagrangian coupling (17) is
included into the model (16). If the parameter $\beta_{2}$ is nonzero,
equations (16) do not follow from the least action principle for any
functional with higher derivatives. The variational principle with auxiliary
fields still exists, even if the higher derivative model is non-Lagrangian. In
the last case, the theory (16) admits consequent quantization, and
establishing relationship between symmetries and conserved quantities.

The concept of interaction consistency for not necessarily Lagrangian theories
has been developed in [53]. This paper tells us that the non-Lagrangian
interaction is consistent if the non-linear theory admits the same number of
(i) gauge symmetries, (ii) gauge identities, and (iii) physical degrees of
freedom as a free model. All these fact are easily verified. At first, the
equations of motion (16) are preserved by the Yang-Mills gauge symmetry,

$\phantom{\frac{1}{2}}\delta_{\zeta}\mathcal{A}_{\mu}=D_{\mu}\zeta\,,\qquad\delta_{\zeta}\mathbb{T}_{\mu}=[\,\zeta\,,\mathbb{T}_{\mu}]\,,\phantom{\frac{1}{2}}$
(18)

with $\zeta=(\zeta{}^{a}(x),a=1,\ldots,n)$ being the gauge transformation
parameter. The free model (17) is preserved by the standard gradient gauge
symmetry, $\delta_{\zeta}A{}^{a}{}_{\mu}=\partial_{\mu}\zeta^{a}$. As it is
required, the gauge symmetry (18) is given by the deformation of the gradient
gauge symmetry of free model. The important difference is that the gauge
symmetry (18) is non-abelian. So, the inclusion of interaction leans to the
model with non-abelian gauge symmetry. We have an obvious set of gauge
identities between the equations (16),

$\phantom{\frac{1}{2}}\mathbb{D}{}_{\mu}\mathbb{T}{}^{\mu}=0\,,\qquad\mathbb{D}_{\mu}=D_{\mu}+\frac{\alpha{}_{3}}{C(\beta;\alpha)}[\beta{}_{1}\mathcal{F}{}_{\mu}+\beta{}_{2}\mathcal{G}{}_{\mu},\cdot]\,.\phantom{\frac{1}{2}}$
(19)

Again, the leading term of the gauge identity is given by the free
contribution. This agrees with the concept of interaction consistency. At the
final step of analysis, we verify that the physical degrees of freedom number
is preserved by the coupling. Equation (8) of paper [53] expresses the
physical degrees of freedom number via the orders gauge symmetries, gauge
identities, and equations of motion in the involutive of dynamics. The systems
(3), (16) are involutive, and they have equal orders of gauge symmetries,
gauge identities and equations of motion. So, they have to poses one and the
same physical degrees of freedom number. All the above implies that the non-
Lagrangian interaction (16) is consistent for the general values of the
parameters $\beta,\alpha$.

The theory (16) admits a symmetric conserved tensor of second rank of the
following form:

$\begin{array}[]{c}\displaystyle\boldsymbol{\Theta}{}^{\mu\nu}{}(\beta;\alpha)=\text{tr}\Big{\\{}\beta_{2}\alpha_{3}(\mathcal{G}{}^{\mu}\mathcal{G}^{\nu}-\frac{1}{2}g^{\mu\nu}\mathcal{G}{}_{\rho}\mathcal{G}{}^{\rho})+\beta_{1}\alpha_{3}(\mathcal{G}{}^{\mu}\mathcal{F}{}^{\nu}+\mathcal{G}{}^{\nu}\mathcal{F}{}^{\mu}-\\\\[14.22636pt]
\displaystyle-g^{\mu\nu}\mathcal{G}{}_{\rho}\mathcal{F}{}^{\rho})+(\beta_{1}\alpha_{2}-\beta_{2}\alpha_{1})(\mathcal{F}{}^{\mu}\mathcal{F}{}^{\nu}-\frac{1}{2}g^{\mu\nu}\mathcal{F}{}_{\rho}\mathcal{F}{}^{\rho})\Big{\\}}\,,\end{array}$
(20)

with $\alpha_{k},k=1,2,3$ and $\beta_{l},l=1,2$ being the model parameters.
The divergence of the quantity
$\boldsymbol{\Theta}{}^{\mu\nu}{}(\beta;\alpha)$ reads

$\partial_{\nu}\boldsymbol{\Theta}{}^{\mu\nu}{}=-\text{tr}(\epsilon^{\mu\nu\rho}(\beta{}_{1}\mathcal{F}{}_{\nu}+\beta{}_{2}\mathcal{G}{}_{\nu})\mathbb{T}_{\rho})\,.$
(21)

Expression (20) is the covariantization of a selected representative in the
conserved tensor series (6),

$\phantom{\frac{1}{2}}\beta{}^{a}{}_{1}=\beta_{1}\,,\qquad\beta{}^{a}{}_{2}=\beta_{2}\,,\qquad
a=1,\ldots,n\,.\phantom{\frac{1}{2}}$ (22)

As is seen, the model (16) represents the class of deformations of free ECS
equations (3) that preserves a selected representative in the series (6)
conserved quantities at the non-linear level. It is important to note that the
other representatives in the series (20) no longer conserve in the non-linear
theory (16). This happens because the parameters $\beta_{k}$, $k=1,2$ in the
conserved tensor (20) are unambiguously fixed by the interaction.

As far as stability is concerned, the $00$-component of tensor (20) is
relevant. The latter reads

$\begin{array}[]{c}\displaystyle\boldsymbol{\Theta}{}^{00}{}(\beta;\alpha)=\text{tr}\Big{\\{}\frac{1}{2}\beta_{2}\alpha_{3}(\mathcal{G}{}^{0}\mathcal{G}^{0}+\mathcal{G}{}^{i}\mathcal{G}{}^{i})+\beta_{1}\alpha_{3}(\mathcal{G}{}^{0}\mathcal{F}{}^{0}+\mathcal{G}{}^{i}\mathcal{F}{}^{i})+\\\\[14.22636pt]
\displaystyle+\frac{1}{2}(\beta_{1}\alpha_{2}-\beta_{2}\alpha_{1})(\mathcal{F}{}^{0}\mathcal{F}{}^{0}+\mathcal{F}{}^{i}\mathcal{F}{}^{i})\Big{\\}}\,.\end{array}$
(23)

The conserved tensor is a bilinear form in (generalized) strengthes
$\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$. The bringing of the quadratic
form (23) to the principal axes reveals that the model is stable if

$\beta_{2}\alpha_{3}>0\,,\qquad C(\beta;\alpha)>0\,.$ (24)

These conditions can be consistent or inconsistent depending on the values of
the model parameters $\alpha_{k}$, $k=1,2,3$ and $\beta_{l}$, $l=1,2$. The
Lagrangian interaction vertex (17) does not meet stability requirements. This
confirms the instability of the variational coupling proposed in [42]. The
stable interactions in the class of models (16) are inevitably non-Lagrangian.
The similar form of stability conditions (3), (16) at the free and interacting
cases implies that the linear and non-linear dynamics is stable or unstable
simultaneously. In the class of stable at the linear level theories, formula
(16) determines a class of non-linear models that preserve a selected bounded
conserved quantity in the series (6).

Now, we can return to the special case $C(\beta;\alpha)=0$, which is excluded
in (16). The conserved quantity (6), (22) is a degenerate quadratic form of
the initial data. The 00-component of the free conserved tensor can be bounded
from below, but its degeneracy implies existence of zero vector(s). The motion
of the system in the degenerate direction can be infinite, and the
corresponding conserved quantity appears to be irrelevant to stability already
at the free level. The formula (16) prevents construction of interacting
theories that preserve the conserved tensors with a semi-definite
00-component.

##  4 Hamiltonian formalism

In this section, we construct constrained Hamiltonian formalism for the
higher-derivative theory (16).

Let us first explain what do we understand by the constrained Hamiltonian
formalism for the system of not necessarily Lagrangian field equations with
higher derivatives. A general fact is that the higher-derivative system can be
reduced to the first order by introduction of extra fields absorbing the time
derivatives of original dynamical variables. The set of original dynamical
variables and extra fields is denoted by
$\\{\varphi^{I}(\boldsymbol{\mathit{x}},t),\lambda^{A}(\boldsymbol{\mathit{x}},t)\\}$,
where $t=x^{0}$ is the time, and $\boldsymbol{\mathit{x}}=(x^{1},x^{2})$
stands for space coordinates. The multi-indices $A,I$ label the phase-space
variables. The first-order equations are said to be Hamiltonian if there
exists a Hamiltonian function
$\mathcal{H}(\varphi^{I},\boldsymbol{\nabla}\varphi^{I},\boldsymbol{\nabla}^{2}\varphi^{I},...)$
($\boldsymbol{\nabla}$ stands for derivatives by the space
$\boldsymbol{\mathit{x}}$) and a Poisson bracket
$\\{\varphi^{I}(\boldsymbol{\mathit{x}}),\varphi^{J}(\boldsymbol{\mathit{y}})\\}$
such that the equations constitute a constrained Hamiltonian system (the dot
denotes the derivative by time), i.e.

$\displaystyle\begin{array}[]{c}\dot{\varphi}^{I}(\boldsymbol{\mathit{x}})=\\{\varphi^{I}(\boldsymbol{\mathit{x}}),\int\mathcal{H}d\boldsymbol{y}\,\\},\quad\theta_{A}(\varphi^{I}(\boldsymbol{\mathit{x}}),\boldsymbol{\nabla}\varphi^{I}(\boldsymbol{\mathit{x}}),\ldots)=0\,;\\\\[14.22636pt]
\displaystyle\mathcal{H}=\mathcal{H}{}_{0}(\varphi^{J}(\boldsymbol{\mathit{x}}),\boldsymbol{\nabla}\varphi^{J}(\boldsymbol{\mathit{x}}),\ldots)+\lambda^{A}(\boldsymbol{x})\theta_{A}(\varphi^{I}(\boldsymbol{\mathit{x}}),\boldsymbol{\nabla}\varphi^{I}(\boldsymbol{\mathit{x}}),\ldots)\,.\end{array}$
(25)

The system is multi-Hamiltonian if there exists a series of Hamiltonian
$\mathcal{H}_{\beta}$ and Poisson brackets
$\\{\varphi^{I}(\boldsymbol{\mathit{x}}),\varphi^{J}(\boldsymbol{\mathit{y}})\\}_{\beta}$
parameterized by constrains $\beta_{1},\ldots,\beta_{k}$ such that determine
one and the same equations of motion for the dynamical fields. The existence
of the constraint Hamiltonian formulation is not guaranteed for a system of
general non-Lagrangian equations of motion. In particular, the Hamiltonian
formulation for the ECS theory may not exist, at least for certain
combinations of the model parameters.

We start the construction of the Hamiltonian formalism from the reduction of
order of equations (16). The space components of the field strength
$\mathcal{F}_{i},i=1,2,$ and generalized field strength of second order
$\mathcal{G}_{i},i=1,2,$ are chosen as extra fields. By construction, they
absorb the first and second time derivatives of space components of the vector
field $\mathcal{A}_{i},i=1,2,$

$\displaystyle\begin{array}[]{c}\mathcal{F}{}_{i}=\epsilon_{ij}(\dot{\mathcal{A}}{}_{j}-\partial_{i}\mathcal{A}_{0}+[\mathcal{A}_{0},\mathcal{A}_{j}])\,,\\\\[14.22636pt]
\displaystyle\mathcal{G}{}_{i}=-\ddot{\mathcal{A}}{}_{i}+\partial_{i}\dot{A}_{0}+\partial_{j}(\partial_{j}\mathcal{A}_{i}-\partial_{i}\mathcal{A}_{j})+\epsilon_{ij}[\epsilon_{jk}\dot{\mathcal{A}}{}_{0}+\delta{}_{jk}(\epsilon_{rs}\partial{}_{r}\mathcal{A}{}_{s}+\\\\[14.22636pt]
\displaystyle+1/2[\mathcal{A}{}_{r},\mathcal{A}{}_{s}],\mathcal{A}{}_{i}]+\epsilon_{ij}[\mathcal{A}{}_{0},\epsilon{}_{jk}\dot{\mathcal{A}}{}_{k}+\delta_{jk}\mathcal{F}{}_{k}]\,.\end{array}$
(26)

Here, $\epsilon_{ij},\epsilon_{12}=1$ is the $2d$ Levi-Civita symbol. The
Latin indices $i,j$ run over the values $1,2$. The summation over repeated at
one and the same level Latin indices is implied. As is seen form equations,
(26) the fields $\mathcal{F}{}_{i}$ absorb the first time derivatives of space
components of original vector field $\mathcal{A}{}_{i}$. The vector
$\mathcal{G}{}_{i}$ involves second time derivatives of $\mathcal{A}{}_{i}$.
The time components $\mathcal{F}_{0},\mathcal{G}_{0}$ of vectors (26) are
functions of $\mathcal{A}{}_{0}$, $\mathcal{A}{}_{i}$, $\mathcal{F}{}_{i}$,
$\mathcal{G}{}_{i}$, with no new combinations of time derivatives being
involved,

$\displaystyle\begin{array}[]{c}\phantom{\frac{1}{2}}\mathcal{F}{}_{0}\equiv\epsilon_{ij}(\partial{}_{i}\mathcal{A}{}_{j}+[\mathcal{A}_{i},\mathcal{A}_{j}])\,,\qquad\mathcal{G}{}_{0}\equiv\epsilon_{ij}(\partial_{i}\mathcal{F}_{j}+[\mathcal{A}_{i},\mathcal{G}_{j}])\,.\phantom{\frac{1}{2}}\end{array}$
(27)

In remaining part of the article, we associate the quantities
$\mathcal{F}_{0}$, $\mathcal{G}_{0}$ with their expressions in terms of the
phase-space variables $\mathcal{A}{}_{0}$, $\mathcal{A}{}_{i}$,
$\mathcal{F}{}_{i}$.

Substituting the extra variables (26) into the system (16), we obtain the
first-order equations for the fields $\mathcal{A}_{i}$, $\mathcal{F}_{i}$,
$\mathcal{G}_{i}$:

$\phantom{\frac{1}{2}}\dot{\mathcal{A}}{}_{i}=\partial_{i}\mathcal{A}_{0}-[\mathcal{A}{}_{0},\mathcal{A}{}_{i}]-\epsilon_{ij}\mathcal{F}_{j}\,;\phantom{\frac{1}{2}}$
(28)
$\phantom{\frac{1}{2}}\dot{\mathcal{F}}{}_{i}=\partial_{i}\mathcal{F}_{0}+[\mathcal{A}{}_{i},\mathcal{F}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{F}{}_{i}]-\epsilon_{ij}\mathcal{G}_{j}\,;\phantom{\frac{1}{2}}$
(29)
$\begin{array}[]{c}\displaystyle\phantom{\sum_{1}^{2}}\dot{\mathcal{G}}{}_{i}=\partial_{i}\mathcal{G}_{0}+\frac{\alpha_{2}}{\alpha_{3}}\epsilon_{ij}\mathcal{G}_{j}+\frac{\alpha_{1}}{\alpha_{3}}\epsilon_{ij}\mathcal{F}_{j}+[\mathcal{A}{}_{i},\mathcal{G}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{G}{}_{i}]+\phantom{\sum_{1}^{2}}\\\\[14.22636pt]
\displaystyle+\frac{\alpha{}_{3}{}^{2}}{C(\beta;\alpha)}[\beta{}_{1}\mathcal{F}{}_{0}+\beta{}_{2}\mathcal{G}{}_{0},\beta{}_{1}\mathcal{F}{}_{i}+\beta{}_{2}\mathcal{G}{}_{i}]\,.\end{array}$
(30)

The evolutionary equations (28), (29), (30) are supplemented by the constraint

$\begin{array}[]{c}\displaystyle\Theta=\epsilon_{ij}(\alpha_{1}\partial_{i}\mathcal{A}{}_{j}+\alpha_{2}\partial_{i}\mathcal{F}{}_{j}+\alpha_{3}\partial_{i}\mathcal{G}{}_{j}+[\mathcal{A}{}_{i},\frac{1}{2}\alpha_{1}\mathcal{A}{}_{j}+\alpha{}_{2}\mathcal{F}{}_{j}+\alpha{}_{3}\mathcal{G}{}_{j}]+\\\\[14.22636pt]
\displaystyle-\frac{\alpha_{3}}{2C(\beta;\alpha)}[\beta_{1}\mathcal{F}{}_{i}+\beta_{2}\mathcal{G}{}_{i},\beta_{1}\mathcal{F}{}_{j}+\beta_{2}\mathcal{G}{}_{j}])\approx
0\,.\end{array}$ (31)

The systems (16) and (28), (29), (30), (31) are equivalent. Evolutionary
equations (28), (29) express the auxiliary fields
$\mathcal{F}{}_{i},\mathcal{G}{}_{i}$ in terms of time derivatives of the
original vector potential $\mathcal{A}_{\mu}$. Solving them with respect to
the unknown $\mathcal{F}_{i}$, $\mathcal{G}_{i}$, we obtain relations (26).
Equations (30), (31) represent the space and time components of the original
higher derivative system (16) where all the higher derivatives of the vector
potential are expressed in terms of extra variables $\mathcal{F}{}_{i}$,
$\mathcal{G}{}_{i}$. Equation (31) serves as constraint, which does not
involve time derivatives. Once it is preserved the time evolution (see the
gauge identity (19)), no secondary constraints are imposed on the fields.

We associate the on-shell Hamiltonian $\mathcal{H}_{0}$ with the 00-component
of the conserved tensor (20). In terms of the phase-space variables
$\mathcal{A}{}_{i}$, $\mathcal{F}{}_{i}$, $\mathcal{G}{}_{i}$, it reads

$\begin{array}[]{c}\displaystyle\mathcal{H}_{0}=\frac{1}{2}\text{tr}\Big{\\{}\beta_{2}\alpha_{3}(\mathcal{G}{}^{0}\mathcal{G}{}^{0}+\mathcal{G}{}^{i}\mathcal{G}{}^{i})+2\beta_{1}\alpha_{3}(\mathcal{G}{}^{0}\mathcal{F}{}^{0}+\mathcal{G}{}^{i}\mathcal{F}{}^{i})+\\\\[14.22636pt]
\displaystyle+(\beta_{1}\alpha_{2}-\beta_{2}\alpha_{1})(\mathcal{F}{}^{0}\mathcal{F}{}^{0}+\mathcal{F}{}^{i}\mathcal{F}{}^{i})\Big{\\}}\,.\end{array}$
(32)

Here, the functions $\mathcal{F}{}^{0}$, $\mathcal{G}{}^{0}$ are defined in
(27). The on-shell Hamiltonian $\mathcal{H}_{0}$ depends on the free model
parameters $\alpha_{k},k=1,2,3$ and coupling constants $\beta_{l},l=1,2$. Off-
shell, The Hamiltonian is a sum of on-shell part (32) and a linear combination
of constraints. We chose the following ansatz for the total Hamiltonian,

$\begin{array}[]{c}\displaystyle\mathcal{H}=\mathcal{H}_{0}+\text{tr}\Big{[}\frac{C(\beta;\alpha)}{\beta{}_{2}\alpha{}_{2}-\beta{}_{1}\alpha{}_{3}}\mathcal{A}_{0}-\frac{\beta_{2}\beta_{1}\alpha_{3}}{\beta_{2}\alpha_{2}-\beta_{1}\alpha_{3}}\mathcal{F}_{0}-\frac{\beta_{2}{}^{2}\alpha_{3}}{\beta_{2}\alpha_{2}-\beta_{1}\alpha_{3}}\mathcal{G}_{0}\Big{]}\Theta\,,\end{array}$
(33)

where $\beta,\alpha$ are constant parameters. The on-shell vanishing terms,
being included into the Hamiltonian, does not contribute to the equations of
motion of gauge-invariant quantities. The equations of motion do alter for
non-gauge-invariant variables. As we are seeking for the Hamiltonian and
Poisson bracket literally reproducing the first-order form (28), (29), (30),
(31) of equations (16), the on-shell vanishing contributions are kept under
control in (33).

The Hamiltonian is well-defined if the parameters $\beta,\alpha$ subject to
the following conditions:

$C(\beta;\alpha)\neq 0\,,\qquad\beta_{1}\alpha_{3}-\beta_{2}\alpha_{2}\neq 0.$
(34)

These relations have a clear origin. The first condition in this set insures
that the on-shell Hamiltonian is a non-degenerate quadratic form of the phase-
space variables $\mathcal{A}{}_{i}$, $\mathcal{F}{}_{i}$, $\mathcal{G}{}_{i}$.
This requirement is reasonable because the degenerate Hamiltonian cannot
generate the evolution of all physical degrees of freedom. The second relation
(34) ensures that the numerical factor at the Lagrange multiplier
$\mathcal{A}{}_{0}$ is non-singular. This is necessary to reproduce the
correct gauge transformations for all the dynamical variables. We also note
that the obstructions to existence of Hamiltonian remain valid in the free
limit. This means that the inclusion of interaction by the scheme of Section 3
does not restrict the class models that admit the Hamiltonian formulation. In
the other words, every theory in the class (16) admitting the Hamiltonian
formulation with the Hamiltonian (33) in the free limit, is Hamiltonian at the
interacting level. Hereinafter, we assume that the relations (34) are
satisfied. The Hamiltonian is on-shell bounded from below if the conditions
(24) are satisfied. In this case, we expect to construct the constrained
Hamiltonian formulation with a bounded Hamiltonian. The corresponding quantum
theory has a good chance to be stable.

Now, let us seek for the Poisson bracket between the fields
$\mathcal{A}{}_{i}$, $\mathcal{F}{}_{i}$, $\mathcal{G}{}_{i}$ that brings
equations (28), (29), (30), (31) to the Hamiltonian form (25) with the
Hamiltonian (33). Comparing the right hand sides of relations (28), (29), (30)
and equations (25), we obtain the following system of algebraic equations:

$\phantom{\frac{1}{2}}\Big{\\{}A{}_{i},\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\mathcal{A}_{0}-[\mathcal{A}{}_{0},\mathcal{A}{}_{i}]-\epsilon_{ij}\mathcal{F}_{j}\,;\phantom{\frac{1}{2}}$
(35)
$\phantom{\frac{1}{2}}\Big{\\{}F{}_{i},\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\mathcal{F}_{0}+[\mathcal{A}{}_{i},\mathcal{F}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{F}{}_{i}]-\epsilon_{ij}\mathcal{G}_{j}\,;\phantom{\frac{1}{2}}$
(36)
$\begin{array}[]{c}\displaystyle\phantom{\sum_{1}^{2}}\Big{\\{}G{}_{i},\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\mathcal{G}_{0}+\frac{\alpha_{2}}{\alpha_{3}}\epsilon_{ij}\mathcal{G}_{j}+\frac{\alpha_{1}}{\alpha_{3}}\epsilon_{ij}\mathcal{F}_{j}+[\mathcal{A}{}_{i},\mathcal{G}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{G}{}_{i}]+\phantom{\sum_{1}^{2}}\\\\[14.22636pt]
\displaystyle+\frac{\alpha{}_{3}{}^{2}}{C(\beta;\alpha)}[\beta{}_{1}\mathcal{F}{}_{0}+\beta{}_{2}\mathcal{G}{}_{0},\beta{}_{1}\mathcal{F}{}_{i}+\beta{}_{2}\mathcal{G}{}_{i}]\,.\end{array}$
(37)

The sign $\approx$ means equality modulo constraint (31). The Poisson bracket,
being defined by these equations, depends on five independent arguments: the
free model parameters $\alpha_{3}$, $\alpha_{2}$, $\alpha_{1}$, and coupling
constants $\beta_{2}$, $\beta_{1}$. It has the following form:

$\phantom{\sum_{a}^{b}}\\{\mathcal{G}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{G}{}^{b}{}_{j}(\boldsymbol{y})\\}=\frac{\beta_{1}\alpha_{2}{}^{2}-\beta_{1}\alpha_{1}\alpha_{3}-\beta_{2}\alpha_{2}\alpha_{1}}{\alpha{}_{3}{}^{2}C(\beta;\alpha)}\epsilon_{ij}\delta^{ab}\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,;\phantom{\sum_{a}^{b}}$
(38)
$\phantom{\sum_{a}^{b}}\\{\mathcal{F}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{G}{}^{b}{}_{j}(\boldsymbol{y})\\}=\frac{\beta_{2}\alpha_{1}-\beta_{1}\alpha_{2}}{\alpha_{3}C(\beta;\alpha)}\epsilon_{ij}\delta^{ab}\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,;\phantom{\sum_{a}^{b}}$
(39)
$\phantom{\sum_{a}^{b}}\\{\mathcal{F}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{F}{}^{b}{}_{j}(\boldsymbol{y})\\}=\\{\mathcal{A}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{G}{}^{b}{}_{j}(\boldsymbol{y})\\}=\frac{\beta_{1}}{C(\beta;\alpha)}\epsilon_{ij}\delta^{ab}\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,;\phantom{\sum_{a}^{b}}$
(40)
$\phantom{\sum_{a}^{b}}\\{\mathcal{A}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{F}{}^{b}{}_{j}(\boldsymbol{y})\\}=\frac{-\beta_{2}}{C(\beta;\alpha)}\epsilon_{ij}\delta^{ab}\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,;\phantom{\sum_{a}^{b}}$
(41)
$\phantom{\sum_{a}^{b}}\\{\mathcal{A}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{A}{}^{b}{}_{j}(\boldsymbol{y})\\}=0\,.\phantom{\sum_{a}^{b}}$
(42)

Here,
$\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})=\delta(x^{1}-y^{1})\delta(x^{2}-y^{2})$
is the $2d$ $\delta$-function in the space coordinates. The Poisson bracket
(38), (39), (40), (41), (42) is a covariant generalization of the its free
analog, being derived in [43]. This result is not surprising. The free limit
of expressions (38), (39), (40), (41), (42) is determined by the linear model,
while the Poisson brackets with the polynomial dependence on fields
contradicts the structure of equations (25). In the last case, the Poisson
bracket with the Hamiltonian involves higher powers of fields than the right
hand side of first-order equations (28), (29), (30).

The Poisson bracket (38), (39), (40), (41), (42) is a non-degenerate tensor,
so it has inverse, being a symplectic two-form. The latter defines a
Hamiltonian action functional

$\begin{array}[]{c}\displaystyle
S_{h}=\int\Big{\\{}\text{tr}\Big{(}\frac{C(\beta;\alpha)}{\beta_{2}\alpha_{2}-\beta_{1}\alpha_{3}}\varepsilon_{ij}(\alpha_{1}\mathcal{A}{}_{i}+2\alpha_{2}\mathcal{F}{}_{i}+2\alpha_{3}\mathcal{G}{}_{i})\dot{\mathcal{A}}{}_{j}+\frac{\beta_{1}{}^{2}\alpha_{3}}{\beta_{2}\alpha_{2}-\beta_{1}\alpha_{3}}\times\\\\[14.22636pt]
\displaystyle\times\varepsilon_{ij}(\beta{}_{1}\mathcal{F}{}_{i}+\beta{}_{2}\mathcal{G}{}_{i})(\beta{}_{1}\mathcal{F}{}_{j}+\beta{}_{2}\mathcal{G}{}_{j})\Big{)}-\mathcal{H}\Big{\\}}d\boldsymbol{x}dt\,,\end{array}$
(43)

where $\mathcal{H}$ denotes the Hamiltonian (33). In the case
$\beta_{2}=0\,,\beta_{1}=1$, equation (43) reproduces the Ostrogradski action
for the variational model (17),

$S_{c}=\int\Big{\\{}\text{tr}(\epsilon_{ij}(\alpha_{1}\mathcal{A}{}_{i}+2\alpha_{2}\mathcal{F}{}_{i}+2\alpha_{3}\mathcal{G}{}_{i})\dot{\mathcal{A}}{}_{j}-\epsilon_{ij}\mathcal{F}{}_{i}\dot{\mathcal{F}}{}_{j})-\mathcal{H}{}_{c}\Big{\\}}d\boldsymbol{x}dt\,,$
(44)

where $\mathcal{H}_{c}$ is the canonical Hamiltonian,

$\mathcal{H}{}_{c}=\frac{1}{2}\text{tr}\Big{(}\alpha_{3}(\mathcal{G}{}^{0}\mathcal{F}{}^{0}+\mathcal{G}{}^{i}\mathcal{F}{}^{i})+\alpha_{2}(\mathcal{F}{}^{0}\mathcal{F}{}^{0}+\mathcal{F}{}^{i}\mathcal{F}{}^{i})+\mathcal{A}_{0}\Theta\Big{)}\,.$
(45)

In the free case, the Hamiltonian formulation (43) has been first proposed in
[43]. Later, it has been re-derived in [42] by the direct application of the
Ostrogradski procedure to the action functional (1). The comparison of
conditions (12), (34) implies that the Hamiltonian (33) is on-shell unbounded
for all the Lagrangian interactions. For non-Lagrangian interactions, the
Hamiltonian (33) can be on-shell bounded or unbounded depending on the value
of model parameters. If the Hamiltonian is on-shell bounded, the interacting
theory is stable at the quantum level. To our knowledge, the ECS theory (16)
is a first higher derivative model with a non-abelian gauge symmetry admitting
alternative Hamiltonian formulation with a bounded Hamiltonian. This means
that the concept of stabilization of dynamics by means of alternative
Hamiltonian formalism applies beyond the linear level.

In the free limit, the Hamiltonian (33) and Poisson bracket (38), (39), (40),
(41), (42) depend on the parameters $\beta_{1}$ and $\beta_{2}$, while the
equations of motion (28), (29), (30), (31) do not. This means that the free
ECS theory (1) is a multi-Hamiltonian theory. The general representative the
series of Hamiltonian formulations does not follow from the Ostrogradski
procedure. The reason is that the bounded and unbounded Hamiltonian cannot be
connected by the change of coordinates in the phase space. The number of
entries in the free Hamiltonian series equals $2n$, where $n$ is the number of
color indices. The on-shell Hamiltonian is given by the $00$-component (9) of
the conserved tensor (6). The Poisson bracket between the phase-space
variables $A{}^{a}{}_{i}$, $F{}^{a}{}_{i}$ ,$G{}^{a}{}_{i}$ is determined by
the formula (38), (39), (40), (41), (42), but the parameters $\beta_{1}$,
$\beta{}_{2}$ are replaced by $\beta{}^{a}{}_{1}$, $\beta{}^{a}{}_{2}$ for
each field in the multiplet. At the non-linear level, a selected Poisson
bracket (38), (39), (40), (41), (42) is preserved. If the respective
Hamiltonian is bounded, the Poisson bracket cannot follow from the
Ostrogradski construction. So, the Hamiltonian formalism for the stable non-
linear theories does not follow from the canonical formulation.

In the conclusion of this section, we present the Poisson brackets between the
constraints (38), (39), (40), (41), (42),

$\\{\Theta{}^{a}(\boldsymbol{x}),\Theta{}^{b}(\boldsymbol{y})\\}=\frac{\beta_{2}\alpha_{2}-\beta_{1}\alpha_{3}}{C(\beta;\alpha)}f{}^{abc}\Theta{}^{c}(\boldsymbol{x})\delta(\boldsymbol{x}-\boldsymbol{y}).$
(46)

The Hamiltonian is gauge-invariant,

$\\{\Theta{}^{a}(\boldsymbol{x}),\mathcal{H}{}_{0}(\boldsymbol{y})\\}=0\,,\qquad
a=1,\ldots,n\,.$ (47)

As is seen, the first-order model (43) is gauge theory of special form in the
whole range (34) of admissible values of the parameters $\alpha$, $\beta$.
Equation (43) determines the least action principle for the model. The
quantization of this theory can be performed by means of the well-known
procedures [55]. In all the instances with bounded from below Hamiltonian, we
expect the stability of the model. This means that equations (16) determine a
higher derivative theory, which can be stable at classical and quantum level.
Because of existence of first-order formalism, the model is as good as the
Lagrangian theories. In particular, it admits consequent quantization, and
correspondence between symmetries and conserved quantities. The Hamiltonian
serves as quantity that is associated to the invariance of the model (43) with
respect to the time translations. The Hamiltonian (33) is on-shell bounded if
the conditions (24) are satisfied.

##  5 Resonance case

In case of resonance (positions (2b), (3) of classification on pages 4-5), the
non-linear theory (3) admits the Hamiltonian form of dynamics, but the
Hamiltonian is unbounded in all the instances (conditions (24) are
inconsistent). So, the model has to be considered as unstable. In this
section, we demonstrate that the dynamics of the theory with resonance can be
stabilized by means of inclusion of interaction with extra dynamical field. We
apply the ”Higgs-like” mechanism, which has been first proposed in the context
of study of the ”doubly massless” generalized Podolsky electrodynamics in the
paper [48]. Here, we use it in the theory with non-abelian gauge symmetry for
the first time. We mostly consider the model (16) with the third-order
resonance. We chose the following values for the parameters of free theory
(1): $\alpha_{1}=\alpha_{2}=0$, $\alpha_{3}=-1$. This choice does not restrict
generality because the constant $\alpha_{3}$ accounts the possibility of
multiplication of equations of motion by an overall factor. The wave operator
of free model (3) appears to be cube of the CS operator, so we have a sort of
”triply massless” extended theory.

We start construction of interaction by extending the set of dynamical
variables by a real scalar field $\phi(x)$. The non-linear theory of vector
multiplet $\mathcal{A}{}_{\mu}(x)$ and scalar field $\phi(x)$ is determined by
the equations of motion, which have the following form:

$\begin{array}[]{c}\displaystyle\mathbb{T}_{\mu}=\epsilon_{\mu\nu\rho}\Big{\\{}D{}^{\nu}[(\widetilde{\gamma}{}^{2}\phi^{2}-1)\mathcal{G}{}+\gamma^{2}\phi^{2}\mathcal{F}]{}^{\rho}+\frac{(\widetilde{\gamma}{}^{2}\phi^{2}-1)^{2}}{2\beta_{1}{}^{2}}[\beta_{1}\mathcal{F}^{\nu}+\beta_{2}\mathcal{G}^{\nu},\beta_{1}\mathcal{F}^{\rho}+\beta_{2}\mathcal{G}^{\rho}]\Big{\\}}=0\,,\\\\[14.22636pt]
\displaystyle\mathbb{T}=\Big{\\{}\partial_{\mu}\partial^{\mu}-(m^{2}-\widetilde{\gamma}^{2}\frac{(\beta{}_{1}\mathcal{F}_{\rho}+\beta{}_{2}\mathcal{G}_{\rho})(\beta{}_{1}\mathcal{F}^{\rho}+\beta{}_{2}\mathcal{G}^{\rho})}{\beta_{1}{}^{2}}+\phi^{2}\Big{\\}}\phi=0\,.\end{array}$
(48)

Here, the vectors $\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$ are defined in
(2), and the abbrevialtion
$\widetilde{\gamma}{}^{2}=\gamma{}^{2}\beta_{2}/\beta_{1}$ is used. The
constants $\beta_{2}$, $\beta_{1}$, $\gamma$, $m$ are model parameters, being
real numbers. Throughout the section we assume that $m,\gamma>0$ and
$\beta_{1}\neq 0$. The option $\beta_{1}=0$ is not admissible because the
self-coupling (16) between the vector fields becomes inconsistent in this
case. The value $\beta_{2}$ can be arbitrary real number (positive, negative
or zero). Only positive values of $\beta_{2}$, $\beta_{1}$ lead to stable
couplings. This justifies notation $\widetilde{\gamma}{}^{2}$ for the quantity
$\gamma^{2}\beta_{2}/\beta_{1}$. Without loss of generality, we put a unit
coefficient at $\phi^{3}$-term. An overall factor at the $\phi{}^{3}$-vertex
can be absorbed by the scaling of scalar field, $\phi\mapsto\lambda\phi$ with
the appropriate $\lambda\neq 0$.

Equations (48) have a clear meaning. The first line of the system (48)
describes the motion of the vector multiplet. The linear in the fields term in
equations corresponds to the ”triply massless” ECS theory. The coupling
includes self-interaction term (16), and extra contribution involving scalar
field. It is convenient to think that (48) follows from (16) after the formal
redefinition of the model parameters

$\phantom{\frac{1}{2}}\alpha_{1}=0\,,\qquad\alpha_{2}=\gamma{}^{2}\phi{}^{2}\,,\qquad\alpha_{3}=\widetilde{\gamma}{}^{2}\phi{}^{2}-1\,.\phantom{\frac{1}{2}}$
(49)

The characteristic polynomial (4) for the model reads

$\phantom{\frac{1}{2}}M(\alpha;z)=\gamma^{2}\phi^{2}z^{2}+(\widetilde{\gamma}^{2}\phi^{2}-1)z^{3}\,.\phantom{\frac{1}{2}}$
(50)

If the scalar field is set to nonzero constant from outset, we obtain theory
(16) with the second-order resonance for zero root (case (2a) of
classification on pages 4-5). As it has been explained above, the
corresponding model describes a massless spin-1 field, and a massive spin 1
vector subjected to the self-duality condition. It is stable at the classical
and quantum levels. The second equation in the system (48) describes the
motion of scalar field. This equation includes $\phi^{3}$-coupling, which
ensures the existence of non-zero stationary solution for
$\mathcal{A}_{\mu}=0$. This means that $\phi$ serves as the Higgs field in the
model (48).

The interaction, being defined by equations (48), is consistent. The equations
are preserved by the standard Yang-Mills gauge symmetry (18). The scalar field
is preserved by the gauge transformation. The gauge identity has a slightly
different form,

$\mathbb{D}_{\mu}\mathbb{T}^{\mu}=0\,,\qquad\mathbb{D}_{\mu}=D_{\mu}+\frac{\widetilde{\gamma}{}^{2}\phi{}^{2}-1}{\beta_{1}{}^{2}}[\beta_{1}\mathcal{F}{}_{\mu}+\beta_{2}\mathcal{G}{}_{\mu},\cdot]\,.$
(51)

We note that the gauge generator involves the scalar field explicitly. The
model has the same number of physical degrees of freedom as the free theory
because the orders of equations of motion (48), gauge identities (51), and
gauge symmetries (18) are preserved by coupling. Equations (48) do not follow
from the least action principle from the least action principle for a
functional with higher derivatives unless $\beta_{2}=0$. In the case of
Lagrangian coupling, the action principle reads

$S[\phi(x),\mathcal{A}(x)]=\frac{1}{2}\int\Big{\\{}\mathcal{F}_{\mu}(\mathcal{G}^{\mu}+\gamma{}^{2}\phi^{2}\mathcal{F}^{\mu})+\partial_{\mu}\phi\partial^{\mu}\phi+m^{2}\phi^{2}-\frac{1}{2}\phi^{4}\Big{\\}}d^{3}x\,.$
(52)

The dynamics of the Lagragian theory is unstable. In the case of non-Lagragian
couplings, the dynamics of the model can be stable. In the remaining part of
the section, we address an issue of construction of the constrained
Hamiltonian formalism with a bounded Hamiltonian for the higher derivative
equations (48). Existence of such a formalism implies classical and quantum
stability of the model. To avoid ”tautology” we outline to the most crucial
steps of the construction, while the details are provided above.

We introduce a special notation for the linear combination of generalized
strengths $\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$ entering the free part
equations (48),

$\mathcal{W}{}_{\mu}=(\widetilde{\gamma}{}^{2}\phi{}^{2}-1)\mathcal{G}{}{}_{\mu}+\gamma^{2}\phi{}^{2}\mathcal{F}{}_{\mu}\,.$
(53)

We chose the space components of the vectors $\mathcal{F}{}_{i}$,
$\mathcal{W}{}_{i}$, $i=1,2$ (26) as variables absorbing first and second
derivatives of the original dynamical field $\mathcal{A}{}_{\mu}$. In the
sector of scalar field, we introduce the canonical momentum $\pi=\dot{\phi}$.
We consider $\mathcal{W}{}_{0}$ as a special notation for the combination (53)
of the derivatives of phase-space variables
$\mathcal{A}{}_{\mu},\mathcal{F}{}_{i}$. In terms of the variables $\phi$,
$\pi$ $\mathcal{A}{}_{\mu}$, $\mathcal{F}{}_{i}$, $\mathcal{G}{}_{i}$, the
first-order equations eventually read

$\phantom{\frac{1}{2}}\dot{\mathcal{A}}{}_{i}=\partial_{i}\mathcal{A}_{0}-[\mathcal{A}{}_{0},\mathcal{A}{}_{i}]-\epsilon_{ij}\mathcal{F}_{j}\,;\phantom{\frac{1}{2}}$
(54)
$\phantom{\sum_{1}^{2}}\dot{\mathcal{F}}{}_{i}=\partial_{i}\mathcal{F}_{0}+[\mathcal{A}{}_{i},\mathcal{F}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{F}{}_{i}]-\epsilon_{ij}\frac{\mathcal{W}{}_{j}-\gamma{}^{2}\phi{}^{2}\mathcal{F}{}_{j}}{\widetilde{\gamma}^{2}\phi^{2}-1}\,;\phantom{\sum_{1}^{2}}$
(55)
$\begin{array}[]{c}\displaystyle\dot{\mathcal{W}}{}_{i}=\partial_{i}\mathcal{W}{}_{0}+[\mathcal{A}{}_{i},\mathcal{W}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{W}{}_{i}]+\epsilon_{ij}[\beta_{2}\mathcal{W}{}_{0}-\beta_{1}\mathcal{F}{}_{0},\beta_{2}\mathcal{W}{}_{j}-\beta_{1}\mathcal{F}{}_{j}]\,;\phantom{\sum_{1}^{2}}\end{array}$
(56)
$\begin{array}[]{c}\displaystyle\phantom{\sum_{1}^{2}}\dot{\pi}=\Big{\\{}\partial_{i}\partial_{i}-m^{2}+\widetilde{\gamma}{}^{2}\frac{(\beta_{2}\mathcal{W}{}_{\mu}-\beta_{1}\mathcal{F}_{\mu})(\beta_{2}\mathcal{W}{}^{\mu}-\beta_{1}\mathcal{F}^{\mu})}{(\widetilde{\gamma}{}^{2}\phi{}^{2}-1)^{2}}+\phi^{2}\Big{\\}}\phi\,;\phantom{\sum_{1}^{2}}\end{array}$
(57) $\phantom{\sum_{1}^{2}}\dot{\phi}=\pi\,.\phantom{\sum_{1}^{2}}$ (58)

The evolutionary equations (54), (55), (56), (57), (58) are supplemented by
the constraint

$\begin{array}[]{c}\displaystyle\Theta=\epsilon_{ij}\Big{(}\partial_{i}\mathcal{W}{}_{j}+[\mathcal{A}{}_{i},\mathcal{W}{}_{j}]+\frac{1}{2\beta_{2}{}^{2}}[\beta_{2}\mathcal{W}{}_{i}-\beta_{1}\mathcal{F}{}_{i},\beta_{2}\mathcal{W}{}_{j}-\beta_{1}\mathcal{F}{}_{j}]\Big{)}=0\,.\end{array}$
(59)

The constraint conserves on-shell if the identity (51) is taken into account.
The higher derivative system (48) follows from equations (54), (55), (56),
(59) if all the extra variables are expressed in terms of derivatives of
original dynamical fields $\phi,\mathcal{A}$.

The model (48) admits a second-rank symmetric conserved tensor,

$\begin{array}[]{c}\displaystyle\boldsymbol{\Theta{}^{\mu\nu}}(\beta;\alpha)=\text{tr}\Big{\\{}\frac{\widetilde{\gamma}^{2}\phi^{2}-1}{\beta_{2}}\Big{(}(\beta_{1}\mathcal{F}^{\mu}+\beta_{2}\mathcal{G}^{\mu})(\beta_{1}\mathcal{F}^{\nu}+\beta_{2}\mathcal{G}^{\nu})-\frac{}{2}g^{\mu\nu}(\mathcal{F}{}_{\rho}+\mathcal{G}{}_{\rho})(\mathcal{F}{}^{\rho}+\mathcal{G}{}^{\rho})\Big{)}+\\\\[14.22636pt]
\displaystyle+\frac{\beta_{2}{}^{2}}{\beta_{1}}\Big{(}\mathcal{F}{}^{\mu}\mathcal{F}{}^{\nu}-\frac{1}{2}g^{\mu\nu}\mathcal{F}{}_{\rho}\mathcal{F}{}^{\rho}\Big{)}\Big{\\}}+\partial^{\mu}\phi\partial^{\nu}\phi+g^{\mu\nu}(-\frac{1}{2}\partial_{\rho}\phi\partial^{\rho}\phi-\frac{1}{2}m^{2}\phi^{2}+\frac{1}{4}\phi^{4})\,;\end{array}$
(60)
$\phantom{\sum_{1}^{2}}\partial_{\nu}\boldsymbol{\Theta{}^{\mu\nu}}(\beta;\alpha)=-\text{tr}(\epsilon^{\mu\nu\rho}(\beta_{1}\mathcal{F}{}_{\nu}+\beta_{2}\mathcal{G}{}_{\nu})\mathbb{T}{}_{\rho})+\partial^{\mu}\phi\cdot\mathbb{T}\,.\phantom{\sum_{1}^{2}}$
(61)

The Hamiltonian is given by a sum of the $00$-component of conserved tensor,
being expressed in terms of variables $\phi$, $\pi$, $\mathcal{A}{}_{\mu}$,
$\mathcal{F}{}_{i}$, $\mathcal{W}{}_{i}$ (26), (53), and a constraint term,

$\begin{array}[]{c}\displaystyle\mathcal{H}=\text{tr}\Big{\\{}\frac{(\beta_{2}\mathcal{W}{}^{0}-\beta_{1}\mathcal{F}{}^{0})(\beta_{2}\mathcal{W}{}^{0}-\beta_{1}\mathcal{F}{}^{0})+(\beta_{2}\mathcal{W}{}^{i}-\beta_{1}\mathcal{F}{}^{i})(\beta_{2}\mathcal{W}{}^{i}-\beta_{1}\mathcal{F}{}^{i})}{2(\widetilde{\gamma}^{2}\phi^{2}-1)\beta_{2}}+\\\\[14.22636pt]
\displaystyle+\frac{\beta_{1}{}^{2}}{2\beta_{2}}(\mathcal{F}{}^{0}\mathcal{F}{}^{0}+\mathcal{F}{}^{i}\mathcal{F}{}^{i})-(\mathcal{A}{}_{0}-\beta_{1}\mathcal{F}_{0}+\beta_{2}\mathcal{W}_{0})\Theta\Big{\\}}+\\\\[14.22636pt]
\displaystyle\phantom{\frac{1}{2}}+\frac{1}{2}(\pi\pi+\partial^{i}\phi\partial^{i}\phi)-\frac{1}{2}m^{2}\phi^{2}+\frac{1}{4}\phi^{4}\,.\phantom{\frac{1}{2}}\end{array}$
(62)

Here, the quantities $\mathcal{F}{}_{0},\mathcal{W}{}_{0}$ denote
abbreviations (27), (53), and the numbers $m,\gamma$ are model parameters. The
Hamiltonian is on-shell bounded if

$\widetilde{\gamma}{}^{2}\phi^{2}-1>0\,,\qquad\beta_{2}>0\,.$ (63)

These conditions involve the scalar field $\phi$. Once the initial value of
$\phi$ is a Cauchy data for equations (48), the Hamiltonian cannot be globally
bounded. However, the Hamiltonian is given by a positive definite quadratic
form in the variables $\mathcal{F}$, $\mathcal{G}$ in the range
$|\widetilde{\gamma}\phi|>1$ of values of scalar field. This corresponds to
the case of stability island.

The Poisson bracket between the fields $\phi$, $\pi$, $\mathcal{A}{}_{\mu}$,
$\mathcal{F}{}_{i}$, $\mathcal{G}{}_{i}$ is determined by the following system
of equations:

$\phantom{\frac{1}{2}}\Big{\\{}\mathcal{A}{}_{i},\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\mathcal{A}_{0}-[\mathcal{A}{}_{0},\mathcal{A}{}_{i}]-\epsilon_{ij}\mathcal{F}_{j}\,;\phantom{\frac{1}{2}}$
(64)
$\begin{array}[]{c}\displaystyle\phantom{\frac{1}{2}}\Big{\\{}\mathcal{F}{}_{i},\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\mathcal{F}_{0}+[\mathcal{A}{}_{i},\mathcal{F}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{F}{}_{i}]-\epsilon_{ij}\frac{\mathcal{W}{}_{j}-\gamma{}^{2}\phi{}^{2}\mathcal{F}{}_{j}}{\widetilde{\gamma}{}^{2}\phi{}^{2}-1}\,;\phantom{\sum_{1}^{2}}\end{array}$
(65)
$\begin{array}[]{c}\displaystyle\Big{\\{}G{}_{i},\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\mathcal{W}{}_{0}+[\mathcal{A}{}_{i},\mathcal{W}{}_{0}]-[\mathcal{A}{}_{0},\mathcal{G}{}_{i}]+\frac{\epsilon_{ij}}{\beta_{1}{}^{2}}[\beta_{2}\mathcal{W}{}_{0}-\beta_{1}\mathcal{F}{}_{0},\beta_{2}\mathcal{W}{}_{j}-\beta_{1}\mathcal{F}{}_{j}]\,.\end{array}$
(66)
$\begin{array}[]{c}\displaystyle\phantom{\sum_{1}^{2}}\Big{\\{}\pi,\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\partial_{i}\partial_{i}\phi+\Big{(}m^{2}-\widetilde{\gamma}^{2}\frac{(\beta_{2}\mathcal{W}{}_{\mu}-\beta_{1}\mathcal{F}{}_{\mu})(\beta_{2}\mathcal{W}{}^{\mu}-\beta_{1}\mathcal{F}{}^{\mu})}{(\widetilde{\gamma}^{2}\phi^{2}-1){}^{2}}-\phi^{2}\Big{)}\phi\,;\phantom{\sum_{1}^{2}}\end{array}$
(67)
$\Big{\\{}\phi,\int\mathcal{H}(\boldsymbol{y})d\boldsymbol{y}\Big{\\}}\approx\pi.$
(68)

All the equalities are considered modulo terms that vanish modulo constraint
(59). The solution for the equations (64), (65), (66), (67), (68) reads

$\phantom{\sum_{a}^{b}}\phantom{\sum_{a}^{b}}\\{\mathcal{W}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{W}{}^{b}{}_{j}(\boldsymbol{y})\\}=\\{\mathcal{F}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{W}{}^{b}{}_{j}(\boldsymbol{y})\\}=\\{\mathcal{A}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{F}{}^{b}{}_{j}(\boldsymbol{y})\\}=0\,;\phantom{\sum_{a}^{b}}$
(69)
$\phantom{\sum_{a}^{b}}\phantom{\sum_{a}^{b}}\\{\mathcal{A}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{W}{}^{b}{}_{j}(\boldsymbol{y})\\}=-\\{\mathcal{F}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{F}{}^{b}{}_{j}(\boldsymbol{y})\\}=\frac{1}{\beta_{1}}\epsilon_{ij}\delta^{ab}\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,;\phantom{\sum_{a}^{b}}$
(70)
$\phantom{\sum_{a}^{b}}\phantom{\sum_{a}^{b}}\\{\mathcal{A}{}^{a}{}_{i}(\boldsymbol{x}),\mathcal{F}{}^{b}{}_{j}(\boldsymbol{y})\\}=\frac{\beta_{2}}{\beta_{1}{}^{2}}\epsilon_{ij}\delta^{ab}\delta^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,;\phantom{\sum_{a}^{b}}$
(71)
$\phantom{\sum_{a}^{b}}\\{\phi(\boldsymbol{x}),\pi(\boldsymbol{y})\\}=\delta{}^{(2)}(\boldsymbol{x}-\boldsymbol{y})\,.\phantom{\sum_{a}^{b}}$
(72)

All the Poisson brackets between $\phi$, $\pi$ and $\mathcal{A}{}_{i}$,
$\mathcal{F}{}_{i}$, $\mathcal{G}_{i}$ vanish. The relations (69), (70), (71)
can be obtained from (38), (39), (40), (41) after substitution (49). We note
that the Poisson bracket between the fields $\mathcal{A}{}_{i}$,
$\mathcal{F}{}_{i}$, $\mathcal{W}{}_{i}$ is constant, even though the theory
involves extra fields $\pi,\phi$.

Now, we can demonstrate the phenomenon of dynamics stabilization by means of
”Higgs-like” mechanism proposed in [48]. Equations (48) admit a nonzero
stationary solution

$\phantom{\sum_{a}^{b}}\mathcal{A}(\boldsymbol{x})=\mathcal{F}(\boldsymbol{x})=\mathcal{G}(\boldsymbol{x})=0\,,\qquad\pi(\boldsymbol{x})=0\,,\qquad\phi(\boldsymbol{x})=\pm
m\,.\phantom{\sum_{a}^{b}}$ (73)

Introducing notation $\phi^{\ast}(\boldsymbol{x})=\phi(\boldsymbol{x})\mp m$
and expanding equations (48) in the vicinity of this vacua, we obtain the
linearized equations for the fields

$\begin{array}[]{c}\displaystyle\phantom{\sum_{a}^{b}}\mathbb{T}{}_{\mu}=(m{}^{2}\widetilde{\gamma}{}^{2}-1)\mathcal{K}{}_{\mu}+m{}^{2}\gamma{}^{2}\mathcal{G}{}_{\mu}+\ldots=0\,,\phantom{\sum_{a}^{b}}\\\\[14.22636pt]
\displaystyle\phantom{\sum_{a}^{b}}\mathbb{T}=(\partial_{\mu}\partial^{\mu}+5m{}^{2})\phi{}^{\ast}+\ldots=0\,.\phantom{\sum_{a}^{b}}\end{array}$
(74)

The dots denote the terms that are at least quadratic in the fields. As we
see, the dynamics of the field $\mathcal{A}$ is described by the higher
derivative equation (3) with $\alpha_{3}=m{}^{2}\widetilde{\gamma}{}^{2}-1$,
$\alpha_{2}=m{}^{2}\gamma{}^{2}$, $\alpha_{1}=0$. In this cases, the system
has a second order resonance for zero root, which does not affect the
stability of dynamics, at least at the free level. The on-shell Hamiltonian
reads

$\begin{array}[]{c}\displaystyle\mathcal{H}=\text{tr}\Big{\\{}\frac{1}{2(m{}^{2}\widetilde{\gamma}{}^{2}-1)\beta_{2}}((\beta_{1}\mathcal{W}^{0}-\beta_{1}\mathcal{F}^{0})(\beta_{1}\mathcal{W}^{0}-\beta_{1}\mathcal{F}^{0})+(\beta_{2}\mathcal{W}{}^{i}-\beta_{1}\mathcal{F}{}^{i})(\beta_{2}\mathcal{W}{}^{i}-\beta_{1}\mathcal{F}{}^{i}))+\\\\[14.22636pt]
\displaystyle+\frac{\beta_{1}{}^{2}}{2\beta_{2}}(\mathcal{F}{}^{0}\mathcal{F}{}^{0}+\mathcal{F}{}^{i}\mathcal{F}{}^{i})-(\mathcal{A}{}_{0}-\beta_{1}\mathcal{F}_{0}+\beta_{2}\mathcal{W}_{0})\Theta\Big{\\}}+\\\\[14.22636pt]
\displaystyle+\frac{1}{2}(\pi\pi+\partial{}^{i}\phi{}^{\ast}\partial{}^{i}\phi{}^{\ast}+5m{}^{2}(\phi{}^{\ast}){}^{2})+\ldots\phantom{\frac{1}{2}}\end{array}$
(75)

The dots denote the terms that are at least cubic in the fields. The quadratic
part of the Hamiltonian is on-shell bounded if $m\widetilde{\gamma}>1$. This
means that the dynamics of small fluctuations in the vicinity of vacua (73) is
stable. Once local stability is sufficient for construction of stable quantum
theory, the ”Higgs-like” mechanism can stabilize the dynamics of the non-
linear theory (16) with the third-order resonance.

##  6 Conclusion

In this article, we have studied the issue of stability of the most general
gauge ECS theory of third order for a vector multiplet. We have shown that the
free theory with $n$ dynamical fields admits a $2n$-parameter series of
symmetric conserved tensors of second rank. The series includes a canonical
energy-momentum. Even though the canonical energy of the model is unbounded
from below in all the instances, the other conserved tensors have bounded
$00$-component. The stability conditions for dynamics are satisfied if the
model has no resonance (the characteristic equation has simple real roots, the
spectrum of masses is not degenerate). The only exception is the possible
multiplicity two zero root, which corresponds to the Maxwell mode. We have
confirmed the stability of dynamics by construction of the non-canonical (non-
Ostrogradski) Hamiltonian formalism with the Hamiltonian, being given by the
$00$-component of a selected representative in the conserved tensor series. We
have observed that a bounded Hamiltonian is admitted by the model with no
resonance or order two resonance for the zero root.

At next step of our analysis, the self-interaction vertex with the Yang-Mills
gauge symmetry is proposed for the higher derivative theory. The interaction
preserves a selected conserved tensor of free model, being determined by the
values of coupling constants. The Lagrangian interaction, which has been
previously proposed in [42], leads to unstable dynamics. The couplings that
preserve the conserved tensors with bounded $00$-component are non-Lagrangian.
The non-linear dynamics admits a Hamiltonian form. On-shell, the Hamiltonian
is given by the $00$-component of the conserved tensor. In case of Lagrangian
interacting theory, the Hamiltonian follows from the Ostrogradski procedure.
In the case of non-Lagrangian interactions, we have constructed an alternative
Hamiltonian formalism for a system of higher-derivative equations. The
alternative Hamiltonian is on-shell bounded if a conserved tensor with the
bounded $00$-component is preserved at the interacting level. This observation
demonstrates that the stability of higher derivative dynamics can be
consistent with non-abelian gauge symmetry, even though the canonical energy
of the model is unbounded. To our knowledge the proposed theory (16) is a
first stable higher derivative model with non-abelian gauge symmetry.

In the final part of the paper, we have considered the theory with the third-
order resonance, which can be called the ”triply massless” ECS model. This
theory admits non-abelian interactions with Yang-Mills gauge symmetry, but the
$00$-component of the conserved quantity is unbounded in all instances. To
solve the problem, we apply the ”Higgs-like” mechanism of the paper [48]. An
extra scalar field with the $\phi^{3}$ (at the level of equations of motion)
self-coupling is introduced. The dynamics of the model admits the Hamiltonian
form, while the Hamiltonian has a local minimum for nonzero value of the Higgs
field. The dynamics of small fluctuations in the vicinity of energy minimum is
stable.

Acknowledgments. The authors thank A.A. Sharapov for valuable discussions of
this work and comments on the manuscript. The part of the work concerning
construction of stable interactions in the ECS model (Sections 2,3) is
supported by the Ministry of Science and Higher Education of the Russian
Federation, Project No. 0721- 2020-0033. The construction of Hamiltonian
formalism (Sections 4,5) is worked out with the support of the Russian Science
Foundation grant 18-72-10123 in association with the Lebedev Physical
Institute of RAS. The work of SLL was partially supported by the Foundation
for the Advancement of Theoretical Physics and Mathematics ”BASIS”.

##  Appendix A Uniqueness of consistent interaction vertex

In this Appendix, we demonstrate the fact of uniqueness of interaction vertex
(16) in the class of Poincare-covariant couplings that are polynomial in the
invariants $\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$ without higher
derivatives. We apply the method of inclusion of not necessarily Lagrangian
consistent interactions of article [53].

Let us first explain the concept of interaction consistency in class of non-
Lagragian field theories. The dynamics of theory is determined by a system of
partial differential equations (equations of motion) imposed onto the
dynamical fields $\varphi^{I}(x)$,

$\phantom{\frac{1}{2}}\mathbb{T}_{a}(\varphi^{I}(x),\partial\varphi^{I}(x),\partial^{2}\varphi^{I}(x),\ldots)=0\,.\phantom{\frac{1}{2}}$
(76)

The equations of motion do not necessarily follow from the least action
principle for any functional $S[\varphi(x)]$, so $I$ and $a$ may run over
different sets. For the free model, the left hand side of (76) is supposed to
be linear in the fields. The interactions are associated with the deformations
of system (76) by non-linear terms. The equations of motion for the theory
with coupling are polynomial in the fields,

$\phantom{\frac{1}{2}}\mathbb{T}{}_{a}=\mathbb{T}{}^{(0)}{}_{a}+\mathbb{T}{}^{(1)}{}_{a}+\mathbb{T}{}^{(2)}{}_{a}+\ldots=0\,.\phantom{\frac{1}{2}}$
(77)

Here, $\mathbb{T}^{(0)}$, $\mathbb{T}^{(1)}$, and $\mathbb{T}^{(2)}$ are
linear, quadratic, and cubic in the dynamical variables. Throughout the
section, the system (76) (or, equivalently, (77)) is supposed to be
involutive. The concept of involution implies that equations (76) have no
differential consequences of lower order (hidden integrability conditions).
The ECS equations (16) are involutive.

The defining relations for gauge symmetries and gauge identifies in the system
of partial derivative equations (76) read

$\phantom{\frac{1}{2}}\delta_{\epsilon}\varphi^{I}=R{}^{I}{}_{\alpha}\epsilon^{\alpha},\qquad\delta_{\epsilon}\mathbb{T}_{a}|_{\mathbb{T}=0}=0\,;\phantom{\frac{1}{2}}$
(78) $\phantom{\frac{1}{2}}L{}^{a}{}_{A}\mathbb{T}_{a}\equiv
0\,.\phantom{\frac{1}{2}}$ (79)

Here $R{}^{I}{}_{\alpha},L{}^{a}{}_{A}$ are certain differential operators.
For non-Lagrangian equations, the gauge symmetries and gauge identities are
not related to each other, so the multi-indices $A,\alpha$ are different. In
the perturbative setting, the gauge symmetry and gauge identity generators are
supposed to be polynomial in fields

$R{}^{I}{}_{\alpha}=R{}^{(0)}{}^{I}{}_{\alpha}+R{}^{(1)}{}^{I}{}_{\alpha}+R{}^{(2)}{}^{I}{}_{\alpha}+\ldots\,,\qquad
L{}^{a}{}_{A}=L{}^{(0)}{}^{a}{}_{A}+L{}^{(1)}{}^{a}{}_{A}+L{}^{(2)}{}^{a}{}_{A}+\ldots$
(80)

Here, $R^{(0)},L^{(0}$, $R^{(1)},L^{(1)}$, and $R^{(2)},L^{(2)}$ are field-
independent, linear, and quadratic in the dynamical variables. The interaction
is consistent if all the gauge symmetries and gauge identities of free model
are preserved by the deformation of equations of motion. Equation (8) of paper
[53] provides a simple formula for computation of physical degrees of freedom
number based on the orders of derivatives involved into equations of models,
gauge symmetries, and gauge identities. The free and non-linear theory must
have one and the same number of physical degrees of freedom.

Relations (78), (79) imply the following consistency conditions for
$\mathbb{T}^{(k)}$, $R^{(k)}$, $L^{(k)}$, $k=0,1,2,\ldots$:

$\phantom{\frac{1}{2}}R^{(0)}{}^{I}{}_{\alpha}\partial_{I}\mathbb{T}{}^{(0)}{}_{a}=0\,;\phantom{\frac{1}{2}}$
(81)
$\phantom{\frac{1}{2}}R^{(0)}{}^{I}{}_{\alpha}\partial_{I}\mathbb{T}{}^{(1)}{}_{a}+R^{(1)}{}^{I}{}_{\alpha}\partial_{I}\mathbb{T}{}^{(0)}{}_{a}=0\,;\phantom{\frac{1}{2}}$
(82)
$\phantom{\frac{1}{2}}R^{(0)}{}^{I}{}_{\alpha}\partial_{I}\mathbb{T}{}^{(2)}{}_{a}+R^{(1)}{}^{I}{}_{\alpha}\partial_{I}\mathbb{T}{}^{(1)}{}_{a}+R^{(2)}{}^{I}{}_{\alpha}\partial_{I}\mathbb{T}{}^{(0)}{}_{a}=0\,\phantom{\frac{1}{2}}$
(83)

(the symbol $\partial_{I}$ denotes a variational derivative with respect to
the field $\varphi{}^{I}$);

$\phantom{\frac{1}{2}}L^{(0)}{}^{a}{}_{A}\mathbb{T}{}^{(0)}{}_{a}=0\,;\phantom{\frac{1}{2}}$
(84)
$\phantom{\frac{1}{2}}L^{(0)}{}^{a}{}_{A}\mathbb{T}{}^{(1)}{}_{a}+L^{(1)}{}^{a}{}_{A}\mathbb{T}{}^{(0)}{}_{a}=0\,;\phantom{\frac{1}{2}}$
(85)
$\phantom{\frac{1}{2}}L^{(0)}{}^{a}{}_{A}\mathbb{T}{}^{(2)}{}_{a}+L^{(1)}{}^{a}{}_{A}\mathbb{T}{}^{(1)}{}_{a}+L^{(2)}{}^{a}{}_{A}\mathbb{T}{}^{(0)}{}_{a}=0\,.\phantom{\frac{1}{2}}$
(86)

The equations (81), (84) determine the gauge symmetry and gauge identity
generators $R^{(0)}$, $L^{(0)}$ of the free theory. These quantities are
usually given from outset. Relations (82), (85) determine the first-order
corrections to the equations of motion $T^{(1)}$, gauge symmetry generators
$R^{(1)}$, and gauge identity generators $L^{(1)}$. Relations (82), (85)
determine the second-order corrections to the equations of motion $T^{(2)}$,
gauge symmetry generators $R^{(2)}$, and gauge identity generators $L^{(2)}$.
The procedure of interaction construction can be extended to the third and
higher orders. Once the most general covariant ansatz is applied for
$T^{(k)}$, $k=1,2,\ldots$ the procedure (81), (82), (83), (84), (85), (86),
$\ldots$ allows a complete classification of consistent interactions in a
given field theory. An important subtlety of this procedure is that some
lower-order couplings can be inconsistent at the higher orders of perturbation
theory. The first critical step is the extension of the first-order
(quadratic) interaction vertex at the second-order of perturbation theory.

Relation (3) determines the left hand side of the free ECS equations
$\mathbb{T}{}^{(0)}$,

$\phantom{\frac{1}{2}}\mathbb{T}^{(0)}{}^{a}{}_{\mu}=\alpha_{1}F{}^{a}{}_{\mu}+\alpha_{2}G{}^{a}{}_{\mu}+\alpha_{3}K{}^{a}{}_{\mu}=0\,.\phantom{\frac{1}{2}}$
(87)

The gauge symmetries and gauge identities are defined by the gradient and
divergence operators,

$\phantom{\frac{1}{2}}R{}^{(0)}{}_{\mu}=\partial{}_{\mu}\,,\qquad
L^{(0)}{}{}^{\mu}=\partial{}^{\mu}\,.\phantom{\frac{1}{2}}$ (88)

The free gauge identity (81) and free gauge transformation read

$\phantom{\frac{1}{2}}\partial{}_{\mu}\mathbb{T}{}^{(0)}{}^{\mu}\equiv
0\,,\qquad\delta_{\epsilon}A{}^{a}{}_{\mu}=\partial_{\mu}\zeta{}^{a}\,,\phantom{\frac{1}{2}}$
(89)

where $\zeta$’s are gauge transformation parameters. We consider the Poincare-
covariant interactions such that are expressed in terms of gauge covariants
$\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$, $\mathcal{K}{}_{\mu}$ (14) with
no higher-derivative terms being included into coupling. In this case, the
equations of motion are automatically preserved by the Yang-Mills gauge
symmetry (equations (81), (82), (83) are satisfied). Consistent interaction
vertices of first and second orders are selected by the conditions (85), (86).
We elaborate on this problem below.

We assume that the equations of motion are polynomial in gauge covariants
$\mathcal{F}{}_{\mu}$, $\mathcal{G}{}_{\mu}$, $\mathcal{K}{}_{\mu}$. The
linear term is given by the covariantization of the free equations (3). The
most general covariant first-order interaction vertex without higher-
derivatives reads

$\begin{array}[]{c}\displaystyle\phantom{\frac{1}{2}}\mathbb{T}{}^{(0)}{}_{\mu}+\mathbb{T}{}^{(1)}{}_{\mu}=\phantom{\frac{1}{2}}\\\\[19.91692pt]
\displaystyle=\alpha_{1}\mathcal{F}{}_{\mu}+\alpha_{2}\mathcal{G}{}_{\mu}+\alpha_{3}\mathcal{K}{}_{\mu}+\epsilon_{\mu\nu\rho}\Big{(}\frac{1}{2}k{}_{1}[\mathcal{F}{}^{\nu},\mathcal{F}{}^{\rho}]+k{}_{2}[\mathcal{F}{}^{\nu},\mathcal{G}{}^{\rho}]+\frac{1}{2}k{}_{3}[\mathcal{G}{}^{\nu},\mathcal{G}{}^{\rho}]\Big{)}\,,\end{array}$
(90)

where $k_{l},l=1,2,3$ are constants. The covariant divergence of equations of
motion reads

$\begin{array}[]{c}\displaystyle
D_{\mu}\mathbb{T}^{\mu}=-\frac{1}{\alpha_{3}}[k_{2}\mathcal{F}{}_{\mu}+k_{3}\mathcal{G}{}_{\mu},\mathbb{T}{}^{\mu}]+\\\\[14.22636pt]
\displaystyle+\frac{1}{\alpha_{3}}(\alpha{}_{3}{}^{2}-\alpha_{3}k_{1}+\alpha_{2}{}k_{2}-\alpha_{1}{}k_{3})[\mathcal{F}{}_{\mu},\mathcal{G}^{\mu}]+(k_{2}^{2}-k_{3}k_{1})\varepsilon_{\mu\nu\rho}[\mathcal{F}{}^{\mu},[\mathcal{F}{}^{\nu},\mathcal{G}{}^{\rho}]]\,.\end{array}$
(91)

The gauge identity is satisfied in the first-order approximation (85) if the
coefficients $k_{l},l=1,2,3$ satisfy relation

$\phantom{\frac{1}{2}}\alpha{}_{3}{}^{2}-\alpha_{3}k_{1}+\alpha_{2}{}k_{2}-\alpha_{1}{}k_{3}=0\,.\phantom{\frac{1}{2}}$
(92)

The general solution to this equations reads

$k_{1}=-\frac{\beta{}_{1}{}^{2}\alpha_{3}{}^{2}}{C(\beta;\alpha)}+\alpha{}_{1}\beta{}_{3}\,,\qquad
k_{2}=-\frac{\beta{}_{2}{}\beta{}_{1}\alpha_{3}{}^{2}}{C(\beta;\alpha)}\,\qquad
k_{3}=-\frac{\beta{}_{2}{}^{2}\alpha_{3}{}^{2}}{C(\beta;\alpha)}+\alpha{}_{3}\beta{}_{1}\,,$
(93)

where $\beta_{1}$, $\beta_{2}$, $\beta_{3}$ are coupling parameters. The
parameters $\beta_{1},\beta_{2}$ determine the coupling vertex (16) (two
constants determine a single coupling because the ratio $\beta_{1}/\beta_{2}$
is relevant). The constant $\beta_{3}$ is responsible for another interaction
vertex, which is consistent at the first order of perturbation theory. The
interaction vertex (16) is self-consistent with no higher-order corrections
required for the equations of motion. The other coupling needs cubic in the
fields corrections to the equations of motion. To prove the uniqueness of
interaction (16) we should demonstrate that the ansatz (90), (93) is
inconsistent at the second order of perturbation theory for $\beta_{3}\neq 0$.

The most general second-order covariant interaction vertex reads

$\begin{array}[]{c}\displaystyle\mathbb{T}{}^{(2)}{}_{\mu}=l_{1}[\mathcal{F}{}_{\nu},[\mathcal{F}{}_{\mu},\mathcal{F}{}^{\nu}]]+l_{2}[\mathcal{F}{}_{\nu},[\mathcal{G}{}_{\mu},\mathcal{F}{}^{\nu}]]+l_{3}[\mathcal{F}{}_{\nu},[\mathcal{F}{}_{\mu},\mathcal{G}{}^{\nu}]]+\\\\[19.91692pt]
\displaystyle+l_{4}[\mathcal{G}{}_{\nu},[\mathcal{F}{}_{\mu},\mathcal{G}{}^{\nu}]]+l_{5}[\mathcal{G}{}_{\nu},[\mathcal{G}{}_{\mu},\mathcal{F}{}^{\nu}]]+l_{6}[\mathcal{G}{}_{\nu},[\mathcal{G}{}_{\mu},\mathcal{G}{}^{\nu}]]\,,\end{array}$
(94)

where $l_{p},p=\overline{1,6}$ are constants. The covariant divergence of
equations of motion reads (only cubic terms are written out)

$\begin{array}[]{c}\displaystyle
D_{\mu}\mathbb{T}{}^{\mu}=\frac{1}{\alpha_{3}}\epsilon_{\mu\nu\rho}[\mathcal{S}^{\mu\nu},\mathbb{T}^{\rho}]+\frac{1}{2}C_{1}(k;l)\varepsilon_{\mu\nu\rho}[\mathcal{F}{}^{\mu},[\mathcal{F}{}^{\nu},\mathcal{G}{}^{\rho}]]+\frac{1}{2}C_{2}(k;l)\varepsilon_{\mu\nu\rho}[\mathcal{G}{}^{\mu},[\mathcal{G}{}^{\nu},\mathcal{F}{}^{\rho}]]+\\\\[14.22636pt]
\displaystyle+\frac{l_{1}}{2}[\mathcal{F}{}_{\mu},[\mathcal{F}{}_{\nu},D^{\mu}\mathcal{F}{}^{\nu}+D^{\nu}\mathcal{F}{}^{\mu}]]+\frac{l_{3}-l_{2}}{2}[\mathcal{G}{}_{\mu},[\mathcal{F}{}_{\nu},D^{\mu}\mathcal{F}{}^{\nu}+D^{\nu}\mathcal{F}{}^{\mu}]]+\frac{2l_{2}-l_{3}}{2}[\mathcal{F}{}_{\mu},[\mathcal{G}{}_{\nu},D^{\mu}\mathcal{F}{}^{\nu}+D^{\nu}\mathcal{F}{}^{\mu}]]+\\\\[14.22636pt]
\displaystyle+\frac{l_{5}}{2}[\mathcal{G}{}_{\mu},[\mathcal{G}{}_{\nu},D^{\mu}\mathcal{F}{}^{\nu}+D^{\nu}\mathcal{F}{}^{\mu}]]+\frac{l_{3}}{2}[\mathcal{F}{}_{\mu},[\mathcal{F}{}_{\nu},D^{\mu}\mathcal{G}{}^{\nu}+D^{\nu}\mathcal{G}{}^{\mu}]]+\frac{l_{5}-l_{4}}{2}[\mathcal{F}{}_{\mu},[\mathcal{G}{}_{\nu},D^{\mu}\mathcal{G}{}^{\nu}+\\\\[14.22636pt]
\displaystyle+D^{\nu}\mathcal{G}{}^{\mu}]]+\frac{2l_{4}-l_{5}}{2}[\mathcal{G}{}_{\mu},[\mathcal{F}{}_{\nu},D^{\mu}\mathcal{G}{}^{\nu}+D^{\nu}\mathcal{G}{}^{\mu}]]+\frac{l_{6}}{2}[\mathcal{G}{}_{\mu},[\mathcal{G}{}_{\nu},D^{\mu}\mathcal{G}{}^{\nu}+D^{\nu}\mathcal{G}{}^{\mu}]]\,.\end{array}$
(A20)

Here, the notation is used,

$C_{1}(k;l)=k_{2}{}^{2}-k_{1}k_{3}-3l_{1}-\frac{\alpha_{1}}{\alpha_{3}}l_{4}-\frac{\alpha_{1}}{\alpha_{3}}l_{5}-\frac{\alpha_{1}}{\alpha_{3}}l_{6}\,,\qquad
C_{2}(k;l)=3l_{2}+l_{3}-\frac{\alpha_{2}}{\alpha_{3}}l_{5}-\frac{\alpha_{1}}{\alpha_{3}}l_{6}\,;$
(A21)
$\begin{array}[]{c}\displaystyle\mathcal{S}^{\mu\nu}=l_{3}[\mathcal{F}^{\mu},[\mathcal{F}{}^{\nu},\cdot]]+l_{4}[\mathcal{G}^{\mu},[\mathcal{F}{}^{\nu},\cdot]]+l_{6}[\mathcal{G}^{\mu},[\mathcal{G}{}^{\nu},\cdot]]-l_{4}[[\mathcal{F}^{\mu},\mathcal{F}{}^{\nu}],\cdot]-\\\\[19.91692pt]
\displaystyle-
l_{5}[[\mathcal{F}^{\mu},\mathcal{G}{}^{\nu}],\cdot]-l_{6}[[\mathcal{G}^{\mu},\mathcal{G}{}^{\nu}],\cdot]\,.\end{array}$
(A22)

The interaction is consistent at the second order of perturbation theory if
the right hand side of this expression vanish modulo free equations (3). The
critical observation is that the expressions of the form
$[\mathcal{X}_{\mu}\,,[\mathcal{Y}{}_{\nu}\,,D^{\mu}\mathcal{Z}^{\nu}+D^{\nu}\mathcal{Z}^{\mu}]]$,
where
$\mathcal{Y}{}_{\mu},\mathcal{Z}{}_{\mu},\mathcal{Y}{}_{\mu}=\mathcal{F}{}_{\mu}$
or $\mathcal{G}{}_{\mu}$ represent on-shell independent combinations of fields
and their derivatives. Once they have to vanish, we conclude

$k_{2}{}^{2}-k_{3}k_{1}=0\,,\qquad l_{p}=0\,,\qquad p=\overline{1,6}\,.$ (A23)

The general solution to these equations has the form (93) with $\beta_{3}=0$.
With account of this fact, the interaction (16) is unique in the class of
covariant couplings without higher derivatives.

## References

  * [1] M. Pavsic, Pais-Uhlenbeck oscillator and negative energies, Int. J. Geom. Meth. Mod. Phys. 13 (2016) 1630015.
  * [2] A.V. Smilga, Classical and quantum dynamics of higher-derivative systems, Int. J. Mod. Phys. A 32 (2017) 1730025.
  * [3] N. Boulanger, F. Buisseret, F. Dierick, O. White, Higher-derivative harmonic oscillators: stability of classical dynamics and adiabatic invariants, Eur. Phys. J. C 79 (2019) 60.
  * [4] T.P. Sotiriou, V. Faraoni, _f(R) theories of gravity_ , Rev. Mod. Phys. 82 (2010) 451-497.
  * [5] A. De Felice, S. Tsujikawa, _f(R) theories_ , Living Rev. Rel. 13 (2010) 3\.
  * [6] E.T. Tomboulis, _Renormalization and unitarity in higher derivative and nonlocal gravity theories_ , Mod. Phys. Lett. A 30 (2015) 1540005.
  * [7] A. Belenchia, M. Letizia, S. Liberati, E. Di Casola, _Higher-order theories of gravity: diagnosis, extraction and reformulation via non-metric extra degrees of freedom-a review_ , Rept. Prog. Phys. 81 (2018) 036001.
  * [8] T. Chen, M. Fasiello, E.A. Lim, A.J. Tolley, _Higher derivative theories with constraints: exorcising Ostrogradskis ghost_ , JCAP 1302 (2013) 042.
  * [9] A. Pais and G .E. Uhlenbeck, _On field theories with non-localized action_ , Phys. Rev. 79 (1950) 145-165.
  * [10] B. Podolsky, _A generalized electrodynamics. Part I - non-quantum_ , Phys. Rev. 62 (1942) 68-71.
  * [11] T. Lee and G. Wick, _Negative Metric and the Unitarity of the S Matrix_ , Nucl. Phys. B 9 (1969) 209-243.
  * [12] T. Lee and G. Wick, _Finite theory of quantum electrodynamics_ , Phys. Rev. D 2 (1970) 1033-1048.
  * [13] S. Deser, R. Jackiw, _Higher derivative Chern-Simons extensions_ , Phys. Lett. B 451 (1999) 73-76.
  * [14] H. Weyl, _Gravitation und Elektrizitat_ , Sitzungsberichte der Koniglich Preussischen Akademie der Wissenschaften zu Berlin (1918) 465 - 480.
  * [15] H. Lu, Yi Pang and C.N. Pope, _Conformal Gravity and Extensions of Critical Gravity_ , Phys. Rev. D 84 (2011) 064001.
  * [16] A.A. Nogueira, C. Palechor, A.F. Ferrari, _Reduction of order and Fadeev-Jackiw formalism in generalized electrodynamics_ , Nucl.Phys. B 939 (2019) 372-390.
  * [17] E.C. Manavella, _Quantum field formalism for the higher-derivative nonrelativistic electrodynamics in $1+1$ dimensions_, Int. J. Mod. Phys. A 34 (2019) 1950050.
  * [18] J. Dai, _Stability and consistent interactions in higher derivative matter field theories_ , Eur. Phys. J. Plus 135 (2020) 555.
  * [19] A.V. Smilga, _Benign vs malicious ghosts in higher-derivative theories_ , Nucl. Phys. B 706 (2005) 598-614.
  * [20] M. V. Ostrogradski, _Memoires sur les equations differentielles relatives au probleme des isoperimetretres_. Mem. Acad. St. Petersburg 6 (1850) 385-517.
  * [21] D. M. Gitman, S. L. Lyakhovich, I. V. Tyutin, _Hamilton formulation of a theory with high derivatives_ , Sov. Phys. J. 26 (1983) 61-66.
  * [22] I. L. Buchbinder, S. L. Lyahovich, _Canonical quantisation and local measure of $R^{2}$ gravity_, Class. Quantum Grav. 4 (1987) 1487-1501.
  * [23] J. Kluson, M. Oksanen, A. Tureanu, _Hamiltonian analysis of curvature-squared gravity with or without conformal invariance_ , Phys. Rev. D 89 (2014) 064043\.
  * [24] Y. Ohkuwa and Y. Ezawa, _On the canonical formalism of $f(R)$-type gravity using Lie derivatives_, Eur. Phys. J. Plus 77 (2015) 130.
  * [25] K. Bolonek, P. Kosinski, _Hamiltonian structures for Pais-Uhlenbeck oscillator_ , Acta Phys. Polon. B 36 (2005) 2115.
  * [26] E.V. Damaskinsky and M.A. Sokolov, _Remarks on quantization of Pais-Uhlenbeck oscillators_ , J. Phys. A: Math. Gen. 39 (2006) 10499.
  * [27] C.M. Bender and P.D. Mannheim, _No-ghost theorem for the fourth-order derivative Pais-Uhlenbeck oscillator model_ , Phys. Rev. Lett. 100, (2008) 110402.
  * [28] A. Mostafazadeh, _A Hamiltonian formulation of the Pais-Uhlenbeck oscillator that yields a stable and unitary quantum system_ , Phys. Lett. A 375 (2010) 93-98.
  * [29] M. Avendaco-Camacho, J.A. Vallejo, Yu. Vorobiev, _A perturbation theory approach to the stability of the Pais-Uhlenbeck oscillator_ , J.Math.Phys. 58 (2017) 093501.
  * [30] G.W. Gibbons, C.N. Pope, and S. Solodukhin, _Higher Derivative Scalar Quantum Field Theory in Curved Spacetime_ , Phys. Rev. D 100 (2019) 105008
  * [31] A. Salvio, _Metastability in Quadratic Gravity_ , Phys.Rev.D 99 (2019) 103507.
  * [32] D.S. Kaparulin, S.L. Lyakhovich, A. A. Sharapov, _Classical and quantum stability of higher-derivative dynamics_ , Eur. Phys. J. C 74 (2014) 3072.
  * [33] D.S. Kaparulin, I.Yu. Karataeva, S.L. Lyakhovich, _Higher derivative extensions of 3d Chern-Simons models: conservation laws and stability_ , Eur. Phys. J. C 75 (2015) 552.
  * [34] P.O. Kazinski, S.L. Lyakhovich, A.A. Sharapov, _Lagrange structure and quantization_ , JHEP 0507 (2005) 076.
  * [35] D.S. Kaparulin, S.L. Lyakhovich, A.A. Sharapov, _Rigid symmetries and conservation laws in non-Lagrangian field theory_ , J. Math. Phys. 51 (2010) 082902.
  * [36] D.S. Kaparulin, S.L. Lyakhovich, A.A. Sharapov, _BRST analysis of general mechanical systems_ , J. Geom. Phys. 74 (2013) 164-184.
  * [37] D.S. Kaparulin, S.L. Lyakhovich, A.A. Sharapov, _Stable interactions via proper deformations_ , J. Phys. A: Math. and Theor. 49 (2016) 155204.
  * [38] V.A. Abakumova, D.S. Kaparulin, S.L. Lyakhovich, _Stable interactions in higher derivative field theories of derived type_ , Phys. Rev. D 99 (2019) 045020.
  * [39] S.C. Sararu and M. Tinca, _Udristioiu Quantization of the higher derivative Maxwell-Chern-Simons-Proca model based on BFT method_ , Mod. Phys. Lett. A31(2016) 1650205.
  * [40] M. Ghasemkhani, and R. Bufalo, _Noncommutative Maxwell-Chern-Simons theory: One-loop dispersion relation analysis_ , Phys. Rev. D93 (2016) 085021.
  * [41] R. Avila, J.R. Nascimento, A.Yu. Petrov, C.M. Reyes, and M. Schreck, _Causality, unitarity, and indefinite metric in Maxwell-Chern-Simons extensions_ , Phys. Rev. D101 (2020) 055011\.
  * [42] J. Dai, _BRST deformations and stability in the higher derivative Chern-Simons gauge theory_ , Class. Quantum Grav. 37 (2020) 245011\.
  * [43] V.A. Abakumova, D.S. Kaparulin, S.L. Lyakhovich, _Multi-Hamiltonian formulations and stability of higher-derivative extensions of $3d$ Chern-Simons_, Eur. Phys. J. C78 (2018) 115.
  * [44] D.S. Kaparulin, I.Yu. Karataeva, S.L. Lyakhovich, _Third order extensions of $3d$ Chern-Simons interacting to gravity: Hamiltonian formalism and stability_, Nucl. Phys. B 934 (2018) 634-652.
  * [45] M. Henneaux, Consistent interactions between gauge fields: the Cohomological approach, Contemp. Math. 219 (1998) 93.
  * [46] G. Barnich, F. Brandt, M. Henneaux, _Local BRST cohomology in gauge theories_ , Phys. Rep. 338 (2000) 439.
  * [47] C. Bizdadea, _On the cohomological derivation of topological Yang-Mills theory_ , EPL 49 (2000) 123-129.
  * [48] D.S. Kaparulin, S.L. Lyakhovich, and O.D. Nosyrev, _Resonance and stability of higher derivative theories of derived type_ , Phys. Rev. D 101 (2020) 125004.
  * [49] P. K. Townsend, K. Pilch, and P. van Nieuwenhuizen, _Self-duality in odd dimensions_ , Phys. Lett. B 136 (1984) 38-42 .
  * [50] B. Binegar, _Relativistic field theories in three dimensions_ , J. Math. Phys. 23 (1982) 1511-1517.
  * [51] D.R. Grigore, _The projective unitary irreducible representations of the Poincare group in (1+2)-dimensions_ , J. Math. Phys. 34 (1993) 4172-4189.
  * [52] D.R. Grigore, _Free fields for any spin in (1+2)-dimensions_ , J. Math. Phys. 35 (1994) 6304-6331.
  * [53] D.S. Kaparulin, S.L. Lyakhovich, A.A. Sharapov, _Consistent interactions and involution_ , JHEP 01 (2013) 097.
  * [54] L.D. Faddeev, A.A. Slavnov, Gauge Fields: An Introduction To Quantum Theory, Frontiers in Physics, Vol. 50 (Benjamin/Cummings, Advanced Book Program, 1980).
  * [55] M. Henneaux and C. Teitelboim, Quantization of Gauge Systems (Princeton U.P., NJ, 1992).